Extended Wertheim theory predicts the anomalous chain length distributions of divalent patchy particles under extreme confinement

Colloidal patchy particles with divalent attractive interaction can self-assemble into linear polymer chains. Their equilibrium properties in 2D and 3D are well described by Wertheim's thermodynamic perturbation theory which predicts a well-defined exponentially decaying equilibrium chain length distribution. In experimental realizations, due to gravity, particles sediment to the bottom of the suspension forming a monolayer of particles with a gravitational height smaller than the particle diameter. In accordance with experiments, an anomalously high monomer concentration is observed in simulations which is not well understood. To account for this observation, we interpret the polymerization as taking place in a highly confined quasi-2D plane and extend the Wertheim thermodynamic perturbation theory by defining addition reactions constants as functions of the chain length. We derive the theory, test it on simple square well potentials, and apply it to the experimental case of synthetic colloidal patchy particles immersed in a binary liquid mixture that are described by an accurate effective critical Casimir patchy particle potential. The important interaction parameters entering the theory are explicitly computed using the integral method in combination with Monte Carlo sampling. Without any adjustable parameter, the predictions of the chain length distribution are in excellent agreement with explicit simulations of self-assembling particles. We discuss generality of the approach, and its application range.

Colloidal patchy particles with divalent attractive interaction can self-assemble into linear polymer chains.Their equilibrium properties in 2D and 3D are well described by Wertheim's thermodynamic perturbation theory which predicts a well-defined exponentially decaying equilibrium chain length distribution.In experimental realizations, due to gravity, particles sediment to the bottom of the suspension forming a monolayer of particles with a gravitational height smaller than the particle diameter.In accordance with experiments, an anomalously high monomer concentration is observed in simulations which is not well understood.To account for this observation, we interpret the polymerization as taking place in a highly confined quasi-2D plane and extend the Wertheim thermodynamic perturbation theory by defining addition reactions constants as functions of the chain length.We derive the theory, test it on simple square well potentials, and apply it to the experimental case of synthetic colloidal patchy particles immersed in a binary liquid mixture that are described by an accurate effective critical Casimir patchy particle potential.The important interaction parameters entering the theory are explicitly computed using the integral method in combination with Monte Carlo sampling.Without any adjustable parameter, the predictions of the chain length distribution are in excellent agreement with explicit simulations of self-assembling particles.We discuss generality of the approach, and its application range.
Synthetic colloidal particles suspended in a nearcritical binary liquid mixture (e.g.water and lutidine), attract each other via a solvent mediated critical Casimir force.Through novel synthesis routes these particles can be designed such that they form directed bonds between patches on the surface of neighboring particles 1 .As such patchy particles simultaneously experience thermal motion, their statistical behavior follows the Boltzmann distribution.Hence, they can be viewed as mesoscopic analogs of (carbon) atoms, which can be directly observed via, e.g., confocal microscopy 2 .In this way, patchy particles can act as an experimental model system to explore complex self-assembled structures analogous to molecular architectures, such as chains, rings, and networks [2][3][4][5][6] .
To understand the self-assembly in patchy particle systems, one can of course resort to computer simulations [7][8][9][10][11] , but an attractive alternative is to invoke statistical mechanics which aids to a better theoretical understanding and prediction.One of the classical theories for self-assembly of colloidal particles is the Wertheim thermodynamic perturbation theory (TPT) [12][13][14][15] , later reformulated as Statistical Associating Fluid Theory (SAFT) by Chapman et al. 16 .Wertheim's theory was originally intended as a molecular model 17 , but also works for mesoscopic particles.For divalent patchy particles in two and three dimensions, Wertheim theory is able to predict the polymerization equilibrium in terms of a) Electronic mail: p.g.bolhuis@uva.nlfor example the chain length distribution, with the total particle density and a pair bonding strength as the only input parameters [18][19][20][21][22][23] .For systems with average valencies larger than two, equilibrium properties are predicted using Flory-Stockmayer's polymer theory 24,25 .The location of the percolation point, existence of empty liquids and equilibrium gels were predicted theoretically, confirmed in simulation and validated experimentally [26][27][28][29][30][31] .
However, when there is a mismatch in density of the particles and the suspending solvent, particles will sediment to the bottom of the sample due to gravity.For sufficiently low particle concentration (or volume fraction) and short gravitational height, the system is then confined to a quasi-2D plane; making single layer structures possible.Direct application of Wertheim's theory for divalent particles in 2D or 3D will give an exponential distribution and a large discrepancy between the experiments and theoretical prediction exists, in particular in the monomer (and dimer) density 6,32 .In this work we address this discrepancy.
The origin of the discrepancy is that, under extreme confinement where spherical particles live in a twodimensional (x, y)−plane, the monomers are still able to rotate around their center-of-mass (Fig. 1a).While some monomer orientations, e.g. when their patch points toward the wall, are part of the orientational phase space, their patches are not available for bonding.This renders the system fundamentally different compared to the standard 2D and 3D systems.As a consequence, in order to predict thermodynamic properties, this excess rotational degree of freedom needs to be taking into account 33 .For strong confinement due to gravity, there is an additional anisotropy in the density along the direction perpendicular to the wall.Monomers and small chains still having freedom to translate and rotate against gravity, while for long chains only a small part of the chain (at the end) has this freedom as illustrated in Fig. 1b.
Note that directional assembly under extreme confinement not only occurs in model colloidal patchy particle systems.In chemistry, there are many examples where confinement has been used at its advantage.For example, nanoporous materials with pore shapes and sizes comparable to the typical size of small molecules such as metal-organic-frameworks (MOFs), covalent organic frameworks (COFs), zeolitic imidazolate frameworks (ZIFs) 34,35 , and nanopores composed of for example carbon nanotubes 36 .Another example is that of self-assembled supramolecular structures where the intermolecular non-covalent bonds determine the structure and chemical function, and confinement affects the reactivity 37 .Even in the confined environment of the living cell, where short-ranged, strongly directional hydrogen bonds provide a mode for molecular assembly.Recently, the fabrication of nanoslits with Ångström-scale separation became possible opening an exciting field of nanofluidics that show unusual dynamics, kinetics, and thermodynamics due to the extreme confinement [38][39][40][41] .However, a thorough theoretical understanding of the effect of confinement on for example separation and phase transitions is still lacking 42,43 .
Such highly confined systems could be in principle described theoretically with TPT [44][45][46] .When solving TPT, one has to compute the interaction parameters in the theory.There are two major routes to do this: via an 'integral method' 47 or via classical density functional theory (DFT) [48][49][50] .Solving TPT becomes increasingly complex for inhomogeneous systems due to positional and orientational coupling 51 .While extremely powerful, the DFT-route is currently not able to predict thermodynamic equilibrium for the Wertheim theory in highly confined systems at low temperatures quantitatively.See e.g.Ref. 52, which shows excellent predicted density distributions for tetrapatch particles at large wall separation, but for small wall separations (between 1.18 − 3.02 times the particle diameter) shows discrepancies for three different associating density functionals.In contrast, Ref. 53, following the integral route, predicted densities accurately of spherical dipatch particles in a one-dimensional pore with a width of the particle diameter.
In this paper we take a different approach.We interpret the rotational and translational freedom against the gravitational field as an additional source of entropy.This orientation-positions-dependent entropy effectively reduces the reactivity.We separate the polymerization reactions of the species which gives rise to adapted expressions for the chain length distributions in the Wertheim theory.By computing the interaction parameter via the integral method with Monte Carlo (MC) sampling, we can directly calculate the excess rotational and translational entropy and capture the corresponding equilibrium reaction constants between the species.
To validate our approach, we simulate patchy particles interacting via a simple square well potential combined with different forms of orientation-dependent switching functions under various gravitational strengths.Additionally, we apply the theory to the critical Casimir dipatch colloid particle system, for which the chain length distribution was experimentally studied in Ref. 6 and an accurate effective potential model was developed recently in Ref. 32.The extended Wertheim theory contains no fit parameters, and needs only two input parameters: the species-dependent interaction parameters and the total particle number density ρ.We compare the predicted distribution with the simulated ones, and find excellent quantitative agreement.
The paper is organized as follows: we start with a brief overview of the traditional Wertheim theory that holds both for 2D and 3D.Then, we will introduce the adapted Wertheim theory for the highly confined system in quasi-2D followed by the gravitationally confined systems.Using the quasi-2D systems, we show how we can determine the excess rotational free energy of the monomers and its effect on the chain length distribution.Next, we introduce the external gravitational field which gives also short chains additional entropy and thus higher probability of occurrence and show that the flexibility of the chain also plays a role on the distributions.This effect too can be determined via the integral method giving excellent predictions of the chain length distributions of divalent colloidal particles.Finally, we apply and validate the theory on our accurate patchy particle model interacting via critical Casimir interactions under realistic gravitational conditions.We end with concluding remarks, and a future outlook.

