# Mixing-demixing transition and void formation in quasi-2D binary mixtures on
a sphere

D. Truzzolillo<EMAIL_ADDRESS>Laboratoire Charles
Coulomb (L2C), UMR 5221 CNRS-Universitè de Montpellier, F-34095 Montpellier,
France

###### Abstract

Motivated by observations of heterogeneous domain structure on the surface of
cells and vesicles and by domain formation due to the adsorption of complex
molecules onto composite membranes, we consider a minimal quasi 2D-model to
describe the structure of binary mixtures on the surface of a spherical
particle. We study the effect of miscibility and adsorbing particle (AP)
addition on the mixture structure. We define a new scalar quantity, the
geodesic mixing parameter $\Xi$, through which we detail the effect of
miscibility and the role of preferential affinity of APs with one of the two
components of the mixture, distinguishing unambiguously between mixing and
demixing solely induced by APs. Finally, by inspecting the distributions of
void sizes, we show how void formation is ruled by miscibility and AP-mixture
interactions, which control the transition from exponentially-tailed to fat-
tailed distributions.111The following article has been accepted by _The
Journal of Chemical Physics_. After it is published, it will be found at
https://aip.scitation.org/journal/jcp

††preprint: AIP/123-QED

##  I Introduction

The investigation of the equilibrium properties of particles whose motion is
confined to curved surfaces is surely challenging, but it is essential to
achieve a satisfactory understanding of processes such as the lateral
diffusion of proteins in membranes Malchus and Weiss (2010); Javanainen  _et
al._ (2013); Metzler, Jeon, and Cherstvy (2016), the controlled synthesis of
colloidal aggregatesSavarala  _et al._ (2011), and demixing in lipid
bilayersAnderson (2002); Binder, Barragan, and Menger (2003); Veatch  _et al._
(2008). Confining particles to the surface of a sphere introduces several
constraints to their placement that are absent in Euclidean (flat) space. Both
the curvature and the topology of the sphere play a role in the way particles
can be arranged on it—most prominently, no regular lattice can be fit onto its
surface, as the topology of the sphere requires the presence of 12 pentagonal
disclinations Bowick and Giomi (2009), and introduces structural defects
characterizing glassy fluids at large area fractions Singh, Sood, and
Ganapathy (2020). The compactness of the sphere also requires a careful
treatment: instead of the straightforward thermodynamic limit, configurations
of finite (small) numbers of particles become more relevant, and both the
number of particles and the size of the sphere need to be considered as two
independent parameters Post and Glandt (1986). Among all possible spherical
systems, multicomponent mixtures are of particular importance since they
represent the class of systems including biological membranes in which the
self-organisation of lipids into different domains is critical for many cell
properties such as recognition, signalling, or endocytosis and exocytosis
Salaün, James, and Chamberlain (2004). Demixing and segregation in bulk
systems can proceed via a number of distinct physical mechanisms (such as
spinodal decomposition and nucleation), and it is generally well understood.
However, the impact of spatial confinement and the interaction with external
agents on the structure of multicomponent mixtures are less clear. Confinement
can arise from the presence of external fields, produced, for example, by
substrates or random obstacles, but they can also be imposed by the geometry
of the embedding space. The latter type of confinement is particularly
relevant to biological cells, for which the mobile fluid particles
constituting the cell membrane are constrained to lie on the surface of a
quasi-spherical body. This has boosted the interest for a deep understanding
of multiphase fluids on smooth manifolds, especially on spheres and spheroids,
and so in the past 2 decades many numerical studies have been focused on the
particular problem of the phase separation of binary systems on surfaces and
its dependence on the surface curvature. Tang  _et al._ (2005); Greer,
Bertozzi, and Sapiro (2006); Marenduzzo and Orlandini (2013); Barreira,
Elliott, and Madzvamuse (2011) Biological membranes are indeed described as a
quasi two- dimensional fluids Saffman and Delbruck (1975); Camley and Brown
(2011) in which the different constituents flow within the membrane manifold,
being still partially free to move orthogonally to the membrane tangent plane
and coupled to a less viscous three-dimensional surrounding medium. They
exhibit very often stable domains, whose spatial distribution is important for
many cell propertiesLenne and Nicolas (2009); Lingwood and Simons (2010). One
common view is that the domains on the surface of a cell are a consequence of
an arrested or incomplete phase separation, however it remains to be
established whether the observed states are permanent or metastable in
character. If the heterogeneous domain structure on the cell surface is an
equilibrium state, then some stabilizing mechanism is required; the line
tension present at interfaces between domains would make inhomogeneous phases
energetically unfavorable when compared to a fully phase-separated system.

Many works have suggested that biological membranes are close to a demixing
line in the temperature-composition plane. When the lipid mixture is close to
the miscibility critical pointOnuki (2002); Honerkamp-Smith, Veatch, and
Keller (2009) (the so-called weak-segregation limit) the temperature $T$ is
close to a critical temperature $T_{c}$. In this limit, density fluctuations
give rise to wide interfaces between fluctuating lipid domains, which do not
have a precise shapeDestainville, Manghi, and Cornet (2018).

As a matter of fact, if a membrane is composed of two or more lipids (as most
biomembranes are), the segregation of one of the lipid components in a set of
domains is often encountered, as for example rafts rich in cholesterol and
sphingolipids that have been implicated recently in signal transduction and
membrane trafficking pathwaysAnderson (2002); Binder, Barragan, and Menger
(2003). Veatch and coworkersVeatch  _et al._ (2008) suggested that lateral
heterogeneities in living cells at physiological conditions correspond to
critical fluctuations and that perturbations that alter the phase boundary
could have a large effect on the size, composition, and lifetime of
fluctuations at physiological temperatures. Also, lipid demixing due to the
existence of a miscibility critical point in red blood cell (RBC) membranes
has been highlighted via Langmuir trough and epifluorescence experiments by
Keller and coworkers Keller  _et al._ (1998), who showed that monolayers of
lipids extracted from RBC membranes at room temperature segregate and form
domains above a threshold surface pressure at the air-water interface.

In addition to single membrane properties, the interaction potential between a
pair of cells is also strongly influenced by the distribution and size of the
domains covering its surface. In this sense, cells may be regarded as a
naturally occurring type of “patchy particles,” which are particles with
distinct surface sites generating anisotropic interparticle interactions.
While synthetically fabricated patchy particles have attractive patches
strategically arranged on their surface Bianchi, Blaak, and Likos (2011), the
domains covering a cell emerge as a result of self-organization. Recent
developments in the controlled fabrication of patchy particles have raised
hopes that materials with desired properties may be tailored by prescribing
the number and geometrical arrangement of the patches Bianchi, Blaak, and
Likos (2011). In this framework, confined binary mixtures of active and
passive particles on a sphere have also been studiedAi, Zhou, and Zhang (2020)
very recently, since particle activity coupled to the presence of intrinsic
surface curvature frustrates local order, giving rise to novel phenomena
Shankar, Bowick, and Marchetti (2017); Keber  _et al._ (2014); Sanchez  _et
al._ (2012), such as curvature-induced defects and unforced flow.

On top of the spontaneous formation of domains and defects on spheres, both
loss of miscibility and domain destabilization in spherically confined fluids
can be also caused by the addition of adsorbing agents, like multivalent ions,
polyelectrolytes and polyampholytes whose adsorption-desorption mechanism
represents one the key problem tackled in soft matter physicsCaetano  _et al._
(2020, 2017). Induced-demixing has been widely investigated in the past.
Domain and raft formation in membranes has been found in the presence of
bivalent ions Jacobson and Papahadjopoulos (1975), RNA and DNA strands
Michanek  _et al._ (2010), charged proteins Mbamala, Ben-Shaul, and May
(2005), and oppositely charged polyelectrolytes Macdonald  _et al._ (1998). By
contrast, induced mixing has been documented much less and only very recently
it has been shown that the addition of functionalized nanoparticles Canepa
_et al._ (2020) or thermoplastic polymers can indeed give rise to a
suppression of lipid segregation Bochicchio  _et al._ (2017). The general
mechanisms dictating a mixing induced by adsorbing molecules, polymers or
colloids and the thermodynamic conditions under which this occurs, have not
been explored neither understood in depth. Finally, domain formation and lipid
miscibility plays a crucial role for the membrane permeability. Passive
transport through biomembranes and bilayer lipid membranes is in many cases
controlled by existing small holes (or voids) within them Kotyk and Janáček
(1977); Kashchiev and Exerowa (1983). There are several experimental evidences
that proton, water and potentially drug permeability of mixed lipid membranes
is enhanced by raft formationLaRocca  _et al._ (2013); Ghysels  _et al._
(2019); Gensure, Zeidel, and Hill (2006); Sciolla  _et al._ (2021), since
inter-domain interfaces enhances ions and small molecules intrusion, and that
membranes can be even destabilized by lipid substitutionLaRocca  _et al._
(2013) or multivalent ionsHa (2001). In addition to biological membranes, hole
formation have been used to describe the rupture of bilayer (Newtonian) black
foam films Derjaguin and Prokhorov ; Kashchiev and Exerowa (1980), the
evaporation of liquids through adsorbed monolayers Dickinson (1978), as well
as the thermodynamic equilibrium of monolayers on liquid surfacesStoeckly
(1977).

This said, despite an undisputed relevance of liquid mixtures lying on spheres
and void formation in them, a systematic investigation through simple models
on the role played by miscibility and its coupling to bulk particle adsorption
in determining the mixture structure and in triggering void nucleation is
still lacking.

In this work we introduce a minimal model to study the structure across a
consolute pointWilding, Schmid, and Nielaba (1998) of a symmetric Lennard-
Jones binary mixture confined on a spherical surface by a harmonic potential,
leaving particles free to oscillate along the radial direction. Such a
quasi-2D model is employed first to characterize the structure of bare liquid
mixtures at fixed temperature (here $T^{*}=0.5$) and varying miscibility
parameter ($\alpha_{\scriptscriptstyle AB}$) ruling the demixing of the two
fluid components (A and B). We compute both the total and the partial geodesic
pair correlation functions, and we introduce a geodesic mixing parameter
($\Xi$) through which we describe both the spontaneous demixing and the effect
of adsorbing bulk particle (AP) addition on the mixture structure. We focused
on the role of the AP preferential affinity with one of the two mixture
components and we give a simple structural criterion to distinguish between AP
induced mixing and demixing, showing that the former occurs only for nearly
demixed systems with large spatial fluctuations of domain boundaries. For both
AP-free mixtures and AP-decorated ones, we further investigate void size
distributions. We show that in the presence of spatially fluctuating domains
and inter-domain interfaces, large void formation is not favoured, and that AP
adsorption can both enhance and disfavour the onset of large voids, depending
on the miscibility of the two fluid constituents.

