Equilibrium gels of colloidal particles can be realized through the introduction
of a second species, a linker that mediates the bonds between colloids. A gel
forming binary mixture whose linkers can self-assemble into linear chains while
still promoting the aggregation of particles is considered in this work. The
particles are patchy particles with *f*_{C} patches of type *C* and the linkers are patchy particles with 2 patches of
type *A* and *f*_{B} patches of type B. The bonds between patches of type *A* (*AA* bonds) promote the formation of linear chains of
linkers. Two different ways (model A and model B) of bonding the linkers to the
particles—or inducing branching—are studied. In model A, there is a competition
between chaining and branching, since the bonding between linkers and particles
takes place through *AC* bonds only. In model B, the linkers
aggregate to particles through bonds *BC* only, making chaining
and branching independent. The percolation behavior of these two models is
studied in detail, employing a generalized Flory–Stockmayer theory and Monte
Carlo simulations. The self-assembly of linkers into chains reduces the fraction
of particles needed for percolation to occur (models A and B) and induces
percolation when the fraction of particles is high (model B). Percolation by
heating and percolation loops in temperature–composition diagrams are obtained
when the formation of chains is energetically favorable by increasing the
entropic gain of branching (model A). Chaining and branching are found to follow
a model dependent relation at percolation, which shows that, for the same
composition, longer chains require less branching for percolation to occur.

## I. INTRODUCTION

The theoretical and simulation studies of patchy particle models^{1–4} have opened the way to the concept of
equilibrium gels, a thermodynamically stable phase formed by a percolated network of
particles, not related to phase separation arrest. These works have shown that the
formation of bonds between particles with low functionality leads to the emergence
of a percolated fluid at low temperatures and low densities, which does not phase
separate with a vapor. Experimental evidence of the existence of equilibrium gels
was first found in single component systems, like laponite^{5} (a colloidal clay), DNA nanostars,^{6} and a solution of
Fmoc-diphenylalanine molecules in dimethyl sulfoxide.^{7}

In recent years, significant attention has been given to the study of aggregation in
binary mixtures where the self-assembly of one of the components is mediated and
controlled by the other. Examples span a large variety of soft matter and biological
systems, such as, among many others, protein–protein aggregation,^{8} cross-linking of actin
filaments,^{9} colloidal
dispersions and protein fibrils,^{10} nanoparticles and globular proteins,^{11} microparticles and small soft microgels,^{12} and cell-mediated colloidal
scaffolds.^{13–15} Following a similar line of reasoning, the search for realizing equilibrium gels
with an extra degree of control has led to the study of model systems where bonds
between particles are mediated by another component, linkers.^{16–20} In these works, the linkers
are bifunctional (i.e., they bond two particles) and can be polymers of different
lengths,^{17} DNA strands,^{18,19} or patchy particles.^{16,20} In any case, their
properties (size, shape, and concentration) can be used to control the aggregation
process and the thermodynamics, and they can be adjusted to obtain the proper
conditions for the formation of equilibrium gels. These works establish some
interesting results for linker–particle aggregating systems: A percolated network
appears at low temperatures in a finite range of linker concentrations^{16,20} that depends on the
functionality of particles; the connectivity properties of this network can be
controlled by the amount of linkers;^{16,17,20} it is possible to find single phase percolated
fluids at low densities and temperatures,^{16,17,20} and these densities can be further reduced using
longer linkers.^{17} An experimental
system^{21} formed by
nanoparticles and telechelic polymer chains (the linkers) confirms some of these
predictions since it is found that the gelation of particles is controlled by the
relative concentration of polymeric linkers.

The aim of the present work is to thoroughly investigate the percolation thresholds
of linker–particle patchy models in which linkers can also assemble into linear
chains, using theory and simulation. This is accomplished by introducing bonds
between linkers. While linker–linker bonds promote chaining, the usual
particle–linker bonds lead to branching of those chains and eventually to the
formation of percolated networks. Two types of interplay between chaining and
branching are addressed (in two variants of this model): (i) A competition between
chaining and branching is set by letting the two patches of the linkers bond either
to another linker or to a particle (model A). (ii) Chaining and branching are set
independently by letting the linkers have two patches that bond only to linkers and
other patches that bond only to particles (model B). It will be shown that the
chaining of linkers strongly affects the conditions at which percolation occurs and
that these changes are different in models A and B. It should be emphasized that
percolation is only a necessary condition to obtain equilibrium gels. The
investigation of the phase diagrams of the percolated phases of these models is left
for future work. Still, it is important to stress that in almost all models of
binary mixtures of patchy particles with low valence,^{16,20,22–24} the phase diagrams
(calculated using Wertheim’s theory) always exhibit percolated single phases at low
temperatures for a range of low densities.

The paper is organized as follows: In Sec. II, the patchy particle models are introduced in detail, and the simulation and theoretical methods employed to study the percolation properties of the models are described. In Sec. III, the percolation thresholds obtained for the two models are described and analyzed. Finally, in Sec. IV, the results are discussed and conclusions are drawn.

## II. MODEL