A. First order thermodynamic perturbation theory
Consider a 3D suspension of hard spherical particles or a 2D suspension of hard disks, in which each particle is divalent, i.e dressed with two attractive patches, usually located at opposing poles.Each patch or site is able to make a bond with a site on another particle, resulting in the association of particles into linear chains.Moreover, each site is able to make only one single bond, and each bond is equally likely to form.The aggregation of the monomers into larger clusters can then be viewed as a set of addition reactions: where A 1 stands for a monomer, A 2 for a dimer, and A n for a chain composed of n monomers.This type of re-< l a t e x i t s h a 1 _ b a s e 6 4 = " l 9     activity is also known as isodesmic polymerization where each addition reaction of a monomer is associated with equal amount of free energy 54 .All reactions have an equilibrium constant K, defined through the law of mass action: where [A n ] denotes the concentration of A n , and ρ n the number density of n-mers, i.e., chains of length n.
For this situation one can apply Wertheim 's first-order thermodynamic perturbation theory (TPT1), and calculate the chain length probability distribution 12,18 .The probability of observing a non-occupied binding site is denoted as X.Then the number fraction of chains of size n is on average given by with ρ the number density.This expression is rationalised as follows.A cluster of size n has (n − 1) links.
The probability of forming a link is 1 − X.The probability of forming n − 1 links is thus (1 − X) n−1 .However, as there are two unoccupied reactive sites, that accounts for a factor X 2 .So, for a monomer this reduces to ρ 1 = ρX 2 .
To convince oneself that this is consistent, one can add up all chain lengths, which would have to add up to the total density of particles: The geometric sum adds up to 1/X 2 , which is indeed consistent.
Next, using Eq. 3 and ρ 1 = ρX 2 , the equilibrium reaction constant K from Eq. is rewritten as: Solving for X gives Thus, given a density ρ as well as an equilibrium constant K, Eq. 3 together with Eq. 6 form a complete description of the system.The slope of (the log of) the chain length distribution is: where the latter equality defines the probability X b for binding or, equivalently, the fraction of bound sites X b .
Traditionally, TPT computes the equilibrium constant K via the interaction parameter ∆ 47 , as K ≡ M ∆, with M the number of binding sites or patches per particle.The ∆-parameter represents the (exponential of the) free energy difference of the bonding reaction with respect to the hard particle reference state.It is calculated via an integration over space of the Boltzmann weighted energy averaged over the allowed orientations of the particles, multiplied by the probability of finding a particle at distance r, i.e. the radial distribution function g(r).The interaction parameter ∆ is then: where r is the inter-particle vector of particle α and γ with their orientations Ω α and Ω γ , respectively, g(r) is the pair correlation function of the reference systems e.g.hard spheres, f (r, Ω α , Ω γ ) = exp −βU (r,Ωα,Ωγ ) −1 is the Mayer function, and, finally, Ωα,Ωγ denotes the orientational average of the Mayer function of particle α and γ separated at distance |r| = r.This calculation of ∆ is non-trivial, and depends on the geometry of the setup.
< l a t e x i t s h a 1 _ b a s e 6 4 = " w / 0 o 8 M 0 2. Schematic illustration of the two differences in sampling K2 and K in the left and right picture, respectively.The freely rotating α-particle (with dotted curved arrows indicate rotation axes) is radially sampled uniformly in the integration area V (blue area) around the γ-particle which is positioned at (0, 0).Due to symmetry we can reduce the sampling volume to a hemispherical shell around the γ-particle.