##  II Model

The binary mixture studied in this paper consists of fluids ($A$ and $B$) made
of spherical particles of the same size, $\sigma_{AA}=\sigma_{BB}$, and at
concentrations of 50% each. Particles of the same type interact through a
truncated 12–6 Lennard-Jones (LJ) potentialDíaz-Herrera, Moreno-Razo, and
Ramírez-Santiago (2004):

$\begin{cases}U_{LJ}(r_{ij})=4\varepsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right],&\mbox{if
}r<3$$\sigma_{\scriptscriptstyle AA}$$\\\ 0,&\mbox{if
}r>3$$\sigma_{\scriptscriptstyle AA}$$\end{cases}$ (1)

with a mixing rule

$\sigma_{\scriptscriptstyle AB}=\frac{1}{2}(\sigma_{\scriptscriptstyle
AA}+\sigma_{\scriptscriptstyle BB}),\quad\varepsilon_{\scriptscriptstyle
AB}=\alpha_{\scriptscriptstyle AB}\varepsilon_{\scriptscriptstyle AA},$ (2)

where $\varepsilon_{\scriptscriptstyle AA}=\varepsilon_{\scriptscriptstyle
BB}$, and $\alpha_{\scriptscriptstyle AB}$ is the parameter controlling the
miscibility of the two fluids. Notice that when, $\alpha_{\scriptscriptstyle
AB}=0$, we obtain two independent LJ fluids while in the opposite case,
$\alpha_{\scriptscriptstyle AB}=1$, the system reduces to a single LJ fluid.
By choosing, $0<\alpha_{\scriptscriptstyle AB}<1$, the attractive part of the
interaction between dissimilar particles becomes weaker than that of the $AA$
and $BB$ interactions, favouring demixingDíaz-Herrera, Moreno-Razo, and
Ramírez-Santiago (2004). Each particle vibrates around the ideal surface of a
sphere of radius $R_{s}$ since they are further confined via a parabolic
potential:

$U_{c}(r_{i})=\frac{1}{2}\kappa_{c}(r_{i}-R_{s})^{2},$ (3)

where $\frac{\kappa_{c}\sigma_{\scriptscriptstyle AA}^{2}}{2k_{B}T}=10$ and
$r_{i}$ is the distance between particle $i$ and the center of the confining
sphere. To check the influence of the system size on demixing, we simulated
three systems characterized by different particle numbers
($N_{p}=250,500,1000$) and the same area fraction occupied by the particles
$\phi=\frac{N_{p}\sigma_{\scriptscriptstyle AA}^{2}}{16R_{s}^{2}}=0.5$,
corresponding to a reduced density $\rho^{*}=\frac{4\phi}{\pi}=0.63$. The
choice of the system size has been made to ensure that the particle-to-sphere
size ratio $\sigma_{\scriptscriptstyle AA}/2R_{s}$, ranging here from 0.044 to
0.089, was comparable to the bilayer-to-vesicle size ratios encountered in
small lipid vesicles Sciolla  _et al._ (2021), the latter being $\approx
0.05$. The attractive interaction between the confined particles ($\leq
2k_{B}T$, see section III) and the energy scale
($\frac{\kappa_{c}\sigma_{\scriptscriptstyle AA}^{2}}{2}$) of the confining
potential have been chosen to reproduce the same order of magnitude
respectively of the in-plane lipid-lipid interactions (few $k_{B}T$
unitsKomura  _et al._ (2004) for homogeneous systems and possibly less than
$k_{B}T$ for dissimilar lipidsAlmeida, Pokorny, and Hinderliter (2005)) and
the bending energy of membranes per lipid (from few tens to hundreds of
$k_{B}T$ per molecule Peliti (1991)). The first allows molecules to diffuse
and structure itself as a fluid and its variation determines whether demixing
occurs at room temperature if more than one lipid type is present. The second
gives rise to membrane stability. The magnitude of the confining potential has
been further selected so as to preserve equipartition of kinetic energy during
all the simulation runs and avoid flying ice-cube artifacts.

Free unconfined APs interacting with the mixture via the potential (1) with LJ
diameter $\sigma_{\scriptscriptstyle PP}=\sigma_{\scriptscriptstyle AA}$ have
been further added to mixtures composed by $N_{p}=500$ particles in order to
study the effect of physical adsorption on the structure of the confined AB-
system. The interaction between APs and the mixtures is set by equation 1,
where miscibility between APs and A-type particles has been fixed by setting
their miscibility parameter to $\alpha_{\scriptscriptstyle AP}=1$, while that
between APs and B-type particles has been tuned by varying their mutual
miscibility parameter (hereinafter $\alpha_{\scriptscriptstyle BP}$) in the
range $1-0.025$. In this way the interaction between the adsorbing particles
and the confined mixtures is such as to be preferential (for
$\alpha_{\scriptscriptstyle BP}<\alpha_{\scriptscriptstyle AP}$) or, let’s
say, symmetrical (for $\alpha_{\scriptscriptstyle
BP}=\alpha_{\scriptscriptstyle AP}$) while the total adsorption energy
decreases when $\alpha_{\scriptscriptstyle BP}$ is lowered. In such a way the
simulated systems mimic the weak physical adsorption of small molecules whose
adsorption energies do not exceed few $k_{B}T$ units. The surface of the
confining sphere has been further made impenetrable to bulk APs via the purely
repulsive potential:

$U_{LJ}(r_{\scriptscriptstyle PS})=4\varepsilon_{\scriptscriptstyle
PS}\left[\left(\frac{2R_{s}+\sigma_{\scriptscriptstyle
PP}}{r_{\scriptscriptstyle PS}}\right)^{12}\right],$ (4)

where $\varepsilon_{\scriptscriptstyle PS}/(k_{B}T)=2\cdot 10^{-4}$ and
$r_{PS}$ is the center-to-center distance between an AP and the sphere
confining the AB mixture.

##  III Simulation details

We have carried out extensive MD simulations to investigate the structural
properties of this model binary mixture as a function of
$\alpha_{\scriptscriptstyle AB}$ and $\alpha_{\scriptscriptstyle BP}$. The
equations of motion were integrated using a velocity-Verlet algorithm with
reduced time step $\Delta t^{*}=\frac{\Delta t}{\sigma_{\scriptscriptstyle
AA}}\sqrt{\frac{\varepsilon_{\scriptscriptstyle AA}}{m_{A}}}=5\cdot 10^{-6}$.
$m_{A}$ is the particle mass for species $A$ and $m_{A}=m_{B}=m_{P}$. We fixed
in all simulations the reduced temperature at
$T^{*}=\frac{k_{B}T}{\varepsilon_{AA}}=0.5$, where 2D Lennard-Jones fluids are
in a condensed liquid state Barker, Henderson, and Abraham (1981); Holian and
Grady (1988). For the chosen set of miscibility parameters
($\alpha_{\scriptscriptstyle AB}$,$\alpha_{\scriptscriptstyle BP}$) we
equilibrated the systems using the Berendsen’s thermostat with a large
coupling constant $\tau_{T}=2\cdot 10^{3}\Delta t$ to obtain a stable
trajectory in equilibrium. Equilibration runs have been performed for $2\cdot
10^{8}$ time steps and production runs have been carried out in the NVE
ensemble for $5.2$ $\cdot 10^{7}$ time steps during which total and kinetic
energy drifts were absent. To exclude the presence of important ballistic
collective motions of the confined particles due to an eventual rigid rotation
of mixtures on the sphere, and energy transfer from high to low frequency
modes, we have further performed long NVT simulations for $N_{p}=1000$ and
$N_{p}=500$ at $\alpha_{\scriptscriptstyle AB}=1$ to compute the mean angular
square displacements: a pure diffusive behavior has been observed (see
Appendix A). To minimize correlations between measurements we calculated
structural quantities every $\Delta t_{s}=5\cdot 10^{{\color[rgb]{0,0,0}5}}$
time steps (see Appendix D). Simulations including $N_{\scriptscriptstyle
AP}=1000$ free bulk particles have been carried out only for mixtures with
$N_{p}=500$. The choice of such an excess of bulk particles has been inspired
by the experimental system recently investigated by the author and coworkers
Sciolla  _et al._ (2021) in which lipid phase separation in mixed liposomes
has been observed in excess of Isoniazide, one of the primary drugs used in
the tuberculosis treatment. An investigation on the effect of bulk particle
concentration, temperature and mixture composition are currently in progress
and will be the subject of a future publication. Bare mixtures confined on the
sphere have been simulated in a cubic box of size $L=10R_{s}$ with no extra
periodic boundary conditions, with the sphere being at the center of the box.
We underline here that the main characteristic of such systems is to be finite
but without boundaries, with periodic boundary conditions being not an
artifact of the simulations, but an essential feature of it. Hence, in this
case there is no need for further assumptions or corrections: the complete
system is simulated. Conversely, we carried out simulations including bulk APs
in a box of size $L=4.2R_{s}$, correspondent to a reduced density for the bulk
particles $\rho^{*}=\frac{N_{\scriptscriptstyle
AP}\sigma_{CC}^{3}}{(L^{3}-4/3\pi R_{s}^{3})}=0.029$, with cubic periodic
boundary conditions applied only to free APs.

##  IV Bare mixtures

The structure of bare mixtures has been studied by inspecting the total $g(s)$
and the partial $g_{m}(s)$ pair correlation functions defined as:

$g(s)=\frac{1}{2N_{p}^{2}sin(s)}\frac{dn(s)}{ds}$ (5)
$g_{m}(s)=\frac{1}{2}(g_{\scriptscriptstyle AB}(s)+g_{\scriptscriptstyle
BA}(s))=\\\
\frac{1}{2}\left[\frac{1}{2N_{A}^{2}sin(s)}\frac{dn_{\scriptscriptstyle
AB}(s)}{ds}+\right.\\\
\left.\frac{1}{2N_{B}^{2}sin(s)}\frac{dn_{\scriptscriptstyle
BA}(s)}{ds}\right]$ (6)

where $s\in[0,\pi]$ is the normalized geodesic distance on the unit sphere,
$dn(s)$ is the average number of particles at geodesic distance comprised
between $s$ and $s+ds$ from one test particle and $dn_{\scriptscriptstyle
AB}(s)$ ($dn_{\scriptscriptstyle BA}(s)$) is the average number of particles
of type B (A) at geodesic distance comprised between $s$ and $s+ds$ from a
test particle of type A (B). For homogeneous mixtures
($\varepsilon_{\scriptscriptstyle AA}=\varepsilon_{\scriptscriptstyle
BB}=\varepsilon_{\scriptscriptstyle AB}$) $g_{m}(s)=g(s)$, while the largest
difference between the two functions is attained when the two fluids are fully
separated on the sphere in a Janus-like configuration.
Figure 1 shows snapshots of equilibrated samples of mixtures for $N_{p}=1000$
with different mixing parameters $\alpha_{\scriptscriptstyle AB}$ ranging from
1 down to 0.05. A continuous transition from fully miscible to Janus-like
patterns characterizes the demixing transition with an intermediate regime
$0.5\leq\alpha_{\scriptscriptstyle AB}\leq 0.8$ where large spatial
interfacial fluctuations and domains emerge.

Figure 1: Snapshots of binary 50:50 AB-mixtures ($N_{p}=1000$) confined to
move on a spherical surface with radius of curvature
$R_{s}=11.18\sigma_{\scriptscriptstyle AA}$ and different miscibility
parameters $0.05\leq\alpha_{\scriptscriptstyle AB}\leq 1$, as indicated in the
figure. All snapshots refer to equilibrium states.

By computing both $g(s)$ (Figure 2) and $g_{m}(s)$ (Figure 3) we follow both
the structure of the whole fluid and the demixing of the two components on the
sphere. Though demixing does not affect remarkably the local structure, we
unambiguously observe a weak melting of the mixture induced by the formation
of large domains, as signalled by the global minimum of the height of the
first peak of $g(s)$ as a function of $\alpha_{\scriptscriptstyle AB}$. The
decrease of the first peak reveals that the average local density of the fluid
in the immediate vicinity of any particle decreases with respect to the
homogeneous system when fluctuating domains coexist (inset of Figure 2), since
a complete structural demixing is not achieved, the number of weak AB bonds is
still large and particles are on average less bound. On the contrary, when the
Janus-like configuration is attained and small spatial fluctuations
characterize the interface between the A-rich and the B-rich phases, the
number of AB contacts decreases and the average particle-particle distance in
the fluid goes back to about the one obtained for $\alpha_{\scriptscriptstyle
AB}=1$.

Figure 2: Total pair correlation functions $g(s)$ for $N_{p}=1000$ and
selected values of the miscibility parameter
$0.2\leq\alpha_{\scriptscriptstyle AB}\leq 1$. Inset: height of the first peak
$g(s_{max})$ for the whole set of miscibility parameters investigated in this
work.

The partial correlation functions $g_{m}(s)$ (Figure 3) follow the demixing of
the two components: when particles are fully mixed $g_{m}(s)=g(s)$, while when
$\alpha_{\scriptscriptstyle AB}\ll 1$ and full demixing occurs, the
probability to find an AB contact drastically decreases (small $s$ region) and
it increases progressively until it reaches its maximum at angular distance
$s=\pi$. This is accompanied by the disappearence of the correlation peaks of
$g_{m}(s)$ when $\alpha_{\scriptscriptstyle AB}$ decreases, as shown in the
inset of Figure 3 for representative values of $\alpha_{\scriptscriptstyle
AB}$.

Figure 3: Partial pair correlation functions $g_{m}(s)$ for $N_{p}=1000$ and
selected values of the miscibility parameter
$0.1\leq\alpha_{\scriptscriptstyle AB}\leq 1$. Inset: Same data of the main
panel for $0\leq s\leq 0.8$.

Since confining spherical surfaces are a closed manifold with uniform
curvature, we can describe in a nutshell the demixing by defining a geodesic
mixing parameter as the square distance between $g(s)$ and $g_{m}(s)$ on
$[0,\pi]$:

$\Xi=\|g(s)-g_{m}(s)\|^{2}=\int_{0}^{\pi}[g(s)-g_{m}(s)]^{2}ds.$ (7)

$\Xi$ is zero when the two components of the mixtures are fully mixed and
larger than zero to an extent depending on the structural demixing. In Figure
4 we show $\Xi(\alpha_{\scriptscriptstyle AB})$ obtained for the three
simulated system sizes ($N_{p}=250,500,1000$). $\Xi(\alpha_{\scriptscriptstyle
AB})$ smoothly passes from very low values compatible with zero (given the
simulation noise) for $\alpha_{\scriptscriptstyle AB}=1$ to larger values,
that depends weakly on $\alpha_{\scriptscriptstyle AB}$ when the latter is
lower than $\approx 0.2$ for each system size, pointing out that radius of
curvature $R_{s}$ does not impact remarkably the structural demixing
transition apart from a very small shift of the critical miscibility parameter
$\alpha_{\scriptscriptstyle AB}^{c}$ defined as the point where
$\frac{d^{2}\Xi}{d\alpha_{\scriptscriptstyle AB}^{2}}=0$, that we extract as
detailed in Appendix B. Actually, we expect that compositional fluctuations
start to play a role for very small radius of curvatures impacting more the
angular particle distribution, favouring mixing.

Figure 4: Geodesic mixing parameter $\Xi$ as a function of the miscibility
parameters $\alpha_{\scriptscriptstyle AB}$ for the three different system
sizes $N_{p}=250,500,1000$, and fixed area fraction $\phi=0.5$. Dashed line
are fits obtained via equation 15 (Appendix B).

To gain a complementary insight on the structure of the mixtures and their
stability against void formation we computed the normalized large-void size
distribution $P(A_{h}/\sigma_{\scriptscriptstyle AA}^{2})$ defined as

$P(A_{h}/\sigma_{\scriptscriptstyle
AA}^{2})=\frac{N(A_{h}/\sigma_{\scriptscriptstyle
AA}^{2})}{\sum\limits_{h=1}^{\infty}N(A_{h}/\sigma_{\scriptscriptstyle
AA}^{2})},\quad\mbox{with\hskip
5.69046pt}A_{h}\geq\frac{\sigma_{\scriptscriptstyle AA}^{2}}{4}$ (8)

where $A_{h}$ is the area of a void with size $h\sigma_{\scriptscriptstyle
AA}^{2}/4$ and $N(A_{h}/\sigma_{\scriptscriptstyle AA}^{2})$ is the number of
voids with normalized size $A_{h}/\sigma_{\scriptscriptstyle AA}^{2}$. We
therefore restrict the search of voids to those having area larger than a
particle radius square, since we are not interested to the distribution of
small interstitial spaces between particles but rather to those voids that are
compatible with inclusion of adsorbing particles or molecules in real systems.
Indeed the minimal void-to-particle size ratio taken into consideration in
this work ($A_{h}/\sigma_{\scriptscriptstyle AA}^{2}=1/4$) is of the same
order of magnitude of that characterizing water molecules and lipid cross
section in membrane lipidsIsraelachvili and Mitchell (1975). To compute each
void area $A_{h}$ we applied a spherical grid to the confining sphere, where
$N_{g}$ grid points are almost-equidistant and have been generated following
the scheme proposed inDeserno , with a lattice spacing
$a\simeq\sigma_{\scriptscriptstyle AA}/2$. We consider as a void a subset of
$n_{g}$ grid points such that: i) each point of the subset is far from any
particles by a geodesic distance larger than $3/4\sigma_{\scriptscriptstyle
AA}$; ii) a void has $n_{g}>1$ if each point of the subset have at least one
grid point distant one lattice spacing and belonging to the same subset
(cluster of points). The area of the void is then well approximated by
$A_{h}=n_{g}\sigma_{\scriptscriptstyle AA}^{2}/4$, with a minimum measurable
value equal to $\sigma_{\scriptscriptstyle AA}^{2}/4$. In figure 5 we report
the normalized distribution of void sizes for selected values of
$\alpha_{\scriptscriptstyle AB}$. The reentrant behavior of the distribution
tails goes along with the already mentioned fluidization of the system (Figure
2) when jagged domain boundaries characterize the mixtures. In this regime,
since $AB$ contacts are maximized and particles are on average less bound, the
formation of large voids is hampered. To describe exhaustively this behavior
we have fitted the void distributions with a phenomenological function of the
following type:

$P(A_{h}/\sigma_{\scriptscriptstyle
AA}^{2})=P_{0}+\frac{a}{1+\left(b\frac{A_{h}}{\sigma_{\scriptscriptstyle
AA}^{2}}\right)^{c}}\cdot e^{-A_{h}/A_{c}}$ (9)

