Coupled equations describing diﬀusion and cross-diﬀusion of tracer particles in hard-sphere suspen- sions are derived and solved numerically. In concentrated systems with strong excluded volume and viscous interactions the tracer motion is subdiﬀusive. Cross diﬀusion generates transient perturbations to the host-particle matrix, which aﬀect the motion of the tracer particles leading to nonlinear mean squared displacements. Above a critical host-matrix concentration the tracers experience clustering and uphill diﬀusion, moving in opposition to their own concentration gradient. A linear stability analysis indicates that cross diﬀusion can lead to unstable concentration ﬂuctuations in the suspension. The instability is a potential mechanism for the appearance of dynamic and structural heterogeneity in suspensions near the glass transition. 1 Introduction In complex ﬂuids and colloidal suspensions the diﬀusion of tracer particles often cannot be described by classical theories of Brownian motion [1,2]. In concentrated suspensions the mean squared displacement of the tracer particles is a nonlinear function of time, and the probability distribution function is non-Gaussian, owing to hydrodynamic and spatial heterogeneity in the host matrix [1,3 –13]. In order to explore such eﬀects within the generic context of hard-sphere suspensions, experimental methods such as dynamic light scattering [3] and particle tracking [6] have been used extensively over the past several decades, as well as molecular dynamics simulations [13] and statistical mechanical models [1,2]. At high concentrations near the glass transition long-range viscous interactions between the particles make accurate simulations and experiments challenging, and there is uncertainty about the origin of the spatial and temporal heterogeneity within the host matrix [6,12–16]. In this work a continuum model of anomalous diﬀusion in hard-sphere suspensions is developed, in which the diﬀerential equations governing diﬀusion and cross-diﬀusion of tracer particles in concentrated suspensions are solved numerically. The model simulates the development of spatiotemporal heterogeneity within the host matrix caused by the initial tracer distribution, and the eﬀect of this on the motion of the tracer particles. This yields a macroscopic view of the coupling eﬀects between the host matrix and tracer particles, and the onset of anomalous tracer diﬀusion. In Section 2 the ﬂux equations describing diﬀusion and cross-diﬀusion in suspensions are brieﬂy reviewed, and in Section 3 the governing conservation equations are described and put in dimensionless form. The equations are solved in Section 4 for the case of a Gaussian tracer pulse placed within an initially uniform host matrix, as depicted in ﬁgure 1. The macroscopic shape and time behaviour of the perturbation to the host matrix caused by the tracer particles is simulated, as well as the feedback eﬀect on the tracer motion. Section 4.2 contains a discussion of the results and their relation to light scattering and particle tracking experiments. Finally, in Section 5 a linear stability analysis of the governing equations is undertaken, showing that near the glass transition an initially uniform suspension is marginally stable to the growth of normal mode concentration ﬂuctuations. For God so loved the world, that He gave His only begotten Son, that whosoever believeth in Him should not perish, but have everlasting life. For God sent not His Son into the world to condemn the world; but that the world through Him might be saved. John 3:16-17 1 2 x x Figure 1: Schematic of a Gaussian pulse of hard-sphere tracer particles centred at x= 0within a uniform suspension of larger host particles. 2 Diﬀusion and cross-diﬀusion in concentrated suspensions The system to be studied is illustrated in ﬁgure 1, in which a uniform suspension of spherical host particles of radius R1contains a Gaussian pulse of smaller tracer particles of radius R2within the pore space. The ﬂux equations describing diﬀusion and cross-diﬀusion of the host and tracer particles in the one-dimensional system in ﬁgure 1 can be written as J1=D11@n1 @xD12@n2 @x; (1a) J2=D21@n1 @xD22@n2 @x; (1b) wheren1andn2are the concentrations of the host and tracer particles (number of particles per unit volume of mixture) [17]. In quiescent systems with zero volume velocity the host particle ﬂux isJ1=n1u1, whereu1is the average host particle velocity, and J2=n2u2is the ﬂux of tracer particles. In equations (1) D11andD22are the Fickian diﬀusion coeﬃcients of the host and tracer particles, respectively, while D12andD21are cross-diﬀusion coeﬃcients. Phenomenological expressions for the diﬀusion and cross-diﬀusion coeﬃcients as functions of concentration have been obtained previously [18], and are brieﬂy summarized below. 2.1 Diﬀusion coeﬃcients In the dilute tracer limit n2!0, the diﬀusion and cross-diﬀusion coeﬃcients in concentrated hard-sphere suspensions can be written in the approximate form D11=D0 1K^1; D 12=D0 1K= ; (2a,b) D21=D0 2(^`K^1 ); D 22=D0 2 ; (2c,d) whereD0 i=kBT=(6Ri0)is the Stokes-Einstein diﬀusivity of particle i,kBis Boltzmann’s constant,Tis temperature, and 0is the viscosity of the suspending ﬂuid [18]. In (2 a)Kis the dimensionless permeability (viscous mobility) of the host-particle matrix and ^1=  1=(kBT) = d( 1Z)=d 1, where 1= (@=@n1)Tis the derivative of the host-matrix osmotic pressure  with respect to concentration; 1=n1v1is the host-particle volume fraction, v1=4 3R3 1is the volume of a host particle and Z( 1)is the host-matrix compressiblity factor. In (2 b)is the reﬂection coeﬃcient of the host matrix and is the equilibrium partition coeﬃcient of the tracer particles between the pore space and the bulk ﬂuid. In (2 c) is a preferential interaction coeﬃcient accounting for excluded volume eﬀects and ^`= 2`=(2 1);where 2=n2v2is the tracer-particle volume fraction, =R2=R1is the ratio of particle sizes and `= ( +1)= is a viscous cross-diﬀusion factor between the tracer and host particles. Finally, in (2 d) is the diﬀusive tortuosity factor accounting for viscous hindrance to tracer diﬀusion within the pore space. 3 For spherical particles approximate expressions for the friction and compressibility factors of the host-particle matrix are K= (1 1)6:55andZ=1 +a 1 1 1= p; (3a,b) wherea= 41= pand p= 0:64is the volume fraction at random close packing. More accurate expressions for KandZare available [18]; however, equations (3) have a relatively simple form similar to that used recently by Worster et al.[19] and will suﬃce for qualitative purposes. In the dilute limit 1 1equations (3) reduce to the exact results K= 16:55 1andZ= 1+4 1[20,21]. GivenZ,^1can be obtained as ^1=d( 1Z) d 1=1 + 2a 1a 2 1= p (1 1= p)2: (4) Approximate expressions for the partition coeﬃcient and preferential interaction coeﬃcient in hard-sphere suspensions are =eb12 1and =n2b12v1= 2b12=3; (5a,b) where=R2=R1is the size ratio of the tracer and host particles, b12= (1+)3is a thermodynamic coupling coeﬃcient, 2=n2v2is the tracer volume fraction and v2=4 3R3 2is the tracer particle volume. In the limit 1!0and!0;(5a) shows that the partition coeﬃcient is equal to the void fraction, = 1 1:A more general expression for is obtained in the Appendix, where it is shown that equation (5 a) gives qualitatively correct results for all 1when.1: Power law constitutive equations for the diﬀusive tortuosity factor and reﬂection coeﬃcient  in hard sphere suspensions are = (1d)a and= 1(1d)a; (6a,b) where d=3 1 cd(1 1)andcd= 3 g (1 g)(1 +2) (1 +): (7a,b) Heredis the ratio of the tracer size R2to the locally averaged pore size within the suspension, Rd=cdRp, whereRp= 3 1=[R1(1 1)]is the physical pore size [22]. The quantity cdis a scaling factor accounting for ﬂuctuations in the pore size and g= 0:58is the volume fraction at the hard-sphere glass transition. The scaling factor ensures that !0=in the large tracer limit R2!1;where( 1)is the dynamic viscosity of the host suspension. When d1the tracer particles are trapped in the pores and = 0while= 1:The exponents a andain (6) are a = 1cd 3anda=1cd 3; (8a,b) where 1= 2:5=(:22 +)and1=b12(1 + 3+2) +2=(1 +3). Finally, given and, the cross-diﬀusion factor can be obtained as `=(+ 1) and ^`= 2` 2 1: (9) With equations (3)–(9), the diﬀusion coeﬃcients in (2) can be determined as functions of the volume fractions 1and 2and the particle size ratio . 3 Governing equations In ﬁgure 1 a Gaussian tracer pulse of width dand concentration n2=n0 2ex2=2d2is placed at t= 0 within a uniform matrix of the host particles at initial concentration n0 1. The ﬂux equations (1) can be combined with conservation of particle number, @ni @t+@niui @x= 0 (i= 1;2); (10) 4 to yield two coupled diﬀusion equations of the form @n1 @t=@ @x D11@n1 @x+D12@n2 @x ; (11a) @n2 @t=@ @x D21@n1 @x+D22@n2 @x : (11b) Upon replacing the number densities with volume fractions using the relation i=nivi, and introducing the dimensionless length and time scales, ^x=x=dand^t=tD0 2=d2, equations (11) can be written in dimensionless form as @ 1 @^t=@ @^x ^D11@ 1 @^x+3^D12@ 2 @^x ; (12a) @ 2 @^t=@ @^x 3^D21@ 1 @^x+^D22@ 2 @^x ; (12b) where ^Dij=Dij=D0 iand=R2=R1=D0 1=D0 2. 3.1 Initial and boundary conditions At^t= 0the host suspension is of uniform volume fraction 0 1giving the initial conditions 1= 0 1and 2= 0 2e^x2=2(1<^x<1): (13) The boundary conditions are 1! 0 1and 2!0 (^x!1 ): (14) Equations (12)–(14) depend on three dimensionless parameters – the initial volume fractions 0 1; 0 2;and the particle size ratio : 4 Results Equations (12)–(14) were solved using Matlab with a method of lines routine developed by Scheisser and Griﬃths [23] (the Matlab ﬁles are included as Supplemental Material). Figure 2 shows proﬁles of i(^x;^t)for the case = 0:75with initial volume fractions corresponding to a concentrated host matrix ( 0 1= 0:59) and dilute tracer distribution ( 0 2= 0:01). The solid curves in ﬁgure 2aare the computed tracer volume fraction proﬁles at ^t= 0:1(blue) and ^t= 100(red), while the dashed curves are the Gaussian proﬁles obtained by setting the cross-diﬀusion terms to zero (^D12=^D21= 0). In the absence of cross-diﬀusion equation (12 b) has the analytical Gaussian solution 2(^x;^t) = 0 2p 1 + 2 ^D22texp ^x2=2 1 + 2 ^D22t ; (15) which also gives the dashed curves in ﬁgure 2 a. At ^t= 0:1the solid and dashed curves overlap, while at ^t= 100the tracer concentration proﬁle departs from the Gaussian solution. The physical origin of the non-Gaussian diﬀusion at ^t= 100is illustrated in ﬁgure 2 b, which shows the host-matrix volume fraction proﬁles computed from equation (12 a) at^t= 0:1(blue), ^t= 1(orange) and ^t= 100(red). Collective motion of the host particles, induced by the initial tracer pulse through the cross-diﬀusion term in (12 a), leads to a transient perturbation of the host-matrix volume fraction in the vicinity of ^x= 0:This 1(^x;^t)proﬁle in turn aﬀects the tracer particle’s motion via the cross-diﬀusion term in (12 b). 5 Sheet1 Page 1-10-8-6-4-202468100.00E+002.00E-034.00E-036.00E-038.00E-031.00E-021.20E-02 (a)0:012 0:010 0:008 0:006 0:004 0:002 0 2 108642 0 2 4 6 8 10 ^x Sheet1 Page 1-10-8-6-4-202468105.88E-015.90E-01 (b) 10:591 0:590 0:589 0:588 108642 0 2 4 6 8 10 ^x Figure 2: Plots of (a) tracer volume fraction 2(^x)and (b) host matrix volume fraction 1(^x) at times ^t= 0:01(blue), ^t= 0:1(orange) and ^t= 100(red) for the case 0 1= 0:59, 0 2= 0:01 and= 0:75. The dashed curves in (a) are the Gaussian 2proﬁles obtained in the absence of cross-diﬀusion. 4.1 Tracer mean squared displacement Figure 3 shows the tracer mean squared displacement (MSD), calculated as h^x(^t)2i=ZL L^x2 P(^x;^t)P(^x;0) d^x; (16) where 2Lis the domain size and P(^x;^t) = 2(^x;^t)=( 0 2p 2)is the tracer probability density function (thep 2factor ensures thatR1 1P(^x;^t)d^x= 1). The dashed lines show the linear MSD, h^x2i= 2^D22^t;obtained in the absence of cross diﬀusion. At early times ^t.0:1, the perturbation to the host particles has not yet become large enough to aﬀect the tracer particles; the tracer diﬀusion is Gaussian in this regime and the MSDs are linear. As time increases the host matrix becomes more heterogeneous and the tracer MSDs become sublinear, implying transient subdiﬀusive tracer motion. Figure 3 ashows the eﬀect on the MSD of varying the initial host concentration. At volume fractions below 0 10:55the MSDs are approximately linear at all times, similar to observations of anomalous yet Brownian motion [5,6]. At higher host volume fractions above 0 1= 0:59, however, the MSDs become increasingly sublinear, similar to experimental observations of tracer subdiﬀusion near the glass transition [3,4,8,10]. 6 Sheet1 Page 10.010.11100.00010.0010.010.1 101 102 103 104h^x2i 0.01 0.1 1.0 10 ^t(a) Sheet1 Page 10.010.11100.00010.0010.010.1101 102 103 104h^x2i 0.01 0.1 1.0 10 ^t(b) Sheet1 Page 10.010.11100.00010.0010.010.1101 102 103 104h^x2i 0.01 0.1 1.0 10 ^t(c)Normal diﬀusion & Subdiﬀusion! Clustering! Figure 3: (a) Mean squared displacement h^x2icomputed from equation (16) for the case 0 2= 0:01; = 0:75and several initial host matrix volume fractions ( 0 1= 0:55(blue); 0:59(orange); 0:60 (red); 0:605(green); 0:607(purple); 0:61(brown)). (b) Mean squared displacements for the case 0 1= 0:55and 0 2= 0:01and several particle size ratios ( = 0:75(blue); 0:85(orange); 0:95(red); 1:0(green); 1:05(purple)). (c) Mean squared displacements for the case 0 1= 0:55; = 1and several tracer concentrations ( 0 2= 0:0005(blue); 0:001(orange); 0:005(red); 0:008(green); 0:01 (purple)). The dashed lines in (a)-(c) are the Gaussian MSDs h^x2i= 2^D22^tobtained in the absence of cross-diﬀusion eﬀects. 4.1.1 Clustering Above 0 10:607the tracer particles near the origin exhibit clustering; they move at times ^t>1 in a direction opposite to their own concentration gradient (uphill diﬀusion) and the MSD becomes negative (brown curve in ﬁgure 3 a). Similarly to cross-diﬀusion in molecular systems [24], the uphill diﬀusion is caused by the interplay between the normal diﬀusion term and the cross-diﬀusion term on the right-hand-side of the tracer conservation equation (12 b). For relatively low host volume fractions 1the Fickian ^D22term dominates and @ 2=@^t<0everywhere in the system. At host concentrations above the glass transition g= 0:58, however, the tracer diﬀusivity ^D22approaches zero and the cross-diﬀusion ^D21term begins to dominate, so that @ 2=@^t >0near the origin. Another way of interpreting the phenomenon is via the colloidal reﬂection coeﬃcient given by equation (6 a), which approaches 1 as the host volume fraction increases above g[18]. Therefore the host particles have a ﬁltration eﬀect on the tracer particles, pushing them out of the high 1regions into the less-concentrated zone at the origin. This eﬀect of the colloidal reﬂection coeﬃcient on the tracer diﬀusion provides a potential physical mechanism for the ﬁltration behaviour observed in Brownian dynamics simulations of concentrated hard-sphere suspensions by Ariza and Puertas [25]. Clustering of the tracer particles has also been observed in the molecular dynamics simulations of Làzaro-Làzaro et al.[26], and experimentally by Sentjabrskaja et al.[10] – these works show normal Gaussian tracer diﬀusion at low host concentrations leading to subdiﬀusion and ﬁnally clustering of the tracer particles at host-matrix concentrations near the glass transition. 4.1.2 Eﬀect of varying the tracer size and concentration Figure 3billustrates the eﬀect of varying the tracer/host size ratio =R2=R1in a suspension with 0 1= 0:55and 0 2= 0:01. For small tracer particles (  < 0:8) the particles can diﬀuse throughout the pore space; the non-Fickian cross-diﬀusion term in (12 b) is small relative to the Fickian diﬀusion term, and the tracer MSDs are linear. As the tracer size increases the diﬀusion becomes more signiﬁcantly aﬀected by the host matrix; the non-Fickian term in (12 b) begins to control the tracer motion while the Fickian term approaches zero. The tracer MSD becomes increasingly sublinear, leading again to uphill diﬀusion and clustering at a critical size ratio c1. A similar transition sequence also occurs by varying the initial tracer concentration 0 2:Figure 3c shows the eﬀect of increasing the tracer concentration in a suspension with 0 1= 0:55and= 1, so that the tracer particles are of the same size as the host particles (intra-diﬀusion). In this case also there is a transition from normal diﬀusion at suﬃciently low tracer concentrations to subdiﬀusion and ﬁnally clustering as the tracer concentration increases above 0 2= 0:01. 7 Together these results are suggestive of a dynamic phase transition initiated by a subcritical instability of the uniform suspension. That is, in the = 1case the Gaussian tracer pulse plays the role of a small concentration perturbation to the initially uniform host matrix. For suﬃciently dilute tracer concentrations or suﬃciently low host-matrix concentrations 0 1 gthe perturbation decays with time and the suspension is stable. However, at a critical host concentration near to or above the glass transition gthe perturbation tends to grow with time. (The total volume fraction T= 1+ 2always decreases with time and is stable.) The possibility that near gthe host matrix is unstable to concentration ﬂuctuations is explored further in Section 5. 4.2 Discussion The results in ﬁgure 2 depend on the phenomenological expressions for the Dijdescribed in Section 2.1, and in order to have conﬁdence in them it is necessary to test the expressions for the diﬀusion coeﬃcients at the high concentrations studied. Previously it has been shown that the expression (2 d) for the tracer Fickian diﬀusivity D22is in good quantitative agreement with experimental data at volume fractions from dilute up to very near the glass transition g= 0:58 and for several values of the size ratio ; also the predicted value of 1at whichD22!0is in agreement with experiment both for the case of stationary host particles and mobile particles [18]. The equations for the cross-coeﬃcients D12andD21have not yet been tested for hard spheres as cross-diﬀusion data are not presently available. Measured cross-diﬀusion coeﬃcients in polymer solutions that interact mainly via excluded-volume eﬀects, however, show a qualitatively similar concentration dependence [27 –30]. One signiﬁcant shortcoming of the expressions for D12andD21 is the assumption that the equilibrium partition coeﬃcient can be written in the form of equation (5a), which is strictly valid only in the dilute limit 1 1[18]. However, it is shown in the Appendix that in the case 1the equation gives qualitatively correct results at concentrations up to the close-packed limit. The anomalous diﬀusion results in ﬁgure 2 are nevertheless a very robust feature of the cross- diﬀusion model. Qualitatively similar results are obtained, for example, using the more tractable expressions D11=D0 1(1 1= p)2; D 12=aD0 1 1; (17a,b) D21=bD0 2 2; D 22=D0 2(1 1= g)2: (17c,d) In these equations D11is in a form used by Davis and Russel [31] while the form of D22was suggested by Rallison [32], and D12andD21witha=b= 3:5were obtained by Batchelor [17] in the case= 1. These equations satisfy the necessary condition for anomalous diﬀusion to occur within the model – the tracer Fickian diﬀusivity D22approaches zero as 1! gwhileD11; D12 andD21remain ﬁnite. 4.2.1 Dynamic light scattering The tracer mean squared displacements in ﬁgure 3 are qualitatively similar to those observed in experimental studies and molecular dynamics simulations of anomalous diﬀusion in hard-sphere suspensions [3,4,6,8,10,13,26,33]. The physical situations, however, are similar but not equivalent. Equations (11) simulate the evolution of a macroscopic concentration distribution within the suspension, while particle tracking experiments and light scattering studies typically explore the motion of individual particles or microscopic concentration ﬂuctuations in a macroscopically uniform system [3,6]. For normal Gaussian diﬀusion the MSD from either approach is the same; however, in anomalous systems this is no longer the case [34]. Related to this issue is the diﬀerent interpretations of gradient diﬀusion experiments and dynamic light scattering results [35–39]. In dynamic light scattering and particle tracking measurements the particles are observed to undergo Gaussian diﬀusion at early times, reﬂecting cage-rattling, and non-Gaussian diﬀusion at later times, via hopping from cage to cage [3,32]. The continuum model developed here does not capture particle-scale eﬀects such as cage rattling, but simulates the evolution of a macroscopic particle concentration distribution. The similarity of the tracer diﬀusion results to light scattering measurements of the MSD suggests there may nevertheless be an underlying connection between 8 cage dynamics and the development of a macroscopic host-matrix perturbation [16]. The eﬀect of concentration perturbations on the stability of a uniform suspension is explored further in the next section. 5 Linear stability analysis In linear stability analysis a uniform suspension with initial concentrations n0 1andn0 2is perturbed by small sinusoidal ﬂuctuations in the form of normal modes, such that n1=n0 1+a1e!tsinkxandn2=n0 2+a2e!tsinkx; (18a,b) wherekand!are the wave number and growth rate of the perturbation and a1;a2are the amplitudes [35,40,41]. If the growth rate !is less than zero the perturbations decay with time and the suspension is stable, while if ! >0the perturbations grow and the system is unstable; when!= 0the suspension is at marginal stability. For suﬃciently small perturbations ( a1;a2 n0 1;n0 2) the diﬀusion coeﬃcients Dijcan be taken as constants [42,43]. Inserting (18) into (11) then gives !a1=k2a1D11k2a2D12and!a2=k2a1D21k2a2D22: (19a,b) Eliminating the amplitudes a1anda2from (19) and solving the resulting quadratic equation for ! gives the solutions !C=k2DL Cand!S=k2DL S; (20) where DL C=1 2h (D11+D22) +p (D11+D12)24(D11D22D12D21)i ; (21a) and DL S=1 2h (D11+D22)p (D11+D22)24(D11D22D12D21)i ; (21b) are the long-time collective and self diﬀusion coeﬃcients measured by dynamic light scattering experiments [35,37,43]. Setting the growth rate !Sto zero in equation (20) shows that marginal stability occurs when DL S= 0:Light scattering experiments showing that DL S!0as the glass transition is approached [3,44] suggest that an initially uniform suspension is therefore only marginally stable to the growth of concentration perturbations. RequiringDL S0for stability in (21 b) gives D11D22D12D210; (22) which is equivalent to the well-known thermodynamic equilibrium stability condition [35,45,46], but here represents a kinetic instability that has been obtained independently of thermodynamic considerations. A kinetic instability could potentially play a role in the formation of the normal modes and long-range concentration ﬂuctuations observed in hard-sphere suspensions near the glass transition [4,10,26,47,48]. 6 Conclusions Phenomenological equations describing diﬀusion of spherical tracer particles in concentrated hard- sphere suspensions have been solved numerically. A pulse of tracer particles added to an initially homogeneous suspension perturbs the host matrix via cross diﬀusion. Depending on the initial host concentration, several distinct transitions of the tracer diﬀusion occur. At low concentrations, the diﬀusion is Fickian and Gaussian. At higher concentrations the diﬀusion becomes nonGaussian, while the tracer mean squared displacement remains approximately linear. Next the mean squared displacement becomes increasingly sublinear, and eventually a kinetic transition occurs in which the tracer particles move in opposition to their own concentration gradient. The origin of the anomalous diﬀusion is cross-diﬀusion eﬀects between the host suspension and tracer particles, which become more signiﬁcant as the tracer Fickian diﬀusivity decreases near the glass transition. 9 At high concentrations the host suspension acts like a dynamic membrane, ﬁltering the tracer particles out of the more crowded regions of the host matrix. Similar transitions also occur by varying the size or concentration of the tracer particles, at a ﬁxed host volume fraction. The theory reproduces several qualitative aspects of anomalous diﬀusion, including non-Gaussian distribution functions, subdiﬀusive mean squared displacements and clustering near the glass transition. The perturbations to both the tracer and host particle concentration proﬁles have been simulated, giving a macroscopic perspective on the physical mechanisms involved. A linear stability analysis suggests that the glass transition is marginally stable to concentration perturbations, yielding a possible mechanism for the onset of normal-mode like dynamic and structural heterogeneity. Acknowledgements This work was made possible by the love and support of Jesus, Lord of heaven and earth and Saviour of all who believe. Believe on the Lord Jesus Christ and thou shalt be saved. Acts 16:31. I would seek unto God, and unto God would I commit my cause: Which doeth great things and unsearchable; marvellous things without number. Job 5:8,9. A Tracer partition coeﬃcient in concentrated suspensions In previous work an expression for the partition coeﬃcient of tracer particles in a hard-sphere suspension was obtained in the form of equation (5 a) [18]. This equation is strictly valid only in the dilute limit, but in practice can be applied at higher concentrations [30,49]. Here, a more general approach is used to obtain an expression for in concentrated suspensions. In order to introduce necessary concepts a brief review of the previous derivation is ﬁrst given below. Figure 4 shows an osmotic equilibrium set-up for a binary hard-sphere suspension. The suspension in compartment A has particle number densities n1andn2;and total osmotic pressure T; thesuspensionisinequilibriumacrossasemi-permeablepartitionwithareservoir(compartment B)containingonlytheﬂuidsolventand tracerparticlesatnumberdensity nrandosmoticpressure . In the McMillan-Mayer theory of suspensions the solvent is assumed incompressible and composed of molecules much smaller than the hard-sphere colloidal particles, so that thermodynamically the system can be treated like a binary hard-sphere gas, with the total osmotic pressure Tof the suspension playing the same role as the total pressure of the gas [17,50,51]. In this case the Gibbs-Duhem equation of the mixture can be written as, at constant temperature, n1d1+n2d2=dT; (23) wherei(n1;n2)is the chemical potential per particle of component i. Similarly, the Gibbs-Duhem equation applied to the suspension of tracer particles in compartment B is nrdr=d; (24) whereris the chemical potential per particle in the reservoir. At equilibrium the tracer chemical potentials are equal and 2=r: A.1 Osmotic pressures A.1.1 Host particles The host-particle osmotic pressure in A is deﬁned as the diﬀerence between the pressures of compartments A and B, so that  = T[18,50]. In the tracer limit 2!0,depends only onn1and can be written as  =n1kBTZ; (25) whereZ( 1)is the hard-sphere compressibility factor given approximately by (3 b). While (3 b) suﬃces for the qualitative results obtained in Section 4, here a more accurate equation of state for 10 A B f;nrg {T;n1;n2} Figure 4: Illustration of osmotic equilibrium in a bidisperse hard-sphere suspension. The membrane separating compartments A and B is permeable to the ﬂuid (blue) and smaller tracer particles but impermeable to the larger host particles. Because of excluded volume eﬀects, the equilibrium tracer concentration n2in compartment A is less than the concentration nrin the reservoir compartment B; =n2=nris the equilibrium partition coeﬃcient. Tis the total osmotic pressure of the mixture in A, while is the osmotic pressure of the tracer suspension in B. Zis used in order to obtain an accurate equation for :At host volume fractions below 10:55; the Carnahan-Starling equation Z=1 + 1+ 2 1 3 1 (1 1)3(26) agrees very closely (to within 1%) with the results of molecular dynamics simulations [52,53]. However, equation (26) becomes inaccurate as 1increases above the glass transition g= 0:58. At volume fractions above 0:55the compressibility factor is more accurately given by the Speedy equation of state Z= 2:67=(1 = p)[54]. The following empirical equation Z=1 +a 1+ 32 3 155 4 1 (1 1= p); (27) wherea= 41= p;agrees to within 5% with the results of molecular dynamics simulation data at volume fractions up to 0.55, and at higher volume fractions is of comparable accuracy to the Speedy equation [18]. More accurate expressions for Zdeveloped by LeFevre [55], Liu [56] or Paricaud [57] for volume fractions above 0.55 could also be used, but equation (27) has a relatively simple form and is of suﬃcient accuracy for the present purposes. A.1.2 Tracer particles The osmotic pressure of the tracer particles in the reservoir compartment B can be written as =nrkBTZ( r); (28) whereZ( r)is the hard-sphere compressibility factor evaluated at r=v2nr:Equation (28) generalizesthecorrespondingequationfor usedpreviously[18]byallowingthetracerconcentration in the reservoir to be non-dilute. At high host concentrations 1! p, excluded volume eﬀects cause the tracer concentration rin the reservoir to be non-dilute, even though 2in compartment A is dilute. For the case = 1in the close-packed limit 1! p;the reservoir concentration must also approach close packing r! pin order to maintain equilibrium. This phenomenon is discussed in more detail below. A.2 Equilibrium partition coeﬃcient A.2.1 Dilute limit In the dilute limit 2 1 1the chemical potential of the tracer particles in compartment A of ﬁgure 4 can be written as 2=0 2+kBTln 2+ex 2 (29) 11 where0 2is a constant and ex 2=kBTb12 1 (30) is the excess chemical potential caused by interparticle interactions [58]. In the limit r 1the tracer chemical potential in compartment B is (for an incompressible ﬂuid) r=0 2+kBTln r: Equating the chemical potentials at local equilibrium, 2=r;gives the partition coeﬃcient as = 2 r=eb12 1; (31) as obtained previously [18]. Using Widom’s insertion method Lekkerkerker [50] derived a more general expression for the case  1valid at semi-dilute concentrations in the form = (1 1)e(Ay+By2+Cy3); (32) wherey= 1=(1 1); A= 3+ 32+3,B= (9=2)2+ 33andC= 33. A.2.2 Concentrated suspensions In concentrated suspensions the tracer chemical potential in compartment A is 2=0 2+kBTln 2+ex 2 (33) where the excess tracer chemical potential is a more complicated function than equation (30). By considering the work required to insert a tracer particle of radius R2into a uniform sea of host particles of radius R1and volume fraction 1, Heyes and Santos [59] derived an expression for ex 2, valid at arbitrary host-particle concentrations 1, in the form ex 2=kBT(c0+c1+c22+c33): (34) The coeﬃcients in (34) are functions of 1such that c0=ln(1 1); c 1=3 1 1 1; c 2=c1+ 3a2; c 3=Z11 3c12 3c2;(35) whereZis the compressibility factor, a2=Z 1 0Z 1 d 1; (36) andZ=Z(1 + 2 1)=(1 1)2:In the dilute limit 1 1, equation (34) gives ex 2= kBT(1 +)3 1;and agrees with (30). Integrating (36) with equation (27) for Zgives1 a2= 1:311 1+ 1:024 2 1+ 11:733 3 11:721ln(1 1= p)(c0+c1): (37) The chemical potential of the particles in compartment B can be obtained by integrating equation (24), using (28) for , giving r=0 2+kBTZ r 01 rd( rZ) d rd r: (38) InsertingZ( r) = (1 +a r+ 32 3 r55 4 r)=(1 r= p)into (38) and integrating gives2 r=0 2+kBTln r+ex r; (39) 0To ﬁrst order in 2,ex 28kBT 2+in 2(T; 1). 1Alternatively, using the Carnahan-Starling equation (26) gives a2=c0 1=(1 1)2: 2Using the Carnahan-Starling equation for Zin (38) gives ex r=kBT r(89 r+ 3 2 r)=(1 r)3; to ﬁrst order in rthis becomes ex r= 8kBT r. 12 where ex r kBT= 2:621 r+ 3:072 2 r+ 46:933 3 r1:721ln(1 r= p) +2:689 r (1 r= p):(40) To ﬁrst order in requation (40) gives ex r= 8kBT r. Equating the chemical potentials in (33) and (39) gives the partition coeﬃcient as  2 r=e(ex 2ex r)=kBT: (41) Equation (41) is a nonlinear equation for that can be solved numerically using Matlab’s fzero function (Matlab ﬁles included as Supplementary Material). Figure 5 shows equations (41) (solid blue curves), (32) (dashed red curves) and (31) (dash-dotted green curves) plotted versus 1for the case 2= 0:001and several values of :When= 0:03, Lekkerkerker’s [50] expression (32) is in good quantitative agreement with the full model at all concentrations up to the close-packed limit 1= p. Equation (31) gives qualitatively correct results for , though somewhat overestimates the actual value. Sheet1 Page 100.10.20.30.40.50.61E-050.00010.0010.010.11 105104103102101100 10 0:1 0:2 0:3 0:4 0:5 0:6= 1=:03 = 3 " min= 2= p Figure 5: The equilibrium partition coeﬃcient plotted verses 1for the case 2=:001and several values of . The blue solid curves are from equation (41); the red dashed curves are from Lekkerkerker’s expression (32); the green dot-dashed curves are from (31). For= 1, equation (31) gives qualitatively correct results but again overestimates the actual value of . Lekkerkerker’s expression (32) is in good agreement with (41) at concentrations up to 10:3;at higher concentrations 1>0:4it gives a value for that is several orders of magnitude too small as 1! p. As 1increases in the = 1case the partition coeﬃcient approaches the minimum possible value for hard spheres, min= 2= p, shown as the dotted black line in ﬁgure 5. In the close-packed limit 1! pthe tracer particles in compartment A are jammed in between the host particles, and in order to maintain equilibrium the tracer particles in compartment B must also approach the close-packing concentration r! p. For hard spheres the maximum possible concentration in the reservoir is p, and the minimum possible value of the partition coeﬃcient is therefore min= 2= p. For1the full solution (41) approaches this minimum as 1increases. An approximate explicit expression for that has the correct limits is = (1 2= p)eb12 1+ 2= p: (42) For the case = 3both equations (31) and (32) break down beyond the dilute limit, and either (41) or (42) should be used in concentrated suspensions when >1. 13