B. Wertheim theory in quasi-2D
The treatment in Sec I A assumes that all association reactions follow identical statistics.However, the situation is slightly different when confining the chain formation to a plane by e.g. two walls with particle diameter separation (Fig. 1A).In principle, the above described TPT/SAFT framework also applies in that case, except for one crucial difference in the assumption about the reactivity.Upon binding of single particles (free monomers) an excess rotational entropy is lost.Therefore, the first reaction in the series where two free monomers react to form a dimer (Eq.9a) is fundamentally different from the others where only one free monomer reacts with an existing cluster of size n > 1 (Eq.9b-9c): which leads to an increased monomer concentration as observed in the chain length distribution.
By classifying not just one type of bonding reaction, but two types of reactions (see Fig. 2), with corresponding free energies and equilibrium constants, we derive the extended Wertheim theory in quasi-2D.For the two addition reactions, the constants K 2 and K are given by Rewriting the first equation gives while the next addition reaction yields Continuing along this line, for n ≥ 2 it follows that Note that the above equations express the number densities of the clusters.So adding all densities multiplied by the chain lengths will give the imposed total particle density ρ Taking out a factor K 2 /K from the sum leads to To make this a tractable sum, we add and subtract a term K2 K ρ 1 , yielding When K 2 = K, the first term vanishes on the rhs of Eq. 16, recovering Eq. 4 in the original TPT1.In this we recover the fraction/probability of bound sites X b = 1 − X = Kρ 1 (see Eq. 7).
The geometric sum can now be evaluated, giving Next we divide by ρ, and define the monomer fraction Applying TPT1 amounts to solving Eq. 18 for the unknown monomer fraction X 1 , which in turn sets the entire distribution ρ n in Eq. 16.This requires knowledge of the equilibrium constants K and K 2 that are directly proportional to ∆ in Eq. 8, and follow from evaluating the integral.The difference between the calculations of the ∆'s corresponding to K and K 2 will be explained later in Sec.II D.
For the chain length probabilities P n , the chain densities are normalized with the chain number density ρ c ≡ n ρ n , which, using Eq. 13, results, analogous to Eq. 17, in so that: In this way, also the average chain length follows:   x < l a t e x i t s h a 1 _ b a s e 6 4 = "

C. Wertheim in a gravitational field
The above description for quasi-2D confinement also holds for an infinitely short gravitational height, in which translation away from the confining wall is strongly suppressed.If the gravitational field is not so strong, the particles are able to levitate to the order of the gravitational height.In turn, this translational freedom affects both the reactivity of particle association, and the free energy.
As a consequence, there may now be multiple reaction constants: where the reaction constants K 2 = K 3 = K k = K may not be equal with each other (see Fig. 3 for an illustration).Only beyond a certain chain length k (Eq.21c-21d), the equilibrium constant can be considered to settle.
The density ρ is written as outlined in the previous section: The second sum in this expression converges to: Similarly, the total number of clusters/chains in the system is: where the infinite sum converges to: Again, for given density ρ and reaction constants K 2 , . . ., K k and K, the only unknown is the monomer fraction X 1 = ρ 1 /ρ (Eq.22 and 23).As Eq. 22 is a higher order polynomial, it should be (it is) solved numerically.

A. General patchy particle pair potential
To test our extension of TPT1 under strong confinement, we simulate several systems with a variety of potentials, from simple toy systems to more accurate ones.The general expression for the pair interaction V pair between two patchy particles i and j with orientation Ω i and Ω j , respectively, and separated by a distance r ij , is where V rep denotes an isotropic repulsive potential and V p ik ,p jl patch-patch attractive interaction 32 .The position of each patch in the particle reference frame is given by n p unit patch vectors p, which point from the particle's center to the center of the patch (Fig. 4).The min function gives the minimum energy of the set of all possible patch-patch combinations and mimics the fact that we restrict our particles to form only one bond per particle pair.For our systems, the range and width of the patch interaction is relatively small, so that this condition is easily fulfilled.
The attractive patch-patch potential is defined as 32 where V attr (r ij ) is an isotropic attractive potential.As the patches turn away from each other, the patch-patch x j w J 5 U 1 J x U j B i x / Q d 0 q 6 d s Q H t V c 5 E u 1 0 6 x a M 3 J q e J A / K E p I s J q 9 N M H U 9 1 Z s X + l n u q c 6 q 7 T e j v Z L l 8 Y i V G x P 7 l m y n / 6 1 O 1 S A x w q m s Q   1

. Square well pair potential
The simple toy systems employs a hard sphere with diameter σ and a square well attraction (Fig. 5).The square well potential is together with a square well attraction: where δ = 0.005σ, β = 1/k B T the inverse temperature with k B the Boltzmann constant, β ∈ [−20, −5] corresponding to a reduced temperature T * = −1/β ∈ [0.05, 0.20].Note that this square well is rather narrow.The switching function S (Ω) is defined as: • conical, also known as the Kern-Frenkel potential 55,56 : • smooth 57 : • or linear: where θ is the angle between the interparticle vector r and the patch vector p, and θ p is the patch size (Fig. 4).
As shown in the inset in Fig. 5, the patch size θ p was varied from 10, 20, and 30 • for the S KF , S sm. and S lin.switch functions, respectively.