where $A_{c}$ is the exponential characteristic cutoff area of the
distribution, $P_{0}$ is an offset set to zero for distribution showing
exponential tails and larger than zero when large voids are created by bulk
particle adsorption (see section V). Finally $a$, $b$ and $c$ characterize the
power law regime observed for small void sizes. The choice of the fitting
function has not been made by chance. Boutreux and De Gennes Boutreux and de
Geennes (1997) first postulated exponential decays for void size distributions
in granular matter and based on this, Caglioti and coworkers Caglioti  _et
al._ (1999) related quantitatively the characteristic dynamical properties for
the same systems to quasistatic quantities, e.g. free-volume and
configurational entropy. In our mixtures, however thermal (and stress)
fluctuations are important and the void formation should also be related to
the probability of spontaneous rupture events of a given size. Rupture events
often show power-law size distributions Amitrano (2012) with an exponential
cut-off due to the finite size of the system, while shifted gamma functions
with exponential tails have proven to be appropriate to describe the
statistics of void sizes in packed protein cores and jammed packings of amino-
acid-shaped particles Treado  _et al._ (2019) and attractive emulsion
dropletsJorjadze  _et al._ (2011). These Boltzmann-type distribution tails are
those maximizing configurational entropy in packed systems Jorjadze  _et al._
(2011). Our results are in line with these previous findings and point to void
distribution tails with two distinct regimes: a power-law decay followed by an
exponential cut-off. Finally, the deviation from a simple power-law at very
low void areas ($\lim_{A_{h}\to 0}P(A_{h}/\sigma_{\scriptscriptstyle
AA}^{2})=P_{0}+a$) reflects the fact that the void size is limited by the
interstitial space between LJ sphere, so that the distribution must converge
to 0 as $A_{h}\rightarrow 0$Treado  _et al._ (2019). Such a deviation is found
systematically for all the sets of parameters employed in our simulations,
hence the need to use equation 9 to fit satisfactorily all our data. A finite
non-zero value for $P_{0}$ further points to a non-exponential (fat) tail of
the distribution. This will be the case of some of the mixtures in presence of
APs. The inset of Figure 5-A shows $A_{c}/\sigma_{\scriptscriptstyle AA}^{2}$
together with the probability to have unit voids
($A_{h}/\sigma_{\scriptscriptstyle AA}^{2}=1/4$). The two quantities show two
non-monotonic trends with opposite concavity, corroborating the scenario in
which demixing on small confining spheres is structurally a smooth process
with the formation of small (large) voids being enhanced (suppressed) in the
presence of fluctuating domains. Similarly to what we observed for the cut-off
areas, the total (large-)void size
$A_{tot}=\sum\limits_{h=1}^{\infty}N(A_{h}/\sigma_{AA}^{2})$shows a minimum
for mixtures with large fluctuating domain interfaces (Figure 5-B). This
further suggests that the integrity of spherical mixtures is stabilized by
compositional fluctuations.

Figure 5: Panel A: Void size distribution $P(A_{h}/\sigma_{AA}^{2})$ for
different miscibility parameter $\alpha_{\scriptscriptstyle AB}$ as indicated
in the panel. The inset shows the cut-of area $A_{c}$ (black sqaures) defined
in equation 9 and the fraction the smallest voids $P(1/4)$ (blue circles) as a
function of $\alpha_{\scriptscriptstyle AB}$. Panel B: Total void area in
function of $\alpha_{\scriptscriptstyle AB}$. The area takes into account only
voids with size $A_{h}\geq\sigma_{\scriptscriptstyle AA}^{2}/4$.

All in all, the driving forces dictating large void formation are two: the
ability of particles to form long-lived clusters, that is enhanced by
increasing the total cohesive energy (here increasing
$\alpha_{\scriptscriptstyle AB}$) and dewetting, occurring for fully demixed
fluids ($\alpha_{\scriptscriptstyle AB}\ll 1$) where void formation is favored
in the inter-domain regions.

##  V Mixtures with adsorbing particles

The adsorption of particles on the confined mixtures modifies their structure
and, largely affect the mixing-demixing transition. Figures 6, 7 and 8 show
selected equilibrium snapshot with $N_{\scriptscriptstyle AP}=1000$ APs for
three representative cases: fully mixed fluids ($\alpha_{\scriptscriptstyle
AB}=1$), fluctuating domains ($\alpha_{\scriptscriptstyle AB}=0.6$) and fully
demixed fluids ($\alpha_{\scriptscriptstyle AB}=0.05$). and we show also, for
comparison, the configurations in the absence of APs. For fully mixed fluids,
the addition of bulk particles with decreasing $\alpha_{\scriptscriptstyle
BP}$ coefficients gives rise to a progressive induced-demixing on the sphere.
This is due to the large energy loss attained when bulk particles adsorb only
on a condensed domain rich in one of the two species on the sphere that hence
tend to segregate. In this case adsorbed particle can be viewed as energetic
bridges between particles of type A.

Figure 6: Snapshots of binary 50:50 AB-mixtures ($N_{p}=500$) confined to move
on a spherical surface with radius of curvature
$R_{s}=7.90\sigma_{\scriptscriptstyle AA}$, with $\alpha_{\scriptscriptstyle
AB}=1.0$ and 3 selected AP-mixture interaction parameter
$\alpha_{\scriptscriptstyle BP}$. A snapshot taken in absence of AP is given
for reference. All snapshots refer to equilibrium states. By decreasing the
affinity between APs and B-particles (snapshots from left to right), particle
adsorption favours particle demixing and domain formation on the sphere, until
reaching a Janus-like configuration for $\alpha_{\scriptscriptstyle
BP}=0.025$.

In the intermediate regime (Figure 7), where large spatial compositional
fluctuations characterize the confined mixture, particle adsorption induces
contrasting effects, depending on whether the interaction of APs with the
particles belonging to the mixtures is symmetric ($\alpha_{\scriptscriptstyle
BP}=1$) or not. When $\alpha_{\scriptscriptstyle BP}=1$ large domains break
up, pushing the mixture towards mixing. Such a non-trivial behavior is the
result of the subtle balance between the configurational entropy of the APs,
that adsorb uniformly on the sphere when $\alpha_{\scriptscriptstyle BP}=1$,
lowering the line tension between the domains and causing their dissolution,
and the energy gain that takes place when more energetic AA (or BB) bond are
replaced by AB bonds.

Figure 7: Snapshots of binary 50:50 AB-mixtures ($N_{p}=500$) confined to move
on a spherical surface with radius of curvature
$R_{s}=7.90\sigma_{\scriptscriptstyle AA}$, with $\alpha_{\scriptscriptstyle
AB}=0.6$ and 3 selected AP-mixture interaction parameter
$\alpha_{\scriptscriptstyle BP}$. A snapshot taken in absence of AP is given
for reference. All snapshots refer to equilibrium states. Particle adsorption
favours particle mixing first ($\alpha_{\scriptscriptstyle BP}=1$) and then,
by decreasing the affinity between APs and B-particles (snapshots from left to
right), domain formation is progressively promoted on the sphere until
reaching a Janus-like configuration ($\alpha_{\scriptscriptstyle BP}=0.025$)

By contrast, when lower values of $\alpha_{\scriptscriptstyle BP}$ are used
and AP interact preferentially with particle of type A, the mixture adopts a
Janus-like (fully separated) configuration, as already discussed for
$\alpha_{\scriptscriptstyle AB}=1$.
Finally for low miscibility (Figure 8) particle adsorption does not affect
much the mixture structure and only a very weak effect can be observed when
$\alpha_{\scriptscriptstyle BP}$ is progressively lowered.

Figure 8: Snapshots of binary 50:50 AB-mixtures ($N_{p}=500$) confined to move
on a spherical surface with radius of curvature
$R_{s}=7.90\sigma_{\scriptscriptstyle AA}$, with $\alpha_{\scriptscriptstyle
AB}=0.05$ and 3 selected AP-mixture interaction parameter
$\alpha_{\scriptscriptstyle BP}$. A snapshot taken in absence of AP is given
for reference. All snapshots refer to equilibrium states. In this case
particle adsorption does not relevantly impact the initially demixed state of
the confined fluid.

We quantify all these changes on the structure of the confined mixture by
computing the total pair correlation function $g(s)$ and the partial
correlation function $g_{m}(s)$ for the AB-mixtures, and the normalized
structural mixing parameter

$\Xi^{*}=\frac{\Xi_{\scriptscriptstyle AP}}{\Xi_{bare}},$ (10)

where $\Xi_{\scriptscriptstyle AP}$ and $\Xi_{bare}$ are calculated according
to equation 7 in the presence and absence of adsorbing particles,
respectively. The visual inspection of the mixture structure can be
substantiated quantitatively in a first instance by comparing $g_{m}(s)$ with
and without APs.

Figure 9: Partial pair correlation functions $g_{m}(s)$ for symmetric AP-
mixture interaction parameter $\alpha_{\scriptscriptstyle BP}$=1 (panels
A,B,C), large preferential affinity of the APs and A-type particles,
$\alpha_{\scriptscriptstyle BP}$=0.025 (panels D,E,F), and three selected
miscibility parameters $\alpha_{\scriptscriptstyle AB}$ as indicated in the
figure. In all the cases $g_{m}(s)$ of the mixtures in absence of APs is also
reported for reference. Insets show the same data in the range $0\leq s\leq 1$
to zoom on the AP effect on the local structure.

This is shown in figure 9 for $\alpha_{\scriptscriptstyle AB}=1$ (panels A,D),
$\alpha_{\scriptscriptstyle AB}=0.6$ (panels B,E) and
$\alpha_{\scriptscriptstyle AB}=0.3$ (panels C,F). We show the two
representative cases of fully symmetric interactions between the APs and two
components of the mixture ($\alpha_{\scriptscriptstyle BP}=1$, panels A,B,C),
and the most asymmetric case investigated ($\alpha_{\scriptscriptstyle
BP}=0.025$, panels D,E,F). The remarkable effect of APs on the mixtures at
$\alpha_{\scriptscriptstyle AB}=0.6$ catches immediately the eye, with
contrasting relative variations passing from symmetric interactions (Figure
9-B) to strongly asymmetric ones (Figure 9-E). This is not observed for the
fully miscible case ($\alpha_{\scriptscriptstyle AB}=1$) or for Janus-like
mixtures ($\alpha_{\scriptscriptstyle AB}=0.3$), for which a decrease in
($\alpha_{\scriptscriptstyle BP}$) gives rise to a progressive induced-
demixing, the extent of the latter being much more pronounced for
$\alpha_{\scriptscriptstyle AB}=1$. In more detail, when large domain
fluctuation dominates and AP addition induces particle mixing partial
correlation functions show augmented local structure peaks and lower values
(closer to one) at large distance with respect to the bare case. This is the
case $\alpha_{\scriptscriptstyle AB}=0.6$ and $\alpha_{\scriptscriptstyle
BP}=1$ (Figure 9-B). On the other hand when AP addition gives rise to induced
demixing, $g_{m}(s)$ show reduced correlation peaks at small geodesic distance
while it increases at large distance with respect to the bare system. This is
the case $\alpha_{\scriptscriptstyle AB}=0.6$ and $\alpha_{\scriptscriptstyle
BP}=0.025$ (Figure 9-E). For the other values of the miscibility parameter,
namely $\alpha_{\scriptscriptstyle AB}=1$ and $\alpha_{\scriptscriptstyle
AB}=0.3$, the AP addition has a relatively weak effect on the structure of the
mixtures when $\alpha_{\scriptscriptstyle BP}=1$ (Figure 9-A,C), while it
produces an important fluid demixing of the initially miscible system
($\alpha_{\scriptscriptstyle AB}=1$) for strongly asymmetric AP-mixture
interactions ($\alpha_{\scriptscriptstyle BP}=0.025$, Figure 9-D). This,
again, is reflected by the net decrease of the structure peaks at small $s$
and a continuous increase of $g_{m}(s)$ at large $s$.