A binary mixture of *N*_{P} particles and *N*_{L} linkers in a volume V is
considered. Particles and linkers are hard spheres (HSs) of diameter *σ*. The (reduced) density of the system is *ρ* =
(*N*_{L} + *N*_{P})*σ*^{3}/*V* and the composition *x* is the fraction of particles, *x* = *N*_{P}/(*N*_{L} + *N*_{P}). Particles are decorated with *f*_{C} ≥ 3 patches of type C on its
surface, while linkers are decorated with 2 patches of type A and up to *f*_{B} patches of type B (see Fig. 1).

*V*

_{ij}between HSs

*i*and

*j*is

*V*

_{HS}is the HS potential,

*V*

_{αβ}is a spherical-well potential (of energy depth −

*ɛ*

_{αβ}≤ 0 and range

*δ*

_{αβ}) corresponding to the interaction between a patch of type

*α*on HS

*i*and a patch of type

*β*on HS

*j*from the set Γ

_{ij}of possible such pairs,

*r*

_{ij}is the distance between HSs

*i*and

*j*, and $rij\alpha \beta $ is the distance between the center of patch

*α*in HS

*i*and the center of patch

*β*in HS

*j*.

^{25}Essentially, a bond

*αβ*between HSs

*i*and

*j*is established and the potential energy decreases by

*ɛ*

_{αβ}when $rij\alpha \beta <\delta \alpha \beta $. In simulations, the placement of patches over the HS and the ranges

*δ*

_{αβ}of the patch–patch potential are chosen so that it is guaranteed that each patch engages at most in a single bond.

^{25}

The interaction potential defined in (1) depends on the number of patches of different types on each species
(*f*_{C} and *f*_{B}) and on the energy scales *ɛ*_{αβ} between all pairs of types of
patches. The choice of *f*_{C}, *f*_{B} and of the nonzero *ɛ*_{αβ} (i.e., of which types of bonds
can be formed) will define a particular model of a binary mixture of patchy
particles. The types of bonds allowed will then determine the types of
self-assembled structures that emerge on the model. This relation is illustrated in Fig. 2. If the only nonzero energy is *ɛ*_{AA}, then only bonds between
linkers can be formed, and chains of linkers are the only self-assembled structures
expected—Fig. 2(a). The formation of
clusters where a single linker mediates the bonds between two (and only two)
particles^{16,20} is obtained
for a model where only *ɛ*_{AC} ≠ 0 (or only *ɛ*_{AB} ≠ 0 and *f*_{B} = 2)—Fig. 2(b).

The goal of this work is to study the percolation thresholds in systems where the linkers, besides mediating the bonding of particles, can also self-assemble into chains. This will be accomplished by investigating the following models:

Model A: Only

*AA*and*AC*bonds can form, i.e., only*ɛ*_{AC}and*ɛ*_{AA}are nonzero; in simulations (and numerical calculations),*f*_{C}= 4, and several values of the energy ratio*ɛ*_{AA}/*ɛ*_{AC}are considered.Model B: Only

*AA*and*BC*bonds can form, i.e., only*ɛ*_{BC}and*ɛ*_{AA}are nonzero; in simulations (and numerical calculations),*f*_{C}= 4,*f*_{B}= 2, and several values of the energy ratio*ɛ*_{AA}/*ɛ*_{BC}are considered.

In both models, particles and linkers self-assemble into structures formed by chains
of linkers (sequences of *AA* bonds) that may branch when bonded to
particles [through *AC* bonds in model A or *BC* bonds
in model B—see Figs. 2(c) and 2(d)]. However, the interplay between chaining
and branching is different: In model A, a patch A can bond to patches A or to
patches C, setting a competition between chaining and branching [see Fig. 2(c)]; in model B, each type of patch is
only engaged in one type of bond and this competition is absent [see Fig. 2(d)].

It is worth noticing that model A and model B (with *f*_{B} = 2) become the thoroughly
studied linker–particle model of Refs. 16 and 20 when *ɛ*_{AA} = 0. Therefore, the present
work can be seen as an extension of the study of linker–particle models to a case
where the self-assembly of linkers into chains is present.

### A. Simulations

Structural properties of the models were obtained as a function of temperature,
density, and composition with a classic Metropolis Monte Carlo (MC) simulation
in the canonical ensemble using Ref. 26 as reference. A mixture of 1200 particles and linkers in a cubic box with
periodic boundary conditions was set up. Particles and linkers were randomly
placed in the box with random orientations and then moved around until
equilibrium was reached. A move was defined as a simultaneous displacement
between ±0.05*σ* in each direction and a rotation between ±0.1
rad around a random axis, with all quantities being drawn from uniform
distributions. Simulations ran for a minimum of 10^{5} MC steps (each
step is defined as 50 000 attempts to move a particle or a linker). Steadiness
of the bonding probabilities (i.e., of the number of bonds formed in the system)
and the fraction of particles belonging to the largest cluster were used to
assess equilibrium. Cluster size distributions were obtained using the
Hoshen–Kopelman algorithm.^{27}

The particles were decorated with *f*_{C} = 4 patches of type C, placed on their surface as vertices of a tetrahedron. The
linkers were decorated with 2 patches of type A as (opposite) poles, and with *f*_{B} = 2 patches of type B,
placed at the equator and as opposite poles. The radius of the patches (or the
range of the patch–patch potential) was set to $\delta =(5\u221223\u22121)\sigma /2\u22480.119\sigma $,
the maximum value that guarantees that each patch is engaged in a single
bond.

### B. Theory

^{28–30}for mixtures with several types of bonds.

^{20,23}In this theory, closed loops are neglected and the clusters assume a tree-like bonding structure. As a consequence, the particles of a cluster can be separated by levels: A random particle is chosen as level 0; the particles bonded to this are at level 1, and so forth. This approach is briefly reviewed here for the case in which patches of a given type are only present in one of the species of the mixture, as happens in the models under study. The number of patches of type

*γ*that belong to particles of level

*i*and are bonded to particles of the previous level,

*b*

_{i,γ}, follows a recursive relation that can be expressed in a matrix form,

*b*

_{i,γ}. The matrix $T\u0303$ encodes the structure of the clusters and is a function of

*f*

_{γ}, the number of patches of type

*γ*in a particle, and of the probabilities

*p*

_{α→β}that are the fraction of patches of type

*α*that are bonded to patches of type

*β*. Percolation will be absent if the absolute values of all eigenvalues

*λ*of $T\u0303$ are lower than one. If at least one of the eigenvalues has absolute value larger than 1, then the number of bonds increases with increasing level and the system percolates. The percolation threshold is obtained, as a function of

*f*

_{γ}and of

*p*

_{α→β}, by finding the conditions for which 1 is the largest absolute value of all eigenvalues of matrix $T\u0303$.

*p*

_{α→β}are obtained as a function of density, composition, and temperature, through the laws of mass action of Wertheim’s theory,

^{23,30}which provide a connection between percolation and thermodynamics. Formation of bonds between patches of type

*α*and patches of type

*β*can be seen as a chemical reaction whose equilibrium is obtained when

*x*

_{β}is the fraction of particles that contain patches of type

*β*and

*p*

_{α}≡ ∑

_{γ}

*p*

_{α→γ}is the fraction of patches of type

*α*that are bonded. Δ

_{αβ}is interpreted as being the reaction constant for the formation of bonds (

*αβ*)

^{30}or the partition function of these bonds.

^{25}For the interaction potential (1), Δ

_{αβ}is obtained using Wertheim’s first order perturbation theory and a linear approximation for the HS fluid pair correlation function,

^{25}

*k*

_{B}is the Boltzmann constant,

*T*is the temperature,

*v*

_{αβ}is the bonding volume, i.e., the volume that can be explored by one particle when keeping the other particle fixed, without breaking the bond

*αβ*between the two. In most calculations,

*δ*

_{αβ}will be equal to the value used in simulations (i.e.,

*δ*

_{αβ}=

*δ*= 0.119

*σ*). The calculation of the probabilities

*p*

_{α→β}as a function of thermodynamic quantities is completed by using the normalization

*p*

_{α}= ∑

_{β}

*p*

_{α→β}in (3),

*p*

_{α}are obtained from the thermodynamic quantities, and the probabilities

*p*

_{α→β}can then be computed from (3) and introduced in matrix $T\u0303$ of (2). Finally, the eigenvalues

*λ*are obtained, and the percolation threshold is determined as a function of

*ρ*,

*T*, and

*x*for a given model.

#### 1. Model A

*A*and particles have

*f*

_{C}patches of type

*C*. Only bonds

*AA*and bonds

*AC*can form; as a consequence, matrix $T\u0303$ is (see Appendix A) given by

*λ*

_{+}, and one negative

*λ*

_{−}, with

*λ*

_{+}> |

*λ*

_{−}|. The equation for the percolation threshold is obtained by imposing

*λ*= 1 as a root of the characteristic polynomial of (8),

*p*

_{α→β}are related to the bonding probabilities

*p*

_{α}using (3) (recall that composition

*x*is the fraction of particles to which patches of type

*C*belong),

*AA*bonds were allowed

^{16,20}are recovered when

*p*

_{A→A}≡ 0.

#### 2. Model B

*A*and

*f*

_{B}patches of type B, while particles have

*f*

_{C}patches of type

*C*. Only bonds

*AA*and bonds

*BC*can form; as a consequence, matrix $T\u0303$ is (see Appendix A) given by

*λ*= 1 as a root of the characteristic polynomial of (15) (see Appendix B),

*p*

_{α→β}are related to the bonding probabilities

*p*

_{α}using (3),

## III. RESULTS

### A. Model A with *ɛ*_{AA}/*ɛ*_{AC} < 2

The percolation diagram of model A when linkers do not form chains (i.e., when *ɛ*_{AA} = 0) has already been
obtained^{16,20} and is
represented in Fig. 3 (solid line in all
panels). Percolation is present only for a limited range of compositions, 1/7 ≤ *x* ≤ 0.6, at low temperatures—Fig. 3(a). These two limiting compositions can be determined from (9) and (12) with *p*_{A→A} = 0
(and *f*_{C} = 4): when *p*_{C→A} → 1
(i.e., all patches C are bonded), *x* → 1/7; when *p*_{A→C} → 1
(i.e., all patches A are bonded), *x* → 0.6. This has a
transparent physical meaning: When chains of linkers are absent, the system
needs both a minimum amount of particles (1/7) and a minimum amount of linkers
(0.4) to form large clusters. Below *x* = 1/7, there are too many
linkers: Part of them bond to patches *C*, almost fully covering
the particles and preventing significant bonding between two particles (patches *C* are “blocked”); the other linkers are free and do not
contribute to clustering. Above *x* = 0.6, the linkers can
promote the bonding of particles, but their number is not enough to form large
clusters. The introduction of chains of linkers changes qualitatively the
temperature–composition percolation diagram for *x* ≤ 0.6—see Fig. 3(a). In particular, the limit *x* = 1/7 disappears and percolation can be obtained for all
compositions up to *x* = 0.6 at low temperatures. The cause of
this change is energy minimization. A bond *AC* decreases the
energy by *ɛ*_{AC}; the formation of a
bond *AA* from bonds *AC* requires the breaking of
two bonds *AC* (to free two patches *A*), and the
energy variation resulting from this process is
−*ɛ*_{AA} +
2*ɛ*_{AC}. As a consequence, at
low temperatures and when *ɛ*_{AA} <
2*ɛ*_{AC}, patches *A* will bond preferentially to patches *C*.
At low composition *x*, bonds *AC* will saturate
patches *C*, but all the remaining patches *A* can
now bond to form *AA* bonds. Therefore, the energy is minimized
by forming chains of linkers (*AA* bonds) that branch when
connected to particles (*AC* bonds). This structure is percolated
and will always form for 0 < *ɛ*_{AA}/*ɛ*_{AC} < 2 and *x* < 0.6 at sufficiently low temperatures.
Therefore, the self-assembly of linkers into chains opens the possibility of
percolation with fewer particles (for the same amount of linkers).

These theoretical predictions are confirmed with simulations, as shown in Fig. 4. For a series of equilibrated
simulations at fixed (*x*, *ρ*, *T*), the cluster size distribution was recorded at each
10^{3} MC steps, the fraction of particles and linkers that belong
to the largest cluster was determined, and its mean value (the points in Figs. 4 and 8) and standard deviation (the error bars in Figs. 4 and 8) were
calculated. The fraction of particles and linkers that belong to the largest
cluster, *p*_{∞}, is displayed, for *ɛ*_{AA}/*ɛ*_{AC} = 0 and *ɛ*_{AA}/*ɛ*_{AC} = 1 at a fixed low density as a function of *x* for a given low
temperature in Fig. 4(a) and as a function
of temperature for a given low *x* in Fig. 4(b). *p*_{∞} is the order parameter for the percolation transitions, going from 0 (in the
thermodynamic limit) below the percolation threshold to a nonzero value above
the percolation threshold. For the finite systems employed, the sudden increase
of *p*_{∞} signals the onset of
percolation. The simulation results of Fig. 4(a) confirm that the formation of chains of linkers decreases the
fraction of particles needed to obtain percolation: A strong variation of *p*_{∞} is observed at lower *x* for *ɛ*_{AA}/*ɛ*_{AC} = 1. On the other hand, the results depicted in Fig. 4(b) show that for a fixed, low *x*, percolation
is only obtained in the case *ɛ*_{AA}/*ɛ*_{AC} = 1.

The temperature at which, for low *x*, percolation disappears, has
a non-monotonic dependence on *ɛ*_{AA}/*ɛ*_{AC},
as can be seen in Figs. 3(a) and 3(c). This can be understood by recognizing
that there are three ways of disrupting the fully bonded percolated network that
is formed at low temperatures: (a) breaking the chains (or *AA* bonds) that connect the particles; (b) breaking the bonds between the chains and
the particles (or *AC* bonds); (c) a combination of (a) and (b).
These three regimes can be identified by determining the fraction of bonds *AA* and *AC* that do not break when the
system is heated from *T* = 0 to the percolation temperature.
These quantities [calculated at (*ρ*, *x*) = (0.1,
0.1)] are represented, as a function of *ɛ*_{AA}/*ɛ*_{AC},
in Fig. 5:

For low values of

*ɛ*_{AA}/*ɛ*_{AC}(up to 0.5), percolation disappears due to the breaking of chains: Around 60% of bonds*AA*are broken, while*AC*bonds remain unchanged. In this case, the energy cost of breaking*AA*bonds is low, and so the breaking of chains to an extent that leads to the disappearance of percolation may happen at a low temperature.For the larger values of

*ɛ*_{AA}/*ɛ*_{AC}(in the range $\u22481.5$ to 2), the network is disrupted due to the breaking of bonds between chains and particles: 2/3 of the bonds*AC*break, while new bonds*AA*are formed. In this regime, the energy cost of replacing two*AC*bonds by one*AA*bond and two free patches*C*is low, and therefore the percolation threshold occurs at low temperatures.For intermediate values of

*ɛ*_{AA}/*ɛ*_{AC}, percolation disappears through both mechanisms: When*ɛ*_{AA}=*ɛ*_{AC}, around 25% of bonds*AA*and bonds*AC*are broken at the percolation threshold. In this regime, the energy cost of breaking chains is larger than in (a) and the energy cost of replacing bonds*AC*by bonds*AA*and two free*C*patches is larger than in (c). As a consequence, the temperature at which percolation disappears is larger than in those regimes.

The dependence of the percolation threshold on density is illustrated in Figs. 3(b)–3(d). For a fixed intermediate temperature
[*k*_{B}*T*/*ɛ*_{AC} = 0.075 in Fig. 3(b)], percolation is
observed above a minimum density, for a range of compositions *x*. For values of *ɛ*_{AA}/*ɛ*_{AC} up to 1, this minimum density is low and does not change; the lower limit of the
range of compositions, on the other hand, decreases significantly when *ɛ*_{AA}/*ɛ*_{AC} ≈ 1. This effect is evident in Fig. 3(d):
The minimum number of particles needed to obtain percolation for a fixed number
of linkers decreases significantly for those values of *ɛ*_{AA}/*ɛ*_{AC}.
On the contrary, the maximum number of particles for which, for a given number
of linkers, percolation still exists is unaffected by *ɛ*_{AA}/*ɛ*_{AC}.
This behavior is totally different for *ɛ*_{AA}/*ɛ*_{AC} = 1.5: At *k*_{B}*T*/*ɛ*_{AC} = 0.075, percolation is obtained only at high densities [Figs. 3(b) and 3(d)];
percolation takes place at low densities only if the temperature is
significantly lowered—see Fig. 3(c).