The critical Casimir pair potential
The accurate effective critical Casimir potential of the dipatch particles has its own radial dependence and switching function.In previous work, we optimized the potential (Fig. 6), based on physical dimensions of the dipatch particles and theoretical critical Casimir potentials, to reproduce chain length distributions and bending rigidities observed in experiment over a weak to strong interaction range as function of the temperature 32 .
The repulsive potential is given by the electrostatic Yukawa potential with where Z = πσ 2 c Υ is the charge of the particles, σ c the diameter of the charged colloids, Υ the surface charge density, λ B = βe 2 /4π the Bjerrum length of the solvent with is the permittivity of the solvent, e the elementary charge.The screening length, i.e.Debye length, is defined as κ −1 = k B T /e 2 i ρ i where ρ i is the number density of monovalent ions in the solvent 58 .The model FIG. 6.The patchy particle radial potential for dipatch particles composed of Yukawa repulsion V Yukawa (Eq.33) and critical Casimir attraction VC (Eq.35).The potential uses variables w=0.462,Υ=-0.09e/nm 2 , and θ eff p =19.5 • .The inset shows the switching functions that are additionally a function of dT .
The attractive potential is based on the (isotropic) critical Casimir potential V C between two spherical particles: where A and B are fit parameters depending on the wetting set to w = 0.462 and temperature dT ∈ [0.12, 0.22]K and can be found in Appendix B1 in Ref. 32.The wetting w is determined by the interplay between the patch material and the surrounding binary liquid and the temperature is defined as its distance from the phase separation temperature T cx of the binary liquid dT = T − T cx .The switching function was obtained by explicit integration using V Yukawa and V C at an effective patch width of θ eff p = 19.5 • and fitted to the following functional form: where the coefficients c l are given in Table VI. in Ref. 32.

B. External potential
The external potential V ext at quasi-2D confinement prohibits the translation in the z-direction via: which mimics two hard walls separated by the particle diameter.
The external gravitational potential is composed of two terms: a hard wall represented by a steep Lennard-Jones potential V LJ and a gravitational potential V g that depends, among other things, on the mass of the particle.See the appendix B2 of Ref. 32 for more details.The resulting total gravitational potential is: where LJ is an arbitrary (high) value set to 500k B T and the gravitational force F g is varied from -3.85, to -7.70, and to -11.55 k B T /σ.These gravitational forces ranges from 0.5, 1.0, and 1.5 times the gravitational force of the dipatch particle of interest, respectively.Parameters b and z cut are chosen such that both the potential and the force are continuous at z cut .

C. The system's potential energy
The system's potential energy is a sum over all pair potentials and the external field where i and j run over the N colloidal particles.The external potential V ext is either V quasi−2D or V gravity .

D. The calculation of ∆
The integral in Eq. 8 is performed using Monte Carlo integration via: where the average over the radial distribution function and Mayer function takes place in the integration volume V = dr.The monomers are indicated by α and the other reactant is the γ-particle.Each reaction constant K m in the polymerization is thus defined by the orientational and positional distribution of the reactants and a separate computation of the corresponding ∆ m must be done.Correct determination of the averages is key to calculate ∆.There are two options to measure the average: (1) sample homogeneously over space multiplied with the probability distribution or (2) sample from the correct distribution.This applies to both orientational Ωα,Ωγ as well as the translational V parts of the integration.An advantage of Monte Carlo integration is that it is capable of evaluating both averages simultaneously.
The α-particle, and the γ-particle if it is unbound, are free monomers and their orientational distribution is uniform.Giving them random orientations samples the distribution as option (2).The γ-particle's z-positional distribution is described by the (Boltzmann distribution of the) external potential V ext , and is sampled as option (2).If the γ-particle is restricted to bound configurations, its orientational and positional distributions are additionally described by the pair potential.So instead of uniformly sampling orientational and positional space, its configurations are sampled from a chain with MC (Fig. 2b).Note that this MC sampling of the γ-particle's bonding configurations is independent of the MC sampling of ∆ (Eq.41) and may be performed on-the-fly or beforehand.Thus, the contribution of the orientations Ω α and Ω γ on the average is incorporated using option (2).
The contribution of the inter-particle distance r on the average is sampled by placing the α-particle randomly in a hemispherical shell of volume V around the γ-particle.The radial distribution function then gives the probability of finding the α-particle at distance r.Thus, the contribution of the inter-particle distance on the average is incorporated using option (1).

E. Simulation Details 1. Explicit MC sampling of chain length distributions
Systems with the square well radial potential were simulated with MC for N =1000 divalent particles in a square box with periodic boundary conditions of length 51.17σ or 62.67σ to resemble a density of ρ = 0.382 or 0.255 N/σ 2 , respectively.The densities were chosen to yield reasonable chain length distributions, i.e. that probabilities for longer chain were nonzero.The systems with the critical Casimir potential were taken from Ref. 6.These simulations were done in a rectangular box with periodic boundary conditions of dimensions 43.5σ × 60σ × 43.5σ with 666, or 1000 particles corresponding to ρ = 0.255 or 0.382 N/σ 2 , respectively.
Starting from a random starting configuration MC moves were performed to equilibrate and measure the systems as explained in detail in Ref. 32. Depending on the interaction strength and form of switch function, the equilibration consisted of 1 × 10 4 to 6 × 10 4 MC cycles and the measurements of 5 × 10 4 to 2 × 10 5 MC cycles.Each MC cycle consists of 5 × 10 5 single particle (95%) and cluster moves (5%).
The chain length distribution was measured from three independent simulations after each MC cycle by counting the number N n of chains of length n and normalize by the total number of chains in the system yielding P n = N n / n=1 N n .Irrespective of the use of the discontinuous (square well) or continuous (critical Casimir) attractive potential, a bond is defined for a pair of particles if V p ik ,p jl (r ij , Ω i , Ω j ) < 0k B T .