Figure 10: Normalized geodesic mixing parameter $\Xi^{*}$ in function of the
interaction parameter $\alpha_{\scriptscriptstyle BP}$ for the
$\alpha_{\scriptscriptstyle AB}$ values indicated in the figure. Solid lines
are only a guide to the eye.

To build up and visualize more clearly a general scenario based on all our
results we computed $\Xi^{*}(\alpha_{\scriptscriptstyle BP})$ for all the
investigated systems in the presence of AP. This is shown in figure 10 for
$\alpha_{\scriptscriptstyle AB}=1,0.8,0.5,0.3,0.05$. We note that: i) for
large asymmetries of the AP-mixture interaction ($\alpha_{\scriptscriptstyle
BP}\ll 1$) the extent of the induced demixing increases for increasing values
of $\alpha_{\scriptscriptstyle AB}$; ii) the dependence of $\Xi^{*}$ on
$\alpha_{\scriptscriptstyle BP}$ weakens when AB miscibility decreases,
reflecting the fact that the structure is affected by AP addition if particle
segregation is poor; iii) mixtures that are characterized by large fluctuating
domains and interfaces between A-rich domains and B-rich domains show both
induced demixing ($\Xi^{*}>1$) and induced mixing ($\Xi^{*}<1$) depending on
the magnitude of the asymmetry $\alpha_{\scriptscriptstyle BP}$ of the AP-
mixture interactions. This is quite remarkable since many experimental systems
are close to a demixing line and shows that both mixing and demixing are
possible in such systems when particles or complex molecules, other than those
composing the mixtures, are co-suspended.

Figure 11: Void size distribution $P(A_{h}/\sigma_{\scriptscriptstyle
AA}^{2})$ with no preferential affinity ($\alpha_{\scriptscriptstyle BP}=1$)
and varying miscibility parameter $\alpha_{\scriptscriptstyle AB}$ as
indicated in the panels (A,B,C,D,E). Distributions obtained in absence of APs
are also reported for reference. Solid lines are best fits obtained via
equation 9. For each $\alpha_{\scriptscriptstyle AB}$ value one snapshot of an
equilibrium configuration with $\alpha_{\scriptscriptstyle BP}=1$ shows the
biggest visible void. APs are not shown for clarity.

Finally, as done for bare mixtures, we report the void size distributions for
representative systems in the presence of APs and we separate the case
$\alpha_{\scriptscriptstyle BP}=1$ from the rest of the systems since only in
the former case we observe the formation of very large deviations from a
simple exponential decay of the distribution tails. Figure 11 (panels
A,B,C,D,E) shows the distributions $P(A_{h}/\sigma_{\scriptscriptstyle
AA}^{2})$ and five representative configuration snapshots obtained for
$\alpha_{\scriptscriptstyle BP}=1$ and $\alpha_{\scriptscriptstyle
AB}=1,0.6,0.5,0.3,0.05$. The distributions are further compared to those
obtained for bare mixtures. In all the cases, the effect of the AP addition is
to enhance the formation of large voids since particle adsorption increases
the cohesive energy between the particles confined on the sphere. This results
in the suppression of the exponential tails of the void size distributions and
the onset of fat tails, that we attribute to the transition from a regime
where voids are uniformly distributed on the sphere (no APs) to one where one
very large void, whose size fluctuates, dominates. In all such cases we find
tails compatible with a power-law decay ($A_{c}=\infty$) with $P_{0}>1$ for
$\alpha_{\scriptscriptstyle AB}=1$ and $\alpha_{\scriptscriptstyle AB}=0.05$
(Figure 11-A,E). However, such enhancement of large void formation, though
present for each $\alpha_{\scriptscriptstyle AB}$, does depend on the mixture
structure. First of all, a reduced distribution widening is observed for
intermediate $\alpha_{\scriptscriptstyle AB}$, where large compositional
fluctuations characterize the mixtures (see the $\alpha_{\scriptscriptstyle
AB}=0.6$ case in figure 11-B). This is due to the large free energy cost that
a large void would produce, since it would suppress domain fluctuations and
enhance segregation while, in such a regime, the bare mixture is on average
less packed (see Figure 2) and irregular interfaces between domains form.
Thus, the appearance of many voids of intermediate size is favoured at the
expense of one very large void with fluctuating size. Upon decreasing further
the value of $\alpha_{\scriptscriptstyle AB}$ large voids develop again and,
for Janus-like configurations, are spatially located at the interface between
the A-rich domain and the B-rich domain, since the weaker bonds are those
between dissimilar particles. Such interfacial void formation is reflected
again by the extended (non-exponential) tails of the distributions
$P(A_{h}/\sigma_{\scriptscriptstyle AA}^{2})$.

Figure 12: Panels A,B,C,D,E: Void size distribution
$P(A_{h}/\sigma_{\scriptscriptstyle AA}^{2})$ for varying miscibility
parameter $\alpha_{\scriptscriptstyle AB}$ and AP-preferential affinity
($0.025\leq\alpha_{\scriptscriptstyle AB}\leq 0.5$) as indicated in the
figure. Distributions obtained in absence of APs are also reported for
reference. Solid lines are best fits obtained via equation 9 for
$\alpha_{\scriptscriptstyle BP}=0.025$ and for bare mixtures. Bottom panels:
relative variation $\Delta c/c_{bare}$ and $\Delta A_{c}/A_{c}^{bare}$ of the
exponent and the cut-off area defined in equation (9)

For lower $\alpha_{\scriptscriptstyle BP}$ (Figure 12-A,B,C,D,E) the AP
addition affects much less the distributions
$P(A_{h}/\sigma_{\scriptscriptstyle AA}^{2})$, that preserve their exponential
tails. Notwithstanding this, it’s worth noting that AP addition in mixtures at
$\alpha_{\scriptscriptstyle AB}=1$ and $\alpha_{\scriptscriptstyle AB}=0.6$
give rise to contrasting effects: for the fully miscible system, an increase
of the asymmetry of the AP-mixture interaction (i.e. decreasing
$\alpha_{\scriptscriptstyle BP}$) reduces the void sizes for large enough
interaction asymmetry (small $\alpha_{\scriptscriptstyle BP}$), pushing the
distributions towards narrower tails; on the contrary when large irregular
domains characterize the bare mixture, AP addition gives rise always to a
widening of the distribution tails, reflecting the creation of larger voids.
In the first case ($\alpha_{\scriptscriptstyle AB}=1$, Figure 12-A) the void
size reduction is principally due to the fact that, by progressively breaking
the interaction symmetry between particles, APs induce demixing and trigger
the formation of large fluctuating domains. This suppresses, albeit weakly,
large voids akin to the case of the bare mixtures at intermediate
$\alpha_{\scriptscriptstyle AB}$. In the second case
($\alpha_{\scriptscriptstyle AB}=0.6$, Figure 12-B), the enhanced particle
adhesion due to the AP presence dominates, favouring large void formation as
for the aforementioned $\alpha_{\scriptscriptstyle BP}=1$ case. For the other
cases, namely $\alpha_{\scriptscriptstyle AB}<0.6$ (Figure 12-C,D,E), the
effect of APs on $P(A_{h}/\sigma_{\scriptscriptstyle AA}^{2})$, stay very weak
but still detectable. The relative variation of the power-law decay exponent
$\Delta c/c_{bare}$ and the cut-off area $\Delta A_{c}/A_{c}^{bare}$ are
reported in figure 12 (bottom panels), where $\Delta c=(c-c_{bare})$, $\Delta
A_{c}=(A_{c}-A_{c}^{bare})$, and $c_{bare}$ and $A_{c}^{bare}$ are the power-
law exponent and the cut-off area extracted from the distributions in the
absence of APs, respectively. We note that i) the power-law exponent varies
almost to the same extent for any values of $\alpha_{\scriptscriptstyle AB}$
by decreasing the AP-mixture interaction asymmetry (increasing
$\alpha_{\scriptscriptstyle BP}$), pointing in this case to a net reduction of
the small void fraction; ii) the cut-off area variation is maximum for
$\alpha_{\scriptscriptstyle AB}=0.6$ and decreases down to negative values for
$\alpha_{\scriptscriptstyle AB}=1$ and large interaction asymmetry (low
$\alpha_{\scriptscriptstyle BP}$), confirming that the largest relative
variation of the cut-off area, when APs interact preferentially with one
species of particles, occurs for intermediate values of
$\alpha_{\scriptscriptstyle AB}$. This, we stress once more, is due to the
fact that bare mixtures with intermediate $\alpha_{\scriptscriptstyle AB}$
values are characterized the least number of large voids (small cut-off areas)
and the addition of APs with preferential interaction with one of the two
particle types (A here) gives rise to spatially localized (large) voids at the
interface between the formed domains, since AP addition induces particle
demixing as well as an increase of the cohesive energy in the system. The
entire set of best-fit parameters for the data shown in figures 11 and 12 are
in Table 2 (appendix C).

Figure 13: Ratio $g(s_{max})/g(s_{max})|_{bare}$ in function of
$\alpha_{\scriptscriptstyle BP}$ for different miscibility parameters
$\alpha_{\scriptscriptstyle AB}$ as indicated in the figure. Solid lines are a
guide to the eye.

