A cluster theory for a Janus fluid Riccardo Fantoni1,a) National Institute for Theoretical Physics (NITheP) and Institute of Theoretical Physics, Stellenbosch University, Stellenbosch 7600, 2 South Africa 1 0 2 (Dated: 31 January 2012) n a Recent Monte Carlo simulations on the Kern and Frenkel model of a Janus fluid have J 0 revealed that in the vapour phase there is the formation of preferred clusters made 3 up of a well-defined number of particles: the micelles and the vesicles. A cluster ] h theory is developed to approximate the exact clustering properties stemming from p - the simulations. It is shown that the theory is able to reproduce the micellisation m e phenomenon. h c . s PACS numbers: 64.60.-i, 64.70.-p, 64.70.Fx, 64.60.Ak c i s Keywords: Janus particles, Kern-Frenkel model, Monte Carlo simulation, cluster y h Theory p [ 1 v 6 5 2 6 . 1 0 2 1 : v i X r a a)Electronic mail: [email protected] 1 I. INTRODUCTION In the statistical mechanics of fluids1 the liquid state2 is a particularly fascinating one. A liquid is neither a gas nor a solid, but the state where correlations really play an important role. ThepioneeringworkofBerniJ.Alder3 showedthat,becauseoftheabsenceofattractive forces, thehard-spherefluidadmitsonlyasinglefluidphase. Inordertofindtheliquidphase it is sufficient to add an attractive square-well to the pair-potential of the hard-spheres. The resulting hard-sphere square-well fluid admits a bell-shaped gas-liquid coexistence curve4 with the critical point moving at low temperatures and high densities as the attractive well widthdiminishes. RecentlyN.KernandD.Frenkel5 studied, throughcomputerexperiments, a new fluid model made of hard-spheres with patchy square-well attractions. In its simplest version, the single patch case, the model only depends on the surface coverage χ of the patch and the attraction range. Between the two extreme cases χ = 0, the hard-sphere model, and χ = 1, the hard-sphere square-well model, where the particles pair-potential is isotropic, the particles interaction is directional. The χ = 1/2 model is known as the Janus case, as the particle, like the roman God, has two faces of different functionalities. Another important process, which may lead to the manifestation of macroscopic phe- nomena, in certain fluids, is the clustering or association. In 1956, for example, Leon N. Cooper6 found that the stable state of the degenerate electron fluid in a metal is one in which particles of opposite spin and opposite momentum form pairs. It was then under- stood that whereas the electrons in a metal form pairs with relative angular momentum zero, in 3He this would be prevented by the hard core repulsion, and that therefore Cooper pairing had to occur in a state of finite angular momentum. In 1961 A. Lenard7 proved an- alytically that a two-component plasma living in one dimension undergoes a transition from the conducting to the insulating state by the formation of neutral dimers made of a positive and a negative charge. A two-component plasma living in two dimensions is only stable at sufficiently high temperatures8. But if one adds a hard core to the charges it remains stable even at low temperatures where it undergoes the same transition9. The hard core gives rise to anyonic statistics for the quantum fluid living in two dimensions10. In three dimensions the two-component plasma with a hard core, the so called restricted-primitive model, also undergoes the clustering transition at low temperature and low densities11. An example of a one-component Janus fluid undergoing association is the dipolar hard-sphere fluid. Here 2 a particle can be viewed as the superposition of two uniformly charged spheres: a positive one and a negative one12. In their study of the Kern and Frenkel single patch χ = 1/2 Janus case, F. Sciortino et al.13 found that the gas branch of the coexistence curve bends at high densities at low temperatures. Below the critical point, the fluid tends to remain in the gas phase for a larger interval of densities respect to the χ = 1 case. This behaviour is due to the tendency of particles to associate due to the directional attractive component in the pair-potential and form clusters. At low temperatures, these clusters interact weakly amongst themselves because the particles of which they are composed tend to expose the hard-sphere hemisphere on the outside of the collapsed cluster. By studying the clustering properties of the gas phase of the Janus fluid, F. Sciortino et al. discovered that below the critical temperature there is a range of temperatures where there is formation of two kinds of preferred clusters: the micelles and the vesicles. In the former the particles tend to arrange themselves into a spherical shell and in the latter they tend to arrange themselves as two concentric spherical shells. It is important to confront existing cluster theories with these new findings based oncom- puter experiments. In this work the Bjerrum cluster theory for electrolytes, later extended by A. Tani14 to include trimers, has been employed (preliminary results appeared in Ref. 15) for the description of the exact equilibrium cluster concentrations found in the computer experiment of F. Sciortino et al.. The theory is extended to clusters of up to 12 particles in an attempt to reproduce the micellisation phenomenon observed in the simulations around a reduced temperature of 0.27. A different determination of the intra-cluster configurational partition function has been devised in place of the one used by J. K. Lee16. The Kern and Frenkel fluid has been used to describe soft matter17 biological and non- biological materials like globular proteins in solution5,18,19 and colloidal suspensions5,20, or molecular liquids21. Recently there has been a tremendous development in the techniques for the synthesis of patchy colloidal particles22 in the laboratory. These are particles with dimensions of 10 104˚A in diameter, which obey to Boltzmann statistics23. From the realm − of patchy colloidal particles stems the family of Janus particles for their simplicity24. It is possible to create Janus particles in the laboratory in large quantities25 and to study their clustering properties26. The micelles and the vesicles are complex structures observed in the chemistry of sur- 3 factant molecules analogous to those which may be found in the physical biology of the cell27. The paper is organized as follows: in Section II we describe the fluid model, in Section III we present the clustering properties of the fluid found in the Monte Carlo simulations of F. Sciortino et al., the cluster theory is presented and developed in Sections IV and V, in Section VI we compare the numerical results from our approximation to the exact results of F. Sciortino et al., and Section VII is for final remarks. II. THE KERN AND FRENKEL MODEL As in the work of F. Sciortino et al.13 we use the Kern and Frenkel5 single patch hard- sphere model of the Janus fluid. Two spherical particles attract via a square-well potential only if the line joining the centers of the two spheres intercepts the patch on the surfaces of both particles. The pair-potential is separated as follows Φ(1,2) = φ(r )Ψ(nˆ ,nˆ ,ˆr ) , (1) 12 1 2 12 where + r < σ  ∞ φ(r) =  ǫ σ < r < λσ (2)  − 0 λσ < r     and 1 if nˆ ˆr cosθ and nˆ ˆr cosθ 1 12 0 2 12 0 Ψ(nˆ ,nˆ ,ˆr ) =  · ≥ − · ≥ (3) 1 2 12 0 otherwise  where θ is the angular semi-amplitude of the patch. Here nˆ (ω ) are versors pointing from 0 i i the center of sphere i to the center of the attractive patch, with ω their solid angles and i ˆr (Ω) is the versor pointing from the center of sphere 1 to the center of sphere 2, with Ω its 12 solid angle. We denote with σ the hard core diameter and λ = 1+∆/σ with ∆ the width of the attractive well. A particle configuration is determined by its position and its orientation. We will use σ as the unit of length and ǫ as the unit of energy. 4 One can determine the fraction of the particle surface covered by the attractive patch as follows θ χ = Ψ(nˆ ,nˆ ,ˆr ) 1/2 = sin2 0 , (4) h 1 2 12 iω1,ω2 (cid:18) 2 (cid:19) where ... = ...dω/(4π). ω h i R As in the work of F. Sciortino et al.13 we limit ourselves to the Janus case χ = 1/2. III. CLUSTERING PROPERTIES The Janus fluid just described will undergo clustering as there is a directional attractive component in the interaction between its particles. Moreover at low temperatures the col- lapsed clusters are expected to interact weakly with each other. This is responsible for the bending at high density of the low temperature gas branch of the gas-liquid binodal curve recently determined in Ref. 13. Below the critical temperature, in the vapour phase, the appearance of weakly interacting clusters destabilizes the liquid phase in favour of the gas phase. F. Sciortino et al. during their canonical ensemble (at fixed number of particles N, volume V, and temperature T, with ρ = N/V the density) Monte Carlo simulations of the fluid also studied its clustering properties. In particular they used the following topological definition of a cluster: an ensemble of n particles form a cluster when, starting from one particle, is possible to reach all other particles through a path. The path being allowed to move from one particle to another if there is attraction between the two particles. During the simulation of the fluid they counted the number N of clusters of n particles, which n depends on the particles configurations, and took a statistical average of this number. We show in Fig. 1 the results they obtained for ∆ = σ/2 at a reduced density ρσ3 = 0.01 and various reduced temperatures k T/ǫ. From the figure we can see how at a reduced B temperature of 0.27, in the vapour phase, there is the formation of two kinds of preferred clusters: one made up of around 10 particles and one made up of around 40 particles. 5 1 k T/ε=0.50 B k T/ε=0.40 B 0.1 k T/ε=0.30 B k T/ε=0.27 B 0.01 N 0.001 / > n N 0.0001 < 1e-005 1e-006 1e-007 0 10 20 30 40 50 60 70 n FIG. 1. Exact cluster concentrations of the Janus fluid with ∆ = σ/2 at a reduced density ρσ3 = 0.01 and various reduced temperatures k T/ǫ, from the Monte Carlo simulation of F. B Sciortino et al.13. Intheircollapsedshape, expectedatlowtemperatures, theparticlesintheclusterstendto expose their inactive hemisphere on the outside of the cluster, resulting in a weak interaction between pairs of clusters. In the clusters of around 10 particles the particles tend to arrange themselves into a spherical shell, forming a micellar structure. In the clusters of around 40 particles the particles are arranged into two concentric spherical shells, forming a vesicular structure. The aim of the present work is to see if we can approximate the exact equilibrium cluster concentrations found in the simulation using a cluster theory. We will restrict ourselves to clusters made of up to 12 particles to see if the theory is able to reproduce the micellisation phenomenon. The theory is described next. IV. A CLUSTER THEORY FOR JANUS PARTICLES Following Ref. 14, we describe the fluid of N particles undergoing clustering as a mixture of N species of clusters. Clusters of species n = 1,...,N, which we call n-clusters, are 6 made up of n particles. We denote with N the number of clusters of species n and with n ρ = N /V their density. We assume that the chemical potentials of all the cluster species n n are zero (there is no cost in energy in the formation or destruction of a cluster). Then the grand-canonical partition function of the fluid can be written as N Q = ′ 1 qintra NnQ ( N ,V,T) , (5) tot N ! n inter { n} XNn nY=1 n (cid:0) (cid:1) { } where one separates the coordinates and momenta relative to the center of mass of a cluster from the ones of the center of mass so that qintra will be the intra-cluster partition function n of the cluster of species n and Q the inter-cluster partition function where we consider inter the clusters as non identical. The prime indicates that the sum is restricted by the condition that the number of particles of the fluid is N, N nN = N . (6) n Xn=1 We approximate Q assuming that the sum can be replaced by its largest dominant tot contribution. Using the Stirling approximation N! (N/e)N one then obtains ≈ N lnQ N lnqintra (N lnN N ) +lnQ . (7) tot ≈ n n − n n − n inter Xn=1(cid:2) (cid:3) The maximum of lnQ as a function of N on the constraint of Eq. (6) is given by the tot n { } point N where the gradients of lnQ and of the constraint have the same direction. n tot { } Introducing a Lagrange multiplier λ the equilibrium cluster distribution N is then found n { } from the conditions ∂ lnQ +lnλn = 0 , n = 1,2,3,... (8) tot ∂N (cid:12) n (cid:12) Nn=Nn (cid:12){ } (cid:12) The resulting Helmholtz free energy, βF = lnQ , can then be written in terms of tot tot − the intra-cluster free energy, βfintra = lnqintra, and the inter-cluster partition function as n − n follows N N N βF 1 tot = [ρ lnρ ρ ]+ ρ βfintra + ρ lnV lnQ . (9) V n n − n n n n − V inter Xn=1 Xn=1 Xn=1 where β = 1/k T with k Boltzmann constant and ρ = N /V. B B n n We expect the equilibrium cluster concentrations, N /N, to approximate the ones mea- n sured in the simulation, N /N. n h i 7 V. RELATIONSHIP BETWEEN THE CONFIGURATIONAL PARTITION FUNCTIONS We will assume that Eq. (5) also holds at the level of the configurational partition functions Z, as follows N Z = ′ 1 zintra Nn Z ( N ,V,T) . (10) tot N ! n inter { n} XNn nY=1 n (cid:0) (cid:1) { } In the calculation we only work at the level of the configurational partition functions. Since we expect the clusters to be weakly interacting amongst themselves we will ap- proximate the inter-clusters configurational partition function with: i. the ideal gas approx- imation for pointwise clusters and ii. the Carnahan-Starling approximation28 for clusters of diameter σ . A third possibility, that we have not investigated, would be to use the Boubl´ık, 0 Mansoori, Carnahan, and Starling approximation29,30 for clusters of different diameters σ . n We will only work with a limited number n of different cluster species. Since we are c investigating whether the cluster theory is able to reproduce the micellisation phenomenon we will only consider the first n clusters: n = 1,2,3,...,n . And choosing n = 12. c c c We will describe next the two approximations used for the inter-cluster configurational partition function. A. Ideal gas approximation The simplest possibility is to approximate the mixture of clusters as an ideal one so that Z = VNt , (11) inter where N = N is the total number of clusters. t n n P The equations for the equilibrium numbers of clusters are N = λnVzintra , n = 1,2,3,...,n (12) n n c N = nN , (13) n Xn from which we can determine all the concentrations N /N and the Lagrange multiplier n by solving the resulting algebraic equation of order n . The case n = 2 is described in c c Appendix A. 8 B. Carnahan-Starling approximation A better approximation is found if we use as the inter-cluster configurational partition function the Carnahan-Starling expression28 for hard-spheres of diameter σ , 0 η (4 3η ) t t lnZ = N lnV N − , (14) inter t − t (1 η )2 t − where η = (π/6)ρ σ3 is the clusters packing fraction and ρ = N /V their density. t t 0 t t In this case one needs to solve a system of n +1 coupled transcendental equations, c N = λnVzintraG(η ) , i = 1,2,3,...,n (15) n n t c N = nN , (16) n Xn with η = (π/6)ρ σ3, ρ = N /V, N = N , and t t 0 t t t n n P x(8 9x+3x2) G(x) = exp − . (17) (cid:20)− (1 x)3 (cid:21) − In order to search for the correct root of this system of equations it is important to choose the one that is continuously obtained from the physical solution of the ideal gas approximation as σ 0. Giving a volume to the clusters we introduce correlations between 0 → them which will prove to be essential for a qualitative reproduction of the micellisation phenomenon though the cluster theory. The Carnahan-Starling approximation amounts to choosing for the sequence of virial coefficients of the hard-spheres, a general term which is a particular second order polynomial and to determine the polynomial coefficients that approximate the third virial coefficient by its closest integer28. It could be interesting to repeat the calculation using for the inter-cluster partition function the hard-spheres one choosing allbut thefirst virial coefficient equal tozero, tosee ifthat issufficient toreproduce the micellisation phenomena. Note that in order to study the vesicles we would have to solve a system of around 40 coupled equations. We will describe next how do we determine the intra-cluster configurational partition function zintra. n C. The intra-cluster configurational partition function Toestimatetheintra-clusterconfigurationalpartitionfunctionweperformedMonteCarlo simulations of an isolated topological cluster. 9 We determined the reduced excess internal energy per particle of the n-cluster uex = n n Φ(i,j) /(nǫ) (uex = 0 by definition) as a function of the temperature, and then used h i<j i 1 P thermodynamic integration to determine the intra-cluster configurational partition function. We found that the results for uex(T⋆) can be fitted by a Gaussian as follows n uex(T⋆) = a e bnT⋆2 +c , (18) n n − n with T⋆ = k T/ǫ the reduced temperature. B Given the excess free energy of the n-cluster Fex,intra, we can then determine fex,intra = n n βFex,intra/n as follows n β⋆ fex,intra(β⋆) = uex(1/x)dx n Z n 0 e bn/β⋆2 = c β⋆ +a b  − +√π erf b /β⋆2 1  , (19) n n n n p  b /β⋆2 (cid:20) (cid:18)q (cid:19)− (cid:21) n q   with β⋆ = 1/T⋆ and v = πσ3/6 the volume of the n-cluster. Then the intra-cluster config- 0 0 urational partition function is given by zintra = vnexp( nfex,intra) with zintra = v . n 0 − n 1 0 We studied only the first 10 clusters with n = 3,...,12. The dimer being trivial. To this end we started with an initial configuration of two pentagons with particles at their vertexes juxtaposed one above the other. The two pentagons are parallel to the (x,y) plane, have the z axis passing through their centers, and are placed one at z = +σ/2 and the other at z = σ/2. The particles patches all point towards the origin. We formed the clusters with − a lower number of particles by simply deleting particles and the clusters with 11 and 12 particles by adding a particle on the z axis just above the upper pentagon and just below the lower one. We performed the simulations of the isolated cluster at a fictitious reduced density of ρσ3 = 0.05 which ensured a simulation box big enough that the cluster did not percolate through the periodic boundary conditions. We also compared our results for the excess internal energy calculation for the isolated cluster with the results of F. Sciortino et al. for the low density Janus fluid, from which one extracts cluster information by taking all the clusters found with the same number of particles and averaging their properties, as shown in Fig. 2. 10