MC sampling of ∆
For the calculation of the volume average V in Eq. 41 we employ three loops.In the first loop of 10 4 cycles, the configurations of the γ-particle are sampled.For the quasi-2D system, a monomer and dimer configuration are sufficient as there are only two reaction constants.For the gravitational systems, equilibrated chains with length l ≤ 15 are decorrelated with 10 5 single particle MC moves.The position and orientation of the hemisphere of the free site is saved as the γ-particle.In the second loop of 10 2 cycles, the new position of the αparticle is set to r α = r γ + r e where r is a random distance r ∈ [σ, σ + δ] and e a random unit vector pointing to the γ hemi-sphere.The RDF gives the probability of finding the particles at positions r α and r γ in the hard particle reference, more details are in section II E 3. In the third loop of 10 3 cycles, the α-particle is given random orientations and the Mayer function is calculated for the sampling of the V .To improve the calculation of the average, the α-particle was given six sites instead of two.In total there are 2 × 6 possible patch combinations, two from the γ and six from the α-particle, leading to a simple correction of 1/12.The final step to calculate ∆ is to multiply the V with the volume of the hemisphericalshell V .The loops may be repeated one to three times independently, depending on the convergence.

The radial distribution function
The radial distribution function (RDF) of the hard sphere reference fluid is important for the radial component of ∆.The quasi-2D system is radially isotropic in the (x, y)-plane, we thus may use a heuristic RDF g HD (r) of hard disks 59,60 .For the gravitational confined systems, MC simulations are performed to measure the radial distribution of hard spheres at various densities and gravitational fields.The RDF was saved to a file and uses the distance between the particles r = | r 12 |, and the z-coordinates of the particles z 1 and z 2 as variables.A bin width of dr = 0.05σ and dz = 0.01σ in combination with a simple flooring of the bin was sufficient to determine the corresponding RDF during the ∆-calculation.
The reference hard particle diameter d when using the square well potential is simply d = σ as the attractive part of the potential is taken as the perturbation.While for the critical Casimir patchy particle potential, we use the WCA separation to determine d which, in principle, is also a function of the orientation, as V pair is a function of r, Ω i , and Ω j (Eq.26) 61,62 .However, the attractive potential is rather narrow and to avoid a variable reference diameter, d is calculated at maximum attraction, i.e.V pair = V Yukawa (r) + V C (r), and used as a constant.

Numerically solving X1
Given the density and the reaction constants, we determine the monomer fraction X 1 which is bound to X 1 ∈ [0, 1] using the solve-function from the SymPy-Python package 63 .There may be multiple (imaginary) solutions to the equations 18 or 22, and solutions containing a very small imaginary part (Im(X 1 ) < 10 −30 ) were accepted.In case of multiple solutions for X 1 , the smallest one was taken to solve the chain length distribution.

III. RESULTS AND DISCUSSION
To check the predictions of the extended Wertheim theory, we compare predicted and simulated chain length distributions in quasi-2D of systems with a square well radial potential for the three different switching function at a wide range of attractive strengths.The resulting simulated (dots) and predicted (solid lines) chain length distribution are shown in Fig. 7.The top panel shows dipatch particles with a conical switch function S = S KF with patch size θ p =10 • at ρ = 0.382N/σ 2 , the middle panel S = S sm. with θ p =20 • at ρ = 0.255N/σ 2 , and the lower panel a S = S lin. with θ p =30 • at ρ = 0.382N/σ 2 .The hallmark of all distribution is a clear increased monomer density that does not follow the exponential decay of the sequential polymerization, and is well described by the theory.All systems show excellent agreement between model predictions and simulations for all radial interaction strengths, forms of switching function and densities.
These quasi-2D systems allow for determination of the associated excess rotational entropy which can be directly determined from the ratio of chain probabilities.For the initial dimerization, the two monomers lose an excess rotational entropy 2 × F rot and gain a bonding free energy F bond , making exp(−β(F bond − 2F rot )) = P 2 /P 1 = K 2 ρ 1 .For the subsequent polymerization steps, only one monomer loses F rot and F bond is gained, thus exp(−β(F bond − F rot )) = P 3 /P 2 = Kρ 1 .
In Fig. 8, the resulting (exponents) of these free energies are collected for various patch types and sizes, and interaction strength at ρ = 0.382N/σ 2 .For the conical potential with S = S KF , the excess rotational entropy is independent of the interaction strength and only a function of the patch size θ.In that case, the bonding probability is a purely geometrical factor and, as expected, larger patch sizes contain less excess rotational entropy.Non-conical switch functions, namely S sm. and S lin., are dependent on the interaction strength.Here, stronger interaction strengths lead to stiffer chains, and thus to more excess rotational entropy for the monomers.lines) chain length distribution of the square well radial potentials with switching function S = S sm. , and θ p =20 • at ρ = 0.255N/σ 2 with a gravitational force of g=0.5, 1.0 and 1.5 times F g = −7.70kB T /σ.Again, very good agreement between predicted and simulated distributions is observed.Since the predictions included a variational K up to and including chain length l = 8 (Eq.21b-21d), the distributions show a gradual change of the initial slope of the exponential decay, manifesting the gradual change of the reaction constant, as clearly shown in the inset.For g = 1.50, in contrast, the distributions start to resemble more the quasi-2D systems and showing a sharp transition K 2 to K 3 and K 3 ≈ K 4 .This suggest that for higher gravitation one may consider only smaller changes, and fewer reaction constants.Now that we have shown that our extended Wertheim theory is able to incorporate the effects of the gravitational field and the hard wall on the chain length distributions, we can show the role of chain's flexibility on its distribution at finite gravity.The binding rigidity is related to excursions of the bending angle θ due to the thermal fluctuations that are on the order of 1/β .The faster S decays as function of θ (see e.g. the insets of Fig. 5 and 6), the stiffer the bond.When we select two systems with comparable K 2 and K, i.e.F rot and F bond , in quasi-2D, their distributions are similar (Fig. 10a).If these system are in a finite gravitational field instead, their statistics start to diverge significantly from each other (Fig. 10b).We assign this divergence due to the flexibility of the chains: higher stiffness makes the alignment of the particles with the wall more prominent and reduces their reactivity.This example emphasizes the complexity of the orientation-position-dependent re- activity.As a results, there is no direct mapping of the statistics in quasi-2D to finite gravitational strength.Finally, we can apply the extended Wertheim theory to the experimentally relevant system of patchy particles interacting via the critical Casimir force shown in Fig. 11.We can confirm that the approximations done for the RDF, namely using the WCA separation of the repulsion and attraction, and fixing the hard sphere reference diameter d are sufficient.The theory predicts the distribution at various interaction strengths well, except for the strong interaction strengths at dT = 0.14 and 0.12K.At this region, the simulations at dT = 0.14K and ρ = 0.382N/σ 2 show significantly longer chain lengths than predicted.Nematic phases are forming that promote formation of longer chains and slowing down the growth kinetics 6,32 .This is additionally where the simulations have difficulty converging.At dT = 0.12K at both densities, the chains in simulation are still shorter than predicted.Another reason for the difficulty of converging is the relatively narrow switching function in combination with the strong attraction of the critical Casimir potential.The current theory uses an isotropic distribution in the (x, y)−plane and can thus not predict the enhanced reactivity due to the nematic phase.