In line with the large variation of the void distribution observed for
miscible systems the local structure of mixtures is most affected by the AP
addition for $\alpha_{\scriptscriptstyle AB}=1$. Figure 13 shows the height of
the first peak of the total pair correlation function $g(s_{max})$ normalized
by its value in the absence of APs $g(s_{max})|_{bare}$. For
$\alpha_{\scriptscriptstyle AB}=1$ this quantity changes the most for varying
$\alpha_{\scriptscriptstyle BP}$ going from values larger than 1 (AP-induced
condensation) for $\alpha_{\scriptscriptstyle BP}=1,0.8$ to values lower than
1 (AP-induced fluidization) for $\alpha_{\scriptscriptstyle BP}\leq 0.5$. In
the other cases $0.05\leq\alpha_{\scriptscriptstyle AB}\leq 0.6$ the local
structure is less affected by the AP adsorption, though we may note that in
this case $g(s_{max})/g(s_{max})|_{bare}$ is always larger than one for
$\alpha_{\scriptscriptstyle BP}=1$, for which the minimum value is attained at
$\alpha_{\scriptscriptstyle AB}=0.6$, where the void size distribution
resulted less affected by AP addition (Figure 11). Finally, we may note that
the normalized peak height shows, albeit very weakly, a systematic non-
monotonic behavior, suggesting a weak melting of the structure for
intermediate values of $\alpha_{\scriptscriptstyle BP}$. More detailed
simulations would be required to reach more quantitative conclusions on this
aspect.

##  VI Conclusions

In this paper we have studied the mixing-demixing transition of symmetric
mixtures characterized by varying miscibility and confined to move on a sphere
by using a minimal quasi-2D model. By computing the total and the partial pair
correlation functions, and by employing a geodesic mixing parameter $\Xi$, we
have shown that: i) demixing on 3D-spheres is a structurally smooth process,
in which fully mixed and Janus-like configurations are separated by states of
spatially fluctuating domains appearing for intermediate miscibility; ii) void
size distributions show exponential tails preceded by a power-law decay with
large void formation being hampered by large compositional fluctuations, that
instead favor the appearance of small ($\sim$one-particle sized) voids. Such a
scenario does not vary remarkably in the range of system sizes investigated
($R_{s}=11.18\sigma_{\scriptscriptstyle AA},7.90\sigma_{\scriptscriptstyle
AA},5.59\sigma_{\scriptscriptstyle AA}$). iii) Particle adsorption can give
rise to both induced-mixing and induced-demixing, depending on whether
mixtures are characterized by strongly fluctuating inter-domain boundaries or
not. On the one hand, when coexisting domains and fluctuating interfaces
characterize the fluids, uniformly adsorbed particles
($\alpha_{\scriptscriptstyle BP}=1$) enhance miscibility, while the latter is
reduced for a large preferential affinity of APs for one of the two components
of the mixtures ($\alpha_{\scriptscriptstyle BP}\ll 1$). On the other hand,
for fully miscible ($\alpha_{\scriptscriptstyle AB}=1$) or Janus-like mixtures
($\alpha_{\scriptscriptstyle AB}\ll 1$) AP addition mainly produces particle
demixing, whose extent, quantified by the normalized geodesic parameter
$\Xi^{*}$, increases progressively with increasing miscibility, reaching its
apex for $\alpha_{\scriptscriptstyle AB}=1$ and $\alpha_{\scriptscriptstyle
BP}\ll 1$. Void size distributions are also affected by AP addition. A
transition from exponentially-tailed to fat-tailed distributions, where one
single large void dominates over the others, is observed for uniform particle
adsorption, namely where AP-mixture cohesion energy is maximum. The extent of
such a transition also depends on the AB miscibility: fat-tailed distributions
are less pronounced by large domain fluctuations occurring for intermediate
$\alpha_{\scriptscriptstyle AB}$ values, while for large enough AP-mixture
interaction asymmetry, we recover exponentially-tailed distributions, with the
largest tail widening and narrowing obtained respectively for an intermediate
miscibility ($\alpha_{\scriptscriptstyle AB}$=0.6) and fully miscible fluids
($\alpha_{\scriptscriptstyle AB}$=1). Our model, albeit simple, allowed to
discern under which mixing conditions a quasi 2D multicomponent fluid confined
to move on a sphere is more susceptible to incur structural changes due to the
interaction with an external agent. Further studies, based on the same model,
will elucidate the role of mixture composition, temperature and AP
concentration, while the effect of local curvature can be inspected by
changing the confining manifold, passing from spheres to ellipsoids or tori.
We hope that our results can pave the way for more targeted experiments to
investigate the role of complex molecules and particles on the mixing state of
quasi-2D liquids.

###### Acknowledgements.

I wish to thank Dr. Simona Sennato for inspiring discussions. We acknowledge
financial support from the Agence Nationale de la Recherche (Grant
ANR-20-CE06-0030-01; THELECTRA).

## Data Availability Statement

The data that support the findings of this study are all available within the
article.

##  Appendix A Mean squared angular displacement on the sphere

We have performed long NVT simulations ($5\cdot 10^{8}$ and $10^{9}$ time
steps) for the two largest system investigated ($N_{p}=500,1000$) and
$\alpha_{\scriptscriptstyle AB}=1$ to compute the angular mean square
displacements (MSD) $\langle(\theta(t)-\theta_{0})^{2}\rangle_{p}^{2}$ and
$\langle(\varphi(t)-\varphi_{0})^{2}\rangle_{p}$, and to check whether our
model produces diffusive trajectories over large runs. Here $\theta(t)$ and
$\varphi(t)$ represent the polar and the azimuthal angle of each particle
confined on the sphere at time $t>0$, $\theta_{0}$ and $\varphi_{0}$ are their
values at the initial time $t=0$, and the brackets $\langle\cdot\rangle_{p}$
stand for the average over all particles confined on the sphere. Figure 14
shows in a log-log scale both the polar and the azimuthal squared
displacements, for which we observe 2 different saturation values as expected
for confined dynamics. Such values can be computed analytically. The angular
MSD of one Brownian particle on the sphere saturate to values that are
dependent on their initial position ($\theta_{0},\varphi_{0}$) Apaza and
Sandoval (2017):

$\displaystyle\overline{(\theta(t)-\theta_{0})^{2}}=\frac{\pi^{2}-4}{2}+\theta_{0}(\theta_{0}-\pi)$
(11)
$\displaystyle\overline{(\varphi(t)-\varphi_{0})^{2}}=\frac{4\pi^{2}}{3}+\varphi_{0}(\varphi_{0}-2\pi)$
(12)

where $\overline{(\cdot)}$ is the average over different trajectories given an
initial angular position ($\theta_{0},\varphi_{0}$). For a crowded system
particles are distributed over the whole spherical surface and we must compute
the saturation values of the MSD as the angular average of equations 11 and
12, i.e.:

$\displaystyle\langle(\theta-\theta_{0})^{2}\rangle_{p}=\int_{0}^{2\pi}d\varphi_{0}\int_{0}^{\pi}\overline{(\theta-\theta_{0})^{2}}\sin(\theta_{0})d\theta_{0}=\frac{\pi-8}{2}$
(13)
$\displaystyle\langle(\varphi-\varphi_{0})^{2}\rangle_{p}=\int_{0}^{\pi}\sin(\theta_{0})d\theta_{0}\int_{0}^{2\pi}\overline{(\varphi-\varphi_{0})^{2}}d\varphi_{0}=\frac{2\pi^{2}}{3}$
(14)

The values obtained are in very good agreement with our simulations. We also
note that the polar MSD is a linear function of $t$ for small lag times in
accordance with recent theoretical works Apaza and Sandoval (2017); Castañeda-
Priego  _et al._ (2013), while the azimuthal MSD shows a subdiffusive behavior
$\langle(\varphi-\varphi_{0})^{2}\rangle_{p}\sim t^{0.47\pm 0.04}$, that, up
to our knoweldge, has not been observed so far. $\varphi(t)$ in a spherical
set of coordinates isn’t in fact a variable satisfying the Langevin equation
and we attribute this feature to the space-dependent azimuthal diffusivity
Apaza and Sandoval (2017): particles close to the pole $\theta=0$ move faster
in $\varphi$ with respect to particles moving around the equator
$\theta=\pi/2$ for a fixed geodesic displacement. Sub-diffusive motions are
often observed in experiments and simulations in systems characterized by
space-dependent diffusion constants Cherstvy, Chechkin, and Metzler (2014).
This aspect, interesting _per se_ , goes beyond the scope of this paper and
surely deserves a deeper investigation. Finally the displacements obtained for
the two systems confined onto spheres of different radius $R_{s}$ superimpose
when a rescaled time $t/R^{2}$ is considered, in agreement with the general
solution of the diffusion equation for spherical Brownian motion. Yosida
(1949); Castañeda-Priego  _et al._ (2013).

Figure 14: Azimuthal and polar mean squared displacement in function of the
rescaled time $t/R_{s}^{2}$ for $N_{p}=500$ and $N_{p}=1000$. Dashed
horizontal lines indicate the saturation values obtained analytically (Eqs. 13
and 14).

##  Appendix B Fit of the geodesic mixing parameter
$\Xi(\alpha_{\scriptscriptstyle AB})$

We extracted the critical mixing parameter $\alpha_{\scriptscriptstyle
AB}^{c}$ from fits of the geodesic mixing parameter
$\Xi(\alpha_{\scriptscriptstyle AB})$ through the phenomenological function

$\Xi(\alpha_{\scriptscriptstyle AB})=\alpha_{\scriptscriptstyle
AB}^{0}+\Delta\alpha_{\scriptscriptstyle
AB}\tanh\left[(\psi(\alpha_{\scriptscriptstyle AB}-\alpha_{\scriptscriptstyle
AB}^{c})\right]\\\ +A(\alpha_{\scriptscriptstyle
AB}-\alpha_{\scriptscriptstyle AB}^{c})+B(\alpha_{\scriptscriptstyle
AB}-\alpha_{\scriptscriptstyle AB}^{c})^{2}$ (15)

We obtain the critical miscibility parameters $\alpha_{\scriptscriptstyle
AB}^{c}$ reported in Table 1.

$N_{p}$ | $\alpha_{\scriptscriptstyle AB}^{c}$
---|---
250 | 0.40 $\pm$ 0.02
500 | 0.43 $\pm$ 0.04
1000 | 0.49 $\pm$ 0.01
Table 1: Critical miscibility parameter $\alpha_{\scriptscriptstyle AB}^{c}$
obtained by fitting $\Xi(\alpha_{\scriptscriptstyle AB})$ via equation 15.

##  Appendix C Best-fit parameters for $P(A_{h}/\sigma_{\scriptscriptstyle
AA}^{2})$