### B. Model A with *ɛ*_{AA}/*ɛ*_{AC} > 2

The structure of the ground state of model A changes when *ɛ*_{AA} >
2*ɛ*_{AC}. It becomes energetically
more favorable to form bonds *AA* than bonds *AC*,
and all linkers tend to assemble into long chains without connecting to the
particles. Only chains with no branching are formed at low temperatures and so
percolation is not energetically favorable. We found no numerical solutions to
Eqs. (9)–(14) using *δ*_{AC} = *δ*_{AA} = *δ* =
0.119*σ* and *ɛ*_{AA} > 2*ɛ*_{AC}.

However, it is possible to envisage a mechanism that promotes percolation when it
is not energetically favorable to form branched structures. As temperature is
raised, bonds *AA* start to break and to free *A* patches that can bond to *C* patches. The resulting structure
will be determined by the competition between these two types of entropic
defects:^{3,31} free
patches *A* and bonds *AC*. Favoring bonds *AC* promotes branching and the possibility of percolation.
These bonds can be made entropically more favorable by increasing their volume *v*_{AC}. Figure 6 shows the percolation thresholds for *ɛ*_{AA}/*ɛ*_{AC} = 2.1 and values of *δ*_{AC}/*σ* = 0.45,
0.5, 0.6 (i.e., larger values of the bonding volume *v*_{AC}). *δ*_{AA} is kept equal to
0.119*σ*. The temperature–composition diagram of Fig. 6(a) exhibits a percolation loop that
appears only when *δ*_{AC} ≳ 0.4 and
that increases in size when *δ*_{AC} is
increased. This result shows that, at fixed low *ρ* and
intermediate *x*, the temperature increase favors, at first,
“entropic” branching to an extent that makes percolation possible. Further
increase in the temperature promotes breaking of both types of bonds and, at
some point, the vanishing of percolation. Therefore, the nonintuitive phenomena
of percolation (or gelation) by heating^{32} is obtained as a result of an entropic competition
in bond formation. The density of the percolation threshold is represented in Fig. 6(b) as a function of temperature,
at *x* = 0.1: There is a density below which percolation is not
possible; for larger densities, percolation occurs within a range of
temperatures.

The use of *δ*_{AC} >
0.119*σ* in the theoretical framework developed in Sec. II B is, strictly speaking, inconsistent.
Both the approximation for Δ_{AC} in (4) and the tree-like cluster
description employed for the structure assume that each patch can only be
engaged in one bond (single bond per patch condition). This hypothesis is
violated, for the interacting potential given by (1), when *δ*_{αβ} > 0.119*σ*. However, the results represented in Fig. 6 are still physically meaningful. In
fact, flexible linkers (or particles) and mobile patches on particles could give
origin to larger bonding volumes without compromising the single bond per patch
condition. Therefore, the results of Fig. 6 suggest that in systems where linear linker–linker and low valence
linker–particle self-assembly occur, entropic effects related to the flexibility
and mobility of linkers can give rise to unexpected features like percolation
(or gelation) by heating and closed loops in percolation diagrams.

We conclude that model A exhibits two distinct regimes for percolation that are
set solely by the relation between the energy scales *ɛ*_{AA} and *ɛ*_{AC}, irrespective of the
values of *δ*_{AC} and *δ*_{AA}: When *ɛ*_{AC} > *ɛ*_{AA}/2, percolated sates are
always obtained when *T* → 0 and *x* <
2(*f*_{C} −
1)/(3*f*_{C} − 2) (= 0.6 for *f*_{C} = 4); when *ɛ*_{AC} < *ɛ*_{AA}/2, percolated states
are absent when *T* → 0.

### C. Model B

Model B is characterized by the self-assembly of linkers into chains through *AA* bonds and by the branching of these chains when they
connect to the particles through *BC* bonds. Since each type of
patch is only involved in one type of bond, the formation of bonds *AA* does not constrain the formation of *BC* bonds (and vice versa). As a consequence, the structure of the ground state is
independent of *ɛ*_{AA}/*ɛ*_{BC} and corresponds to the maximization of the number of both types of bonds.

The percolation diagrams of Fig. 7 (for *f*_{B} = 2) demonstrate that the
introduction of chain assembly in model B leads to substantial changes. The
temperature–composition diagram of Fig. 7(a) shows that percolation is obtained for all compositions *x* at sufficiently low temperatures. This is a drastic
change from the case *ɛ*_{AA} = 0 and *f*_{C} = 4, where percolation could
only be found for 1/7 < *x* < 0.6, and from model *A*, where percolation was never found for *x* > 0.6. The presence of percolation at extremely low *x* and at
low temperatures means that a vanishingly small amount of particles is enough to
connect, through *BC* bonds, the abundant and long chains of
linkers formed by *AA* bonds, to an extent that leads to the
formation of percolating clusters. On the other hand, the emergence of
percolation at *x* close to 1 and low temperatures means that a
vanishingly small amount of linkers is enough to connect the particles in a
percolated cluster, through the combination of a few bonds *BC* that connect linkers to particles and of bonds *AA* that connect
linkers to linkers and, indirectly, particles to particles.

The density dependence of percolation is shown in Figs. 7(b)–7(d). At a given
temperature
[*k*_{B}*T*/*ɛ*_{BC} = 0.075 in Fig. 7(b)], percolation only
exists above a certain density, and this minimum density decreases with *ɛ*_{AA}/*ɛ*_{BC}.
Above this density, percolation occurs for a range of compositions that
increases with increasing *ɛ*_{AA}/*ɛ*_{BC}.
In practice, when *ɛ*_{AA}/*ɛ*_{BC} is high, percolation occurs for every composition. At fixed composition
[*x* = 0.1 in Fig. 7(c)], the temperature at which percolation occurs for a given density
increases with *ɛ*_{AA}/*ɛ*_{BC}.
At a fixed temperature, the density of particles (linkers) needed to obtain
percolation for a given density of linkers (particles) decreases significantly
with increasing *ɛ*_{AA}/*ɛ*_{BC}—see Fig. 7(d). The results of Fig. 7 show that, in model B, when chaining
of linkers is favored, percolation is controlled mainly by temperature (see the
results for *ɛ*_{AA}/*ɛ*_{BC} = 1.5): Below a threshold temperature (almost constant), percolation occurs for
all compositions and all densities above a very small minimum density.

Some predictions of the theory for model B were tested (and confirmed) by
simulations. Their results are shown in Fig. 8. The introduction of chains of linkers reduces the fraction of
particles at which, at fixed low density and low temperature, percolation
occurs, as shown by the comparison between the results for *ɛ*_{AA}/*ɛ*_{BC} = 0.8 and for *ɛ*_{AA}/*ɛ*_{BC} = 0 in Fig. 8(a). The results shown in Fig. 8(b) confirm the increase of the
temperature at which percolation occurs (for fixed low density and composition)
when chains of linkers are present.

### D. Structural properties at the percolation threshold

The theory presented in Sec. II B can be
developed to obtain some properties of the clusters, within its assumption that
clusters are tree-like. In particular, the percolation threshold Eqs. (9) and (16) can be expressed in terms of
two structural properties: ⟨*n*_{C}⟩,
the mean number of bonds formed by one particle, which quantifies the number of
chains that branch in one particle; and
⟨*ℓ*_{L}⟩, the mean length of
the chains formed by linkers, which quantifies the extent of chaining.

*n*

_{C}⟩ is calculated from the number of patches

*f*

_{C}of the particles and the probability

*p*

_{C}that a patch

*C*is bonded,

*AA*bonds. The number of linkers in sequences with

*n*linkers and

*n*− 1

*AA*bonds is proportional to $pA\u2192An\u22121$.

*n*

_{c}⟩ and ⟨

*ℓ*

_{L}⟩, using (22), (23), and (9) (for model A) or (16) (for model B), to obtain, for model A,

*n*

_{C}⟩ and ⟨

*ℓ*

_{L}⟩ at percolation depends only on composition and not on temperature, density, or the ratio of energy scales. It is possible to obtain percolated clusters with lower values of ⟨

*n*

_{C}⟩ by increasing the length of the chains of linkers. An increase (decrease) in chaining leads to a decrease (increase) in branching. In the original linker–particle model where chains of linkers are absent, it is only possible to control one feature of the onset of percolation since ⟨

*n*

_{C}⟩ at the percolation threshold is defined by composition $(\u27e8nC\u27e92=2fC(1\u2212x)/((fC\u22121)x))$. The models under study have another degree of control: It is possible to create percolated states with more or less branching by changing the propensity to chaining (and vice versa).

The percolation thresholds expressed by Eqs. (24) and (25) can be used to obtain
(⟨*n*_{c}⟩,
⟨*ℓ*_{L}⟩) percolation diagrams. Figure 9 depicts the four
representative diagrams, each obtained at a fixed composition *x*. The thick lines represent the physical limits imposed by the
restrictions: (i) 0 ≤ *p*_{A→A} ≤ 1,
which translates into ⟨*ℓ*_{L}⟩≥ 1; (ii)
0 ≤ *p*_{C} ≤ 1, which is equivalent to
0 ≤ ⟨*n*_{C}⟩ ≤ *f*_{C}; (iii) for model A, *p*_{A→A} + *p*_{A→C} ≤ 1,
which is equivalent to
⟨*n*_{C}⟩⟨*ℓ*_{L}⟩
≤ 2(1 − *x*)/*x*, and, for model B, *p*_{B→C} ≤ 1,
which is equivalent to ⟨*n*_{C}⟩ < *f*_{B}(1 − *x*)/*x*. In model *A* [see Figs. 9(a) and 9(b)], percolation only occurs for a limited range of values
(⟨*n*_{c}⟩,
⟨*ℓ*_{L}⟩):
⟨*n*_{C}⟩ ≥ *f*_{C}/(*f*_{C} − 1) and $\u27e8\u2113L\u27e9\u2264\u27e8\u2113L\u27e9max$,
with $\u27e8\u2113L\u27e9max=2(fC\u22121)(1\u2212x)/(fCx)$. Therefore, larger chains
of linkers in a percolated network can only be obtained if composition *x* is lowered. The percolation diagrams for model A display
another dependence on composition: For
2/(*f*_{C} + 2) ≤ *x* ≤
2(*f*_{C} – 1)/(3*f*_{C} – 2),
the percolated states exhibit limited branching since
⟨*n*_{C}⟩ < *f*_{C} [see Fig. 9(b)]; for *x* >
2(*f*_{C} −
1)/(3*f*_{C} − 2), percolation
disappears, as already seen in Secs. III A and III B. The percolation behavior of
model B is different—see Figs. 9(c) and 9(d): For all compositions, percolated
networks with all possible lengths of chains of linkers can be obtained. The
composition only restricts branching in percolated states: If *x* > *f*_{B}/(*f*_{B} + *f*_{C}), then
⟨*n*_{C}⟩ < *f*_{B}(1 − *x*)/*x* [see Fig. 9(d)].

## IV. CONCLUSIONS AND DISCUSSION

The present work shows that in linker-mediated aggregation of particles, the
introduction of self-assembled chains of linkers gives an extra control over
percolation, compared to the case where chains of linkers are absent. The choice of
the specific interactions between linkers (model A or model B), the strength of
those interactions (value of *ɛ*_{AA}/*ɛ*_{AC} or *ɛ*_{AA}/*ɛ*_{BC}),
and the bonding volumes can be used to tune the composition, temperature, and
density at which percolation occurs and to change the structure of the percolated
clusters.

In model B, where chaining of linkers and branching through bonds between chains and
particles are independent, it is possible to reach percolation at both low and high
concentrations of particles: A vanishingly small number of particles (linkers) is
enough for percolation to occur in systems with a large fraction of linkers
(particles). Model A—where competition between chaining of linkers and branching of
chains (through their bonds to particles) is present—exhibits a more complex
behavior, which depends on the ground state of the model. For *ɛ*_{AA}/*ɛ*_{AC} < 2, branching is energetically favored, and percolation is extended to low
particle composition. The temperature at which, at low compositions, percolation
disappears depends non-monotonically on *ɛ*_{AA}/*ɛ*_{AC},
being maximum for *ɛ*_{AA} ≈ *ɛ*_{AC}. This non-monotonic behavior is
related to the energy costs of disrupting the percolated network, either by breaking
the chains or by breaking the bonds between chains and particles. Percolation occurs
in model A with *ɛ*_{AA} >
2*ɛ*_{AC} only when the entropic gain of
branching is increased. In this case, percolation by heating and closed loops of
percolated phases in temperature–composition diagrams are obtained. The number of
branching bonds per particle is related to the mean length of chains of linkers
through a function that depends only on composition. This means that composition
controls the structure of the percolated clusters: It is possible to obtain clusters
with larger (smaller) chains at the cost of decreasing (increasing) branching. Table I shows a summary of the main
results.

. | Model A
(ɛ_{AC} ≠ 0)
. | Model B
(ɛ_{BC} ≠ 0)
. | ||||
---|---|---|---|---|---|---|

. | . | . | ɛ_{AC} < ɛ_{AA}/2
. | . | . | |

. | ɛ_{AA} = 0
. | ɛ_{AC} > ɛ_{AA}/2
. | Low δ_{AC}
. | High δ_{AC}
. | ɛ_{AA} = 0 f_{B} = 2
. | ɛ_{AA} ≠ 0
. |

Branching | Chaining | Chaining energetically | Branching | |||

energetically | energetically | favorable; branching | energetically | |||

Structural properties | No chaining | favorable | favorable | entropically favorable | No chaining | favorable |

Percolation loop | ||||||

Composition, x | x_{min} ≤ x ≤ x_{max} | x ≤ x_{max} | No percolation | in (x, T) | x_{min} ≤ x ≤ x_{max} | 0 ≤ x ≤ 1 |

Percolation loop | ||||||

Temperature, T | Low T | Low T | No percolation | in (x, T) | Low T | Low T |

. | Model A
(ɛ_{AC} ≠ 0)
. | Model B
(ɛ_{BC} ≠ 0)
. | ||||
---|---|---|---|---|---|---|

. | . | . | ɛ_{AC} < ɛ_{AA}/2
. | . | . | |

. | ɛ_{AA} = 0
. | ɛ_{AC} > ɛ_{AA}/2
. | Low δ_{AC}
. | High δ_{AC}
. | ɛ_{AA} = 0 f_{B} = 2
. | ɛ_{AA} ≠ 0
. |

Branching | Chaining | Chaining energetically | Branching | |||

energetically | energetically | favorable; branching | energetically | |||

Structural properties | No chaining | favorable | favorable | entropically favorable | No chaining | favorable |

Percolation loop | ||||||

Composition, x | x_{min} ≤ x ≤ x_{max} | x ≤ x_{max} | No percolation | in (x, T) | x_{min} ≤ x ≤ x_{max} | 0 ≤ x ≤ 1 |

Percolation loop | ||||||

Temperature, T | Low T | Low T | No percolation | in (x, T) | Low T | Low T |

The simulations performed for model A (for *ɛ*_{AA} <
2*ɛ*_{AC}) and for model B have
confirmed the qualitative differences predicted by the theory for percolation
diagrams obtained when chains of linkers are present. Figures 4 and 8(a) show, through
the size of the error bars, that the cluster size distribution displays much larger
fluctuations in a broader region of composition or temperature when linkers are
dominant (*x* ≲ 0.2) and can bond to each other
(*ɛ*_{AA} ≠ 0). Therefore, the
introduction of self-assembled chains of linkers will require larger systems and
longer simulations if more precise measurements of percolation thresholds are
desired. Moreover, the quantitative comparison between simulation and theoretical
results for model B would have to take into account a fairly strong tendency for
loop formation, which we have checked in a more detailed analysis of clusters in the
simulations of this model.

Percolation is only a necessary condition to obtain an equilibrium gel. The phase
diagrams of the models of this work are needed to know the thermodynamic conditions
at which percolation occurs in single phase regions. These calculations are left for
future work, but some hints can be found in previous studies.^{20,22,23}

The case of model A (or model B with *f*_{B} = 2) with *ɛ*_{AA} = 0 (particle–linker
mixture with no chains) and *f*_{C} = 3 was
investigated in Ref. 20. Percolated liquids,
not affected by phase separation, were found at intermediate compositions, low
densities, and low temperatures (see Figs. 4–6 in Ref. 20). Moreover, in Ref. 23, the phase diagram of model A with *ɛ*_{AA} = *ɛ*_{AC} and *f*_{C} = 3 shows that at high
linker concentrations, the percolated liquids have less tendency to phase separate
with a vapor (see Fig. 6 in Ref. 23). Even in
a mixture of patchy particles with 2 and 4 patches, all equal and all forming bonds,
large regions of a single phase percolated liquid are found.^{22} Therefore, it is much likely that the models
studied in this work will exhibit regions of phase diagrams with a single percolated
liquid phase.

Finally, it should be mentioned that, to complete the study of equilibrium gels, the
onset of mechanical stability in the percolated network has to be addressed.^{33,34} The emergence of rigidity
in equilibrium self-assembling thermodynamic systems is poorly understood.^{35} However, recent developments of a
theory that describes the mechanical properties of ideal reversible polymer networks
at thermodynamic equilibrium may apply to linker–particle systems,^{36} if properly generalized to the
case of a mixture.

## ACKNOWLEDGMENTS

We acknowledge financial support from the Portuguese Foundation for Science and Technology (FCT) under Contract Nos. PTDC/FIS-MAC/28146/2017, LISBOA-01-0145-FEDER-028146, PTDC/FIS-MAC/5689/2020, CEECIND/00586/2017, UIDB/00618/2020, and UIDP/00618/2020.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**M. Gouveia**: Investigation (equal); Software (equal); Visualization
(equal); Writing – review & editing (equal). **C. S. Dias**:
Conceptualization (equal); Investigation (equal); Supervision (equal); Writing –
review & editing (equal). **J. M. Tavares**: Conceptualization
(equal); Investigation (equal); Supervision (lead); Writing – original draft
(lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: MATRIX $T\u0303$

Let us consider the levels of a tree-like cluster (with no loops) formed by
linkers and particles connected through the bonding of its patches as described
in Sec. II B. *b*_{i,γ} is
the number of patches of type *γ* that belong to particles or
linkers in level *i* and that are bonded to the previous level.
Matrix $T\u0303$ in (8) and (15) is obtained from the relation
between the number of patches *b*_{i,γ} of
different levels. Figure 10 represents,
for models *A* and *B*, all the possible ways of
obtaining a patch of level *i* + 1 that is bonded to level *i*, from the patches of level *i* that are
bonded to level *i* − 1.

#### 1. Model A

*AA*and

*AC*are allowed. A patch

*C*at level

*i*+ 1 can only be obtained, with probability

*p*

_{A→C}, from a bond that originates in a patch

*A*of a linker of level

*i*; this linker is connected to level

*i*− 1 through its other patch

*A*. Therefore,

*i*+ 1 can be obtained in two ways: (i) with probability

*p*

_{A→A}, from a bond that originates in a patch A of a linker of level

*i*—this linker is connected to level

*i*− 1 through its other patch A; (ii) with probability

*p*

_{C→A}, from a bond that originates in one of the available

*f*

_{C}− 1 patches of a particle of level

*i*—this particle is connected to level

*i*− 1 by its remaining C patch. As a consequence,

#### 2. Model B

In model B, bonds *AA* and *BC* are allowed and
3 different types of patches are involved in bonds.

*i*+ 1 can be obtained, with probability

*p*

_{C→B}, from a bond that originates in one of the available

*f*

_{C}− 1 patches of a particle of level

*i*; this particle is connected to level

*i*− 1 by its remaining C patch. As a consequence,

*i*+ 1 can be obtained, with probability

*p*

_{A→A}, from a bond that originates in a patch A of a linker of level

*i*; this linker can be bonded to level

*i*− 1 by a patch

*A*or by a patch

*B*(in which case, the 2 patches

*A*of the linker of level

*i*are available to bond to the patch

*A*of level

*i*+ 1). As a consequence,

*i*+ 1 can be obtained, with probability

*p*

_{B→C}, from a bond that originates in a patch B of a linker of level

*i*; this linker can be bonded to level

*i*− 1 by (i) a patch

*A*, in which case the

*f*

_{B}patches B of the linker of level

*i*are available to bond to the patch C of level

*i*+ 1, or (ii) by a patch

*B*, in which case only

*f*

_{B}− 1 patches of the linker of level

*i*are available to bond to the patch

*C*of level

*i*+ 1. As a consequence,

### APPENDIX B: PERCOLATION THRESHOLD OF MODEL B

*λ*= 1 is an eigenvalue of (15), as was done in (16) and in subsequent calculations. The eigenvalues

*λ*of (15) are the solutions of

*z*= (

*f*

_{C}− 1)

*p*

_{B→C}

*p*

_{C→B}and

*p*

_{A}=

*p*

_{A→A}.

We will first show that if

*λ*= 1 is a solution of (B1), this solution is the one with the largest absolute value.If*λ*= 1 is a solution of (B1), then (B1) can be rewritten aswhere(B2)$(\lambda \u22121)(\lambda 2+b\lambda +c)=0,$*b*= 1 −*p*_{A},*b*−*c*= (*f*_{B}− 1)*z*, and*c*= (*f*_{B}+ 1)*p*_{A}*z*. From these equalities,*c*can be expressed as a function of*p*_{A},The other eigenvalues,(B3)$c=(fB+1)(1\u2212pA)pAfB\u22121+(fB+1)pA.$*λ*_{±}, areIf these eigenvalues are not real (i.e., if(B4)$\lambda \xb1=\u2212b\xb1b2\u22124c2.$*b*^{2}− 4*c*< 0), then*λ*= 1 is the only real solution of (B1) and, therefore, is the one with the largest absolute value. If they are real, one can use the relation between*c*,*b*, and*p*_{A}and, taking into account that 0 <*p*_{A}< 1, obtain their absolute values,(B5)$|\lambda +|=1\u2212pA21\u22121\u2212G(fB,pA),$where(B6)$|\lambda \u2212|=1\u2212pA21+1\u2212G(fB,pA),$Since $1\u2212pA2<1/2$, $1\u22121\u2212G(fB,pA)<1$, and $1+1\u2212G(fB,pA)<2$, we conclude that |(B7)$G(fB,pA)=4(fB+1)pA(1\u2212pA)fB\u22121+(fB+1)pA.$*λ*_{+}| < 1 and |*λ*_{−}| < 1. Therefore, if there is an eigenvalue*λ*= 1, the other two, if they exist, have absolute values lower than 1.We will show that if

*λ*= −1 is a solution of (B1), then there is another eigenvalue whose absolute value is larger than 1 and, therefore,*λ*= −1 does not define a percolation threshold.Assume that one of the solutions of (B1) is*λ*= −1, so that (B1) can then be rewritten aswhere(B8)$(\lambda +1)(\lambda 2+b\lambda +c)=0,$*b*= −1 −*p*_{A},*b*+*c*= −(*f*_{B}− 1)*z*, and*c*= −(*f*_{B}+ 1)*p*_{A}*z*. From these equalities,*c*and*z*can be expressed as a function of*p*_{A}(B9)$c=\u2212(fB+1)(1+pA)pAfB\u22121\u2212(fB+1)pA,$Since(B10)$z=1+pAfB\u22121\u2212(fB+1)pA.$*z*> 0, (B10) is only valid if $pA<fB\u22121fB+1$, which means that*λ*= −1 can be a solution of (B1) only when*p*_{A}satisfies this inequality. The other eigenvalues,*λ*_{±}, arewhere(B11)$\lambda \xb1=\u2212b\xb1b2\u22124c2=(1+pA)1\xb11+F(pA,fB)2,$For $pA<fB\u22121fB+1$ [i.e., when(B12)$F(fB,pA)=4(fB+1)pA(1+pA)fB\u22121\u2212(fB+1)pA.$*λ*= −1 can be one of the solutions of (B1)],*F*(*f*_{B},*p*_{B}) > 0 and so*λ*_{±}are always real. The absolute value of*λ*_{+}isSince 1 +(B13)$|\lambda +|=(1+pA)1+1+F(pA,fB)2.$*p*_{A}> 1 and $1+1+F(pA,fB)>2$, |*λ*_{+}| > 1. As a consequence,*λ*= −1 is never the eigenvalue of (15) with maximum absolute value and cannot be used to obtain the percolation threshold.

We conclude that the general definition for the percolation threshold, that is, the thermodynamic conditions for which the maximum absolute value of all eigenvalue of matrix $T\u0303$ is 1, reduces, for model B, to the conditions at which one of the eigenvalues is equal to 1.