IV. CONCLUSIONS
In this work we have extended the Wertheim first-order perturbation theory to describe self-assembly of divalent particles under extreme confinement by introducing additional reaction equilibrium constants that account for the reduction of rotational and translational entropy and bond free energy of the polymerization.In the tested systems, the confinement to a monolayer of particles is created by a gravitational field that leads to sub-diameter gravitational heights and an anisotropy of particle density in the direction perpendicular to a wall.
Explicit calculation of these reaction constants from the interaction potential via the integral method allowed for a prediction of the entire chain length distribution functions that agree excellently with direct simulations of these systems.An essential part of the theory is the radial distribution function of the reference hard particle.For finite gravity, this reference hard particle distribution is computed explicitly, but only once for a certain density.The computation of the interaction parameters ∆ can then be done for all densities simultaneously.
In quasi-2D, we can separate the excess rotational entropy from the bonding free energy; the results show that the patch form, size and interaction strength all play a role on the rotational free energy for non-conical potentials, while for conical potentials only the size matters.Additionally, we illustrate that there is no direct or straightforward mapping of the statistics in quasi-2D onto the gravitational systems as the chain's flexibilitythereby availability of the bonds -defines its reactivity.This complex position-orientation dependent reactivity can be explicitly determined by our method.
As one might expect, due to the formation of nematic phases that are, in fact, also observed in experiments, our approach breaks down.One of the assumptions, that the chains will form isotropically in the (x, y)-plane no longer holds.This situation is beyond the scope of the current work.
The advantage of the approach is that it only needs an approximate form of a reference radial distribution function to allow quantitative prediction of an entire range of densities and does not rely on various forms of associating density functionals.Moreover, the approach allows to understand and explain the anomalous small chain concentration of self-assembly under sedimentation conditions.
Finally, we foresee that for other (molecular) selfassembling systems that are described by Wertheim's theory in bulk, our novel extended theory can be applied to describe the system's behavior in extreme confinement or under an external field, e.g. in nanoslits, as the theory is validated not only for toy models but also for realistic potentials.

(
Dated: 23 June 2022) R I 2 m 2 K p P Y l X T 9 6 8 6 m + S / g P 9 F 8 5 u U 1 C L 6 I Y k b 9 + 8 N 7 s z 4 0 a + F 0 v L 6 m e M s f G J y a n s d G 5 m d m 4 + X 1 h Y b M R h I h i v s 9 A P x Y n r x N z 3 A l 6 X n v T 5 S S S 4 0 3 J 9 f u z e 7 K v 4 c Z u L 2 A u D I 9 m J + F n L u Q q 8 S 4 8 5 k q j z Q r 7 Z 5 q w r e u a u W b 5 d 7 6 y e F 0 p W u p 7 5 M 1 7 w a j S M O + P e e B h I j U z q W c K 3 Z T x + A o + I l 8 Y = < / l a t e x i t > r = (x, y) < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 U J o x p 2 2 o u A y / a R 3 o 9 y F 7 j t H I 0 A = " > A A A C 1 H i c j V H L S s N A F D 2 N r 1 o f j b p 0 E 2 w E V y U p o i 5 L B X F Z w T 6 g L S V J p z U 0 L 5 K J U K o r c e s P u N V v E v 9 A / 8 I 7 Y w p q E Z 2 Q m T P n n n N n 7 l w 7 8 t y E G 8 Z r T l l Y X F p e y a 8 W 1 t Y 3 N o v q 1 n Y z C d P Y Y Q 0 n 9 M K 4 b V s J 8 9 y A N b j L P d a O Y m b 5 t s d a 9 t e x i t s h a 1 _ b a s e 6 4 = " w Z p 5 l 2 6 x d 8 T d O 0 a L 4 X W 2 v 5 C M 5 y 0 = " > A A A C z H i c j V H L S g M x F D 0 d X 7 W + q i 7 d D B b B V Z k W U Z d F N 6 6 k g n 1 I W y S T p j V 0 X s x k h F K 6 9 Q f c 6 n e J f 6 B / 4 U 2 c g l p E M 0 x y c u 4 9 J 7 m 5 b u T J R D n O a 8 5 a W F x a X s m v F t b W N z a 3 i t s 7 z S R M Y y 4 a P P T C u O 2 y R H g y E A 0 l l S f a U S y Y 7 3 q i 5 Y 7 O d b x 1 L + J E h s G 1 G k e i 5 7 N h I A e S M 0 X U D e t 2 C 6 5 N 0 2 2 x 5 J m a 1 X u e 5 a Z 4 1 7 e k B l d + t n M e N K v l y n H 5 6 K p a q p 1 l r c 5 j D / s 4 p H 6 e o I Y L 1 N E g b x + P e M K z d W k p a 2 J N P 1 O t X K b Z x b d h P X w A P n 6 R w w = = < / l a t e x i t > a b < l a t e x i t s h a 1 _ b a s e 6 4 = " w Z p 5 l 2 6 x d 8 T d O 0 a L 4 X W 2 v 5 C M 5 y 0 = " > A A A C z H i c j V H L S g M x F D 0 d X 7 W + q i 7 d D B b B V Z k W U Z d F N 6 6 k g n 1 I W y S T p j V 0 X s x k h F K 6 9 Q f c 6 n e J f 6 B / 4 U 2 c g l p E M 0 x y c u 4 9 J 7 m 5 b u T J R D n O a 8 5 a W F x a X s m v F t b W N z a 3 i t s 7 z S R M Y y 4 a P P T C u O 2 y R H g y E A 0 l l S f a U S y Y 7 3 q i 5 Y 7 O d b x 1 L + J E h s G 1 G k e i 5 7 N h I A e S M 0 X U D e t 2 C 6 5 N 0 2 2 x 5 J