We report below (Table 2) the best-fit parameters obtained by fitting the void
size distributions shown in figure 11 and 12 via equation (9).

$\alpha_{\scriptscriptstyle AB}$ | $\alpha_{\scriptscriptstyle BP}$ | $P_{0}$ | $a$ | $b$ | $c$ | $A_{c}/\sigma_{\scriptscriptstyle AA}^{2}$
---|---|---|---|---|---|---
1 | 1 | $2.95\cdot 10^{-5}$ | 1.95 | 4.68 | 2.24 | $\infty$
0.6 | 1 | 0 | 1.25 | 3.28 | 2.33 | $\infty$
0.5 | 1 | 0 | 0.28 | 6.04 | 2.48 | $\infty$
0.3 | 1 | 0 | 1.53 | 4.08 | 2.12 | $\infty$
0.05 | 1 | $4.5\cdot 10^{-5}$ | 1.40 | 3.64 | 2.32 | $\infty$
1 | 0.5 | 0 | 1.35 | 3.52 | 2.10 | 9.17
0.6 | 0.5 | 0 | 1.40 | 3.56 | 2.07 | 9.24
0.5 | 0.5 | 0 | 1.34 | 3.72 | 2.01 | 8.68
0.3 | 0.5 | 0 | 1.35 | 3.56 | 1.97 | 6.32
0.05 | 0.5 | 0 | 1.44 | 3.88 | 1.95 | 8.92
1 | 0.25 | 0 | 1.38 | 3.6 | 2.06 | 7.66
0.6 | 0.25 | 0 | 1.33 | 3.48 | 1.99 | 5.73
0.5 | 0.25 | 0 | 1.57 | 4.24 | 1.88 | 7.45
0.3 | 0.25 | 0 | 1.52 | 4.13 | 1.90 | 7.75
0.05 | 0.25 | 0 | 1.65 | 4.52 | 1.79 | 6.00
1 | 0.1 | 0 | 1.41 | 3.65 | 2.02 | 5.50
0.6 | 0.1 | 0 | 1.76 | 4.56 | 1.89 | 6.00
0.5 | 0.1 | 0 | 1.70 | 4.56 | 1.78 | 5.00
0.3 | 0.1 | 0 | 1.74 | 4.72 | 1.74 | 5.00
0.05 | 0.1 | 0 | 1.84 | 5.08 | 1.74 | 6.00
1 | 0.025 | 0 | 1.62 | 4.22 | 1.91 | 5.50
0.6 | 0.025 | 0 | 1.55 | 4.14 | 1.90 | 6.25
0.5 | 0.025 | 0 | 1.58 | 4.24 | 1.88 | 6.50
0.3 | 0.025 | 0 | 1.75 | 4.80 | 1.76 | 5.75
0.05 | 0.025 | 0 | 1.49 | 4.16 | 1.80 | 6.00
Table 2: Best fit parameters obtained via equation 9 for the void size
distributions shown in figures 11 and 12.

##  Appendix D Sampling time

The sampling time ($\Delta t_{s}=5\cdot 10^{5}$ time steps) has been chosen by
computing the (spherical) self-intermediate scattering function for the bare
mixture with the lowest energy ($N_{p}=1000$ and $\alpha_{\scriptscriptstyle
AB}=1$) as defined in Singh, Sood, and Ganapathy (2020):

${\color[rgb]{0,0,0}F_{s}(k,t)=\frac{1}{N_{p}}\sum\limits_{j=1}^{N_{p}}\left\langle
P_{kR}\left[\cos\left(\frac{\Delta s_{j}(0,t)}{R}\right)\right]\right\rangle}$
(16)

Where, $P_{kR}$ is the Legendre polynomial with kR being rounded-off to the
nearest integer, $k^{-1}$ is the wavevector corresponding to the geodesic
distance equal to one particle diameter and $\Delta s_{j}(0,t)$ is the
geodesic displacement of particle $j$ over time $t$. We obtain the scattering
function shown in figure 15 that is well fitted by a single exponential decay
with a relaxation time of $8.29\cdot 10^{4}$ time steps that is about 6 time
less than the fixed $\Delta t_{s}$.

Figure 15: Self-intermediate scattering function at a wavevector corresponding
to a geodesic distance equal to one particle diameter for $N_{p}=1000$ and
$\alpha_{\scriptscriptstyle AB}=1$.