FIG. 1 .
FIG. 1. Schematic illustration of possible patchy particle orientations in quasi-2D or under the gravitational field confining the particles close to the wall (striped area).(a) In quasi-2D, i.e. with the translation restricted to the (x, y)−plane, chains cannot rotate around their center-of-mass against the confining wall (red arrow) while monomers can (green arrow).(b) At finite gravitational field, short chains have more freedom to translate against the gravitational force Fg compared to long chains due to the stiffness of the bonds.
r d l P H w A b R m P 7 g = = < / l a t e x i t > z < l a t e x i t s h a 1 _ b a s e 6 4 = " w / 0 o 8 M 0 7 b p c q

e 2 y t 6 6
K v 6 l M y c q 9 r 3 M z v M t b 0 o D t n + O c B o 2 9 s r 1 f r p x V S t V T P e o 8 t r C N X Z r n A a o 4 R g 1 1 8 r 7 C I 5 7 w b J w Y 1 8 a t c f e Z a u S 0 Z h P f l v H w A b g n k T 8 = < / l a t e x i t > y < l a t e x i t s h a 1 _ b a s e 6 4 = " d u u + z V W s R y J L p 2 u Y n w m 3 c 9 O w B 7 8

d 2 S
B m z / H O c s a F U r 9 m H l 4 K x a q p 1 m o 8 5 j B 7 v Y p 3 k e o Y Y 6 G m i S 9 x C P e M K z U T d C I z V u P 1 O N X K b Z x r d l P H w A a F e P 7 A = = < / l a t e x i t > r d l P H w A b R m P 7 g = = < / l a t e x i t > z < l a t e x i t s h a 1 _ b a s e 6 4 = " w / 0 o 8 M 0 7 b p c q

e 2 y t 6 6 FIG. 3 .
FIG.3.A schematic illustration of the calculation of K2 and K4 of particles under a gravitational field.The conformations of the γ-particle are sampled with MC allowing also translations again the gravitational field in the z-direction.

< l a t e x i t s h a 1 _
b a s e 6 4 = " P X b + L o K t W U v d N k C I O C A L y U T z x 4 w = " > A A A C y 3 i c j V H L S s N A F D 2 N r 1 p f V Z d u g k V w V R I p 6 r L o x o 1 Q w T 6 g L W W S T t u h e Z F M h F p d + g N u 9 b / E P 9 C / 8 M 6 Y g l p E J y Q 5 c + 4 5 d + b e 6 0 S e S K R l v e a M h c W l 5 Z X 8 a m F t f W N z q 7 i 9 0 0 j C N H Z 5 3 Q 2 9 M G 4 5 L O G e C H h d C u n x V h R z 5 j s e b z r j c x V v 3 v A 4 E W F w L S c R 7 / p s G I i B c J k k q t W R I y 5 Z L + o V S 1 b Z 0 s u c B 3 Y G S s h W L S y + o I M + Q r h I 4 Y M j g C T s g S G h p w 0 b F i L i u p g S F x M S O s 5 2 d 9 N r U n 0 b W r 3 j I d f 9 N K x a q 9 m 2 l T v K t b 0 o D t n + O c B 4 2 j s n 1 c r l x V S t W z b N R 5 7 G E f h z T P E 1 R x g R r q e o 6 P e M K z c W k k x q 1 x 9 y k 1 c p l n F 9 + W 8 f A B H 6 G S l Q = = < / l a t e x i t > ✓ p < l a t e x i t s h a 1 _ b a s e 6 4 = " m k R d j P X K G B t X + B j l A 9 8 m g r w N + 2 w = " > A A A B 7 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 B I v g q S Q i 6 r H o x W M F W w t t K J v t p F 26 2 c T d i V B K / 4 Q X D 4 p 4 9 e 9 4 8 9 + 4 b X P Q 1 g c D j / d m m J k X p l I Y 8 r x v p 7 C y u r a + U d w s b W 3 v 7 O 6 V 9 w + a J s k 0 x w Z P Z K J b I T M o h c I G C Z L Y S j W y O J T 4 E A 5 v p v 7 D E 2 o j E n V P o x S D m P WV i A R n Z K V W h w Z I r C u 6 5 Y p X 9 W Z w l 4 m f k w r k q H f L X 5 1 e w r M Y F X H J j G n 7 X k r B m G k S X O K k 1 M k M p o w P W R / b l i o W o w n G s 3 s n 7 o l V e m 6 U a F u K 3 J n 6 e 2 L M Y m N G c W g 7 Y 0 Y D s + h N x f + 8 d k b R V T A W K s 0 I F Z 8 v i j L p U u J O n 3 d 7 Q i M n O b K E c S 3 s r S 4 f M M 0 4 2 Y h K N g R / 8 e V l 0 j y r + h d V 7 + 6 8 U r v O 4 y j C E R z D K f h w C T W 4 h T o 0 g I O E Z 3 i F N + f R e X He n Y 9 5 a 8 H J Z w 7 h D 5 z P H y P 8 k A k = < / l a t e x i t > ✓ i < l a t e x i t s h a 1 _ b a s e 6 4 = " x p o i R Y h A 7 i G x Z i e v Z T 6 S o q 4 U m Z c = " > A A A B 7 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 1 G P R i 8 c K 9 g P a U D b b T b t 2 s 4 m 7 E 6 G E / g k v H h T x 6 t / x 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0 T Z x q x h s s l r F u B 9 R w K R R v o E D J 2 4 n m N A o k b w W j m 6 n f e u L a i F j d 4 z j h f k Q H S o S C U b R S u 4 t D j r T 3 0 C t X 3 K o 7 A 1 k m X k 4 q k K P e K 3 9 1 + z F L I 6 6 Q S W p M x 3 M T 9 D O q U T D J J 6 V u a n h C 2 Y g O e M d S R S N u / G x 2 7 4 S c W K V P w l j b U k h m 6 u + J j E b G j K P A d k Y U h 2 b R m 4 r / e Z 0 U w y s / E y p J k S s 2 X x S m k m B M p s + T v t C c o R x b Q p k W 9 l b C h l R T h j a i k g 3 B W 3 x 5 m T T P q t 5 F 1 b 0 7 r 9 S u 8 z i K c A T H c A o e X E I N b q E O D W A g 4 R l e 4 c 1 5 d F 6 c d + d j 3 l p w 8 p l D + A P n 8 w c l g J A K < / l a t e x i t > ✓ j < l a t e x i t s h a 1 _ b a s e 6 4 = " m n X c 1 w g h p 3 R Q 2 Z 0 m 6 6 l 4 h U i J P K c = " > A A A B + H i c b V D L S s N A F L 2 p r 1 o f j b p 0 M 1 g E V y U R U Z d F N y 4 r 2 A e 0 I U y m k 3 b o Z B J m J k I N + R I 3 L h R x 6 6 e 4 8 2 + c t F l o 6 4 G B w z n 3 c s + c I O F M a c f 5 t i p r 6 x u b W 9 X t 2 s 7 u 3 n 7 d P j j s q j i V h H Z I z G P Z D 7 C i n A n a 0 U x z 2 k 8 k x V H A a S + Y 3 h Z + 7 5 F K x W L x o G c J 9 S I 8 F i x k B G s j + X Z 9 G G E 9 C c I s y f 2 M T X P f b j h N Z w 6 0 S t y S N K B E 2 7 e / h q O 2 2 6 q p k S 3 O U v r 5 L u e d O 9 b D r 3 F 4 3 W T V l H F Y 7 h B M 7 A h S t o w R 2 0 o Q M E U n i G V 3 i z n q w X 6 9 3 6 W I x W r H L n C P 7 A + v w B h f 2 T p w = = < / l a t e x i t > p ik < l a t e x i t s h a 1 _ b a s e 6 4 = " P u c K I M I a / 6 a N K 2 d m g l xp T c h T + S I = " >A A A B + H i c b V D L S s N A F L 3 x W e u j U Z d u B o v g q i Q i 6 r L o x m U F + 4 A 2 h M l 0 0 o 6 d T M L M R K g h X + L G h S J u / R R 3 / o 2 T N g t t P T B w O O d e 7 p k T J J w p 7 T j f 1 s r q 2 v r G Z m W r u r 2 z u 1 e z 9 w 8 6 K k 4 l o W 0 S 8 1 j 2 A q w o Z 4 K 2 N d O c 9 h J J c R R w 2 g 0 m N 4 X f f a R S s V j c 6 2 l C v Q i P B A s Z w d p I v l 0 b R F i P g z B L c j 9 7 4 L l v 1 5 2 G M w N a J m 5 J 6 l C i 5 d t f g 2 F M 0 o g K T T h W q u 8 6 i f Y y L D U j n O b V Q a p o g s k E j 2 j f U I E j q r x s F j x H J 0 Y Z o j C W 5 g m N Z u r v j Q x H S k 2 j w E w W M d W i V 4 j / e f 1 U h 1 d e x k S S a i r I / F C Y c q R j V L S A h k x S o v n U E E w k M 1 k R G W O J i T Z d VU 0 J 7 u K X l 0 n n r O F e N J y 7 8 3 r z u q y j A k d w D K f g w i U 0 4 R Z a 0 A Y C K T z D K 7 x Z T 9 a L 9 W 5 9 z E d X r H L n E P 7 A + v w B i Q i T q Q = = < / l a t e x i t > p jl < l a t e x i t s h a 1 _ b a s e 6 4 = " K Y a e 6 L l 1 g x 1 8 a u U n H L j x 1 H e Z 8 c Y FIG.4.A schematic illustration of the interparticle vector r (dotted arrow), patch vectors p on each particle (solid arrows), and the angles θ, and the patch size is defined by the angle θp

FIG. 7 .
FIG. 7.Predicted (solid lines) and simulated (dots) chain length distributions of systems at various interactions strength (β indicated by colors), switch function type S defined in Eq. 30-32, and density ρ [N/σ 2 ].The inset shows the probabilities at small chain lengths.

IG. 8 .
The predicted (exponents of the) excess rotational entropy βFrot and the bonding free energy βF bond of square well radial potentials with various switching functions in quasi-2D confined at ρ = 0.382N/σ 2 .

<
l a t e x i t s h a 1 _ b a s e 6 4 = " w Z p 5 l 2 6 x d 8 T d O 0 a L 4 X W 2 v 5 C M 5 y 0

FIG. 10 .
FIG. 10.The square well radial potentials with switching functions θp = 20 • and S = S lin.(triangles), and θp = 10 • and S = Ssm.(stars) render the same distribution in quasi-2D(a).While when exposed to a gravitational field of g = 0.50, the distribution do not overlap(b).The color coding is the same as in Fig. 8.

FIG. 11 .
FIG. 11.Predicted (solid lines) and simulated (dots) chain length distributions of systems interacting via critical Casimir interactions under a realistic gravitational field at ρ=0.255 and 0.382 N/σ 2 .