## References

  * Malchus and Weiss (2010) N. Malchus and M. Weiss, “Anomalous Diffusion Reports on the Interaction of Misfolded Proteins with the Quality Control Machinery in the Endoplasmic Reticulum,” Biophysical Journal 99, 1321–1328 (2010).
  * Javanainen  _et al._ (2013) M. Javanainen, H. Hammaren, L. Monticelli, J.-H. Jeon, M. S. Miettinen, H. Martinez-Seara, R. Metzler, and I. Vattulainen, “Anomalous and normal diffusion of proteins and lipids in crowded lipid membranes,” Faraday Discuss. 161, 397–417 (2013).
  * Metzler, Jeon, and Cherstvy (2016) R. Metzler, J.-H. Jeon, and A. Cherstvy, “Non-Brownian diffusion in lipid membranes: Experiments and simulations,” Biochimica et Biophysica Acta (BBA) - Biomembranes 1858, 2451–2467 (2016).
  * Savarala  _et al._ (2011) 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).
  * Anderson (2002) R. G. W. Anderson, “A Role for Lipid Shells in Targeting Proteins to Caveolae, Rafts, and Other Lipid Domains,” Science 296, 1821–1825 (2002).
  * Binder, Barragan, and Menger (2003) W. H. Binder, V. Barragan, and F. M. Menger, “Domains and Rafts in Lipid Membranes,” Angew. Chem. Int. Ed. 42, 5802–5827 (2003).
  * Veatch  _et al._ (2008) S. L. Veatch, P. Cicuta, P. Sengupta, A. Honerkamp-Smith, D. Holowka, and B. Baird, “Critical Fluctuations in Plasma Membrane Vesicles,” ACS Chem. Biol. 3, 287–293 (2008).
  * Bowick and Giomi (2009) M. J. Bowick and L. Giomi, “Two-dimensional matter: Order, curvature and defects,” Advances in Physics 58, 449–563 (2009).
  * Singh, Sood, and Ganapathy (2020) N. Singh, A. K. Sood, and R. Ganapathy, “Cooperatively rearranging regions change shape near the mode-coupling crossover for colloidal liquids on a sphere,” Nat Commun 11, 4967 (2020).
  * Post and Glandt (1986) A. J. Post and E. D. Glandt, “Statistical thermodynamics of particles adsorbed onto a spherical surface. I. Canonical ensemble,” The Journal of Chemical Physics 85, 7349–7358 (1986).
  * Salaün, James, and Chamberlain (2004) C. Salaün, D. J. James, and L. H. Chamberlain, “Lipid Rafts and the Regulation of Exocytosis: Lipid Rafts and Exocytosis,” Traffic 5, 255–264 (2004).
  * Tang  _et al._ (2005) P. Tang, F. Qiu, H. Zhang, and Y. Yang, “Phase separation patterns for diblock copolymers on spherical surfaces: A finite volume method,” Phys. Rev. E 72, 016710 (2005).
  * Greer, Bertozzi, and Sapiro (2006) J. B. Greer, A. L. Bertozzi, and G. Sapiro, “Fourth order partial differential equations on general geometries,” Journal of Computational Physics 216, 216–246 (2006).
  * Marenduzzo and Orlandini (2013) D. Marenduzzo and E. Orlandini, “Phase separation dynamics on curved surfaces,” Soft Matter 9, 1178–1187 (2013).
  * Barreira, Elliott, and Madzvamuse (2011) R. Barreira, C. M. Elliott, and A. Madzvamuse, “The surface finite element method for pattern formation on evolving biological surfaces,” J. Math. Biol. 63, 1095–1119 (2011).
  * Saffman and Delbruck (1975) P. G. Saffman and M. Delbruck, “Brownian motion in biological membranes.” Proceedings of the National Academy of Sciences 72, 3111–3113 (1975).
  * Camley and Brown (2011) B. A. Camley and F. L. H. Brown, “Dynamic scaling in phase separation kinetics for quasi-two-dimensional membranes,” The Journal of Chemical Physics 135, 225106 (2011).
  * Lenne and Nicolas (2009) P.-F. Lenne and A. Nicolas, “Physics puzzles on membrane domains posed by cell biology,” Soft Matter 5, 2841 (2009).
  * Lingwood and Simons (2010) D. Lingwood and K. Simons, “Lipid Rafts As a Membrane-Organizing Principle,” Science 327, 46–50 (2010).
  * Onuki (2002) A. Onuki,  _Phase Transition Dynamics_, 1st ed. (Cambridge University Press, 2002).
  * Honerkamp-Smith, Veatch, and Keller (2009) A. R. Honerkamp-Smith, S. L. Veatch, and S. L. Keller, “An introduction to critical points for biophysicists; observations of compositional heterogeneity in lipid membranes,” Biochimica et Biophysica Acta (BBA) - Biomembranes 1788, 53–63 (2009).
  * Destainville, Manghi, and Cornet (2018) N. Destainville, M. Manghi, and J. Cornet, “A Rationale for Mesoscopic Domain Formation in Biomembranes,” Biomolecules 8, 104 (2018).
  * Keller  _et al._ (1998) S. L. Keller, W. H. Pitcher, W. H. Huestis, and H. M. McConnell, “Red Blood Cell Lipids Form Immiscible Liquids,” Phys. Rev. Lett. 81, 5019–5022 (1998).
  * Bianchi, Blaak, and Likos (2011) E. Bianchi, R. Blaak, and C. N. Likos, “Patchy colloids: State of the art and perspectives,” Phys. Chem. Chem. Phys. 13, 6397 (2011).
  * Ai, Zhou, and Zhang (2020) B.-q. Ai, B.-y. Zhou, and X.-m. Zhang, “Binary mixtures of active and passive particles on a sphere,” Soft Matter 16, 4710–4717 (2020).
  * Shankar, Bowick, and Marchetti (2017) S. Shankar, M. J. Bowick, and M. C. Marchetti, “Topological Sound and Flocking on Curved Surfaces,” Phys. Rev. X 7, 031039 (2017).
  * Keber  _et al._ (2014) F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic, and A. R. Bausch, “Topology and dynamics of active nematic vesicles,” Science 345, 1135–1139 (2014).
  * Sanchez  _et al._ (2012) T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, “Spontaneous motion in hierarchically assembled active matter,” Nature 491, 431–434 (2012).
  * Caetano  _et al._ (2020) D. L. Z. Caetano, S. J. de Carvalho, R. Metzler, and A. G. Cherstvy, “Critical adsorption of multiple polyelectrolytes onto a nanosphere: Splitting the adsorption–desorption transition boundary,” J. R. Soc. Interface. 17, 20200199 (2020).
  * Caetano  _et al._ (2017) D. L. Z. Caetano, S. J. de Carvalho, R. Metzler, and A. G. Cherstvy, “Critical adsorption of periodic and random polyampholytes onto charged surfaces,” Phys. Chem. Chem. Phys. 19, 23397–23413 (2017).
  * Jacobson and Papahadjopoulos (1975) K. Jacobson and D. Papahadjopoulos, “Phase transitions and phase separations in phospholipid membranes induced by changes in temperature, pH, and concentration of bivalent cations,” Biochemistry 14, 152–161 (1975).
  * Michanek  _et al._ (2010) A. Michanek, N. Kristen, F. Höök, T. Nylander, and E. Sparr, “RNA and DNA interactions with zwitterionic and charged lipid membranes — A DSC and QCM-D study,” Biochimica et Biophysica Acta (BBA) - Biomembranes 1798, 829–838 (2010).
  * Mbamala, Ben-Shaul, and May (2005) E. C. Mbamala, A. Ben-Shaul, and S. May, “Domain Formation Induced by the Adsorption of Charged Proteins on Mixed Lipid Membranes,” Biophysical Journal 88, 1702–1714 (2005).
  * Macdonald  _et al._ (1998) P. M. Macdonald, K. J. Crowell, C. M. Franzin, P. Mitrakos, and D. J. Semchyschyn, “Polyelectrolyte-induced domains in lipid bilayer membranes: The deuterium NMR perspective,” 76, 13 (1998).
  * Canepa  _et al._ (2020) E. Canepa, S. Salassi, A. L. de Marco, C. Lambruschini, D. Odino, D. Bochicchio, F. Canepa, C. Canale, S. Dante, R. Brescia, F. Stellacci, G. Rossi, and A. Relini, “Amphiphilic gold nanoparticles perturb phase separation in multidomain lipid membranes,” Nanoscale 12, 19746–19759 (2020).
  * Bochicchio  _et al._ (2017) D. Bochicchio, E. Panizon, L. Monticelli, and G. Rossi, “Interaction of hydrophobic polymers with model lipid bilayers,” Sci Rep 7, 6357 (2017).
  * Kotyk and Janáček (1977) A. Kotyk and K. Janáček,  _Membrane Transport_ (Springer US, Boston, MA, 1977).
  * Kashchiev and Exerowa (1983) D. Kashchiev and D. Exerowa, “Bilayer lipid membrane permeation and rupture due to hole formation,” Biochimica et Biophysica Acta (BBA) - Biomembranes 732, 133–145 (1983).
  * LaRocca  _et al._ (2013) T. J. LaRocca, P. Pathak, S. Chiantia, A. Toledo, J. R. Silvius, J. L. Benach, and E. London, “Proving Lipid Rafts Exist: Membrane Domains in the Prokaryote Borrelia burgdorferi Have the Same Properties as Eukaryotic Lipid Rafts,” PLoS Pathog 9, e1003353 (2013).
  * Ghysels  _et al._ (2019) A. Ghysels, A. Krämer, R. M. Venable, W. E. Teague, E. Lyman, K. Gawrisch, and R. W. Pastor, “Permeability of membranes in the liquid ordered and liquid disordered phases,” Nat Commun 10, 5616 (2019).
  * Gensure, Zeidel, and Hill (2006) R. H. Gensure, M. L. Zeidel, and W. G. Hill, “Lipid raft components cholesterol and sphingomyelin increase H+/OH- permeability of phosphatidylcholine membranes,” Biochemical Journal 398, 485–495 (2006).
  * Sciolla  _et al._ (2021) F. Sciolla, D. Truzzolillo, E. Chauveau, S. Trabalzini, L. Di Marzio, M. Carafa, C. Marianecci, A. Sarra, F. Bordi, and S. Sennato, “Influence of drug/lipid interaction on the entrapment efficiency of isoniazid in liposomes for antitubercular therapy: A multi-faced investigation,” Colloids and Surfaces B: Biointerfaces 208, 112054 (2021).
  * Ha (2001) B.-Y. Ha, “Stabilization and destabilization of cell membranes by multivalent ions,” Phys. Rev. E 64, 051902 (2001).
  * (44) B. V. Derjaguin and A. V. Prokhorov, “On the Theory of the Rupture of Black Films,” , 8.
  * Kashchiev and Exerowa (1980) D. Kashchiev and D. Exerowa, “Nucleation mechanism of rupture of newtonian black films. I. Theory,” Journal of Colloid and Interface Science 77, 501–511 (1980).
  * Dickinson (1978) E. Dickinson, “A hard disk fluid model of monolayer permeation and evaporation resistance,” Journal of Colloid and Interface Science 63, 461–471 (1978).
  * Stoeckly (1977) B. Stoeckly, “Equation of state of fatty-acid monolayers on water,” Phys. Rev. A 15, 2558–2562 (1977).
  * Wilding, Schmid, and Nielaba (1998) N. B. Wilding, F. Schmid, and P. Nielaba, “Liquid-vapor phase behavior of a symmetrical binary fluid mixture,” Phys. Rev. E 58, 2201–2212 (1998).
  * Díaz-Herrera, Moreno-Razo, and Ramírez-Santiago (2004) E. Díaz-Herrera, J. A. Moreno-Razo, and G. Ramírez-Santiago, “Wetting phenomenon in the liquid-vapor phase coexistence of a partially miscible Lennard-Jones binary mixture,” Phys. Rev. E 70, 051601 (2004).
  * Komura  _et al._ (2004) S. Komura, H. Shirotori, P. D. Olmsted, and D. Andelman, “Lateral phase separation in mixtures of lipids and cholesterol,” Europhys. Lett. 67, 321–327 (2004).
  * Almeida, Pokorny, and Hinderliter (2005) P. F. Almeida, A. Pokorny, and A. Hinderliter, “Thermodynamics of membrane domains,” Biochimica et Biophysica Acta (BBA) - Biomembranes 1720, 1–13 (2005).
  * Peliti (1991) L. Peliti, ed.,  _Biologically Inspired Physics_, NATO ASI Series, Vol. 263 (Springer US, Boston, MA, 1991).
  * Barker, Henderson, and Abraham (1981) J. Barker, D. Henderson, and F. Abraham, “Phase diagram of the two-dimensional Lennard-Jones system; Evidence for first-order transitions,” Physica A: Statistical Mechanics and its Applications 106, 226–238 (1981).
  * Holian and Grady (1988) B. L. Holian and D. E. Grady, “Fragmentation by molecular dynamics: The microscopic “big bang”,” Phys. Rev. Lett. 60, 1355–1358 (1988).
  * Israelachvili and Mitchell (1975) J. N. Israelachvili and D. Mitchell, “A model for the packing of lipids in bilayer membranes,” Biochimica et Biophysica Acta (BBA) - Biomembranes 389, 13–19 (1975).
  * (56) M. Deserno, “How to generate equidistributed points on the surface of a sphere,” , 1.
  * Boutreux and de Geennes (1997) T. Boutreux and P. de Geennes, “Compaction of granular mixtures: A free volume model,” Physica A: Statistical Mechanics and its Applications 244, 59–67 (1997).
  * Caglioti  _et al._ (1999) E. Caglioti, A. Coniglio, H. J. Herrmann, V. Loreto, and M. Nicodemi, “Cooperative length approach for granular media,” Physica A: Statistical Mechanics and its Applications 265, 311–318 (1999).
  * Amitrano (2012) D. Amitrano, “Variability in the power-law distributions of rupture events: How and why does b-value change,” Eur. Phys. J. Spec. Top. 205, 199–215 (2012).
  * Treado  _et al._ (2019) J. D. Treado, Z. Mei, L. Regan, and C. S. O’Hern, “Void distributions reveal structural link between jammed packings and protein cores,” Phys. Rev. E 99, 022416 (2019).
  * Jorjadze  _et al._ (2011) I. Jorjadze, L.-L. Pontani, K. A. Newhall, and J. Brujic, “Attractive emulsion droplets probe the phase diagram of jammed granular matter,” Proceedings of the National Academy of Sciences 108, 4286–4291 (2011).
  * Apaza and Sandoval (2017) L. Apaza and M. Sandoval, “Brownian self-driven particles on the surface of a sphere,” Phys. Rev. E 96, 022606 (2017).
  * Castañeda-Priego  _et al._ (2013) R. Castañeda-Priego, P. Castro-Villarreal, S. Estrada-Jiménez, and J. M. Méndez-Alcaraz, “Brownian motion of free particles on curved surfaces,” (2013), arXiv:1211.5799 [cond-mat, physics:math-ph] .
  * Cherstvy, Chechkin, and Metzler (2014) A. G. Cherstvy, A. V. Chechkin, and R. Metzler, “Particle invasion, survival, and non-ergodicity in 2D diffusion processes with space-dependent diffusivity,” Soft Matter 10, 1591–1601 (2014).
  * Yosida (1949) K. Yosida, “Brownian Motion on the Surface of the 3-Sphere,” Ann. Math. Statist. 20, 292–296 (1949).