One-dimensional gas of hard needles
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b One-dimensional gas of hard needles
Yacov Kantor ∗ and Mehran Kardar Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA (Dated: November 6, 2018)We study a one dimensional gas of needle-like objects as a testing ground for a formalism thatrelates the thermodynamic properties of “hard” potentials to the probabilities for contacts betweenparticles. Specifically, we use Monte Carlo methods to calculate the pressure and elasticity coefficientof the hard–needle gas as a function of its density. The results are then compared to the samequantities obtained analytically from a transfer matrix approach.
PACS numbers: 02.70.-c 05.70.Ce 62.10.+s 61.20.Gy 61.30.Cz 05.40.-a 02.50.Ey
I. INTRODUCTION
Due to the relative simplicity of the derivation of theirthermodynamic properties, classical one-dimensional(1D) systems are frequently employed as a test bed oftheory and methods for collective behavior in higher di-mensional systems. For instance, the collection of “hardspheres” on a line, sometimes referred to as the Tonksgas [1], has served as an initial step in the study of two-and three-dimensional systems of hard disks/spheres.There is indeed a general method for exact analysis ofa gas of point particles interacting in 1D via potentialsthat depend only on near-neighbor separations [2]. Here,we employ such methods to study a gas of hard needle-shaped objects. Our object is to compare the analyticalresults with those obtained by (Monte Carlo implemen-tation) of a formalism that relates thermodynamic prop-erties (specifically the pressure and elasticity coefficient)of the gas to the probabilities of contact amongst theparticles.“Hard” potentials, which are either zero or ∞ , help toilluminate the geometrical/entropic features of a thermo-dynamic system. Since there is no energy scale arisingfrom such potentials, the temperature T appears onlyas a multiplicative factor in the free energy and variousother thermodynamic quantities, such as pressure andelastic coefficients. Thus the state of the system becomesindependent of T , and only depends upon such featuresas density. The clarity of the geometrical perspective,combined with the simplicity of numerical simulations,has lead to extensive studies of such systems. In fact,simulations with hard potentials date back to the originsof the Metropolis Monte Carlo (MC) method [3], andhave flourished in the decades that have followed (seeRef. [4] and references therein). A typical example ofnon-trivial behavior is the entropically driven first orderphase transition from a liquid to a solid phase [5].Alignments of non-spherically-symmetric moleculeslead to a diversity of phases in liquid crystals [6]. For ∗ Electronic address: [email protected]
FIG. 1: (Color online) Needle-shaped particles of length 2 ℓ and vanishing thickness are free to rotate in two dimensions,with their centers moving along a line. Particle (needle) i ischaracterized by its translational position x i and orientationangle φ i . example, in the nematic phase the molecules have no po-sitional order (like a liquid), while their orientations arealigned to a specific direction. From the early stages re-search into liquid crystal it was realized that the entropic part of the free energy related to non-spherical shapes ofthe molecules, can by itself explain many of the proper-ties of such systems [7]. Not surprisingly, hard potentialswere frequently invoked, and even such simplifications asinfinitely thin disks [8] or rods [9] have provided valu-able insights regarding liquid crystals. The interplay be-tween the rotational and translational degrees of freedomin molecular solids [10] leads to elastic properties that arecoupled to orientational order. How does one computethe elastic response of such systems from first principles?Recently a formalism enabling direct calculation ofelastic properties and stresses of a system of hardnon-spherically-symmetric objects was developed [11]by extending a previously known formalism for hardspheres [12]. Not surprisingly, given that the elastic re-sponse in two and higher dimensions depends on a rankfour tensor, the resulting expressions contain a largenumber of terms. Typical terms correspond to a vari-ety of possible contacts between particles and numerouscomponents of the separations between them (see, e.g.,Eq. 23 in Ref. [11]). Since these expressions are ob-tained after numerous mathematical transformations, itis advisable to subject them to independent tests. In-deed they have been shown to reduce to the known re-sults for isotropic objects, but up to now there had beenno comparison with exact results for non-spherical par-ticles. Here, we consider the statistical mechanics of a1D system of hard needle-like particles rotating in twodimensions with their centers affixed to a 1D line, as de-picted in Fig. 1. The needles are not allowed to intersect,and thus act as “hard” potentials. This model is a par-ticular case of a group models considered by Lebowitz etal. [13] with anisotropic objects in one dimension. Fromthe perspective of complexity, such systems are a slightgeneralization of the Tonks gas, yet they provide non-trivial insights into the interplay of rotational and trans-lational degrees of freedom. The model can be solved ex-actly, and thermodynamic properties, such as elastic co-efficients, can be calculated. We compare the values ob-tained analytically by the transfer matrix method, withthose from MC simulations using the expressions fromRef. [11] adapted to the 1D case.The paper is organized as follows: The model of hardneedles is introduced in Sec. II, and we demonstrate howthe relative orientations of neighbors leads to an effectivehard-potential as a function of their separation. Sec. IIIis devoted to reviewing how elastic properties of a systemcan be characterized, and the expressions for computingelastic coefficients in 1D are presented. The numericaldifficulties associated with evaluation of various quanti-ties by MC simulations are also described. In Sec. IV wepresent the transfer matrix method for solution of themodel. Details of the MC simulation are presented inSec. V, along with comparisons to results obtained bytransfer matrix method. Discussion and additional fea-tures of the model are presented in Sec. VI. II. MODEL
Figure 1 depicts a configuration of our model, consist-ing of needles of length 2 ℓ , with their center positions re-stricted to move on a 1D line. Needle i is characterized byits position x i , and orientation φ i measured with respectto the normal to the line. Since orientations differing by π are indistinguishable, we restrict − π/ ≤ φ i < π/ directors , frequently appear in thedescription of liquid crystals [6].) As ℓ is the only mi-croscopic length scale in the problem, it can be used toconstruct dimensionless parameters. In particular, themean distance between particles a is made dimensionlessby considering a/ℓ , while the density n can be replacedby nℓ .The needles are not allowed to intersect but do notinteract otherwise. Since the particles cannot cross eachother, we number them (left-to-right) along the 1D line,and require that this order is unchanged, i.e. x i −
Shape and size deformations of objects are usually de-scribed by the strain tensor [14]. In 1D this reduces toa scalar quantity η which simply relates the distortedsize of the system L ′ to its original size L via L ′ = L (1+2 η ) [15]. (While this definition is slightly awkwardin 1D, the use of squared distances between points is con-venient because in higher dimensions it clearly separatestrivial changes in geometry caused by rotations, and realdeformations.) In 1D, for small η , the Helmholtz freeenergy can be expanded as F ( η ) L = F (0) L − pη + 12 Cη + . . . , (3)where C is the elastic coefficient of the body. (Note,that the free energy on the left hand side (l.h.s.) of theequation is divided by the undistorted size of the system.)Consequently, p and C can be calculated from the firstand second derivatives of F with respect to η at fixed T .While the elastic properties are more naturally ob-tained from the Helmholtz free energy, we will also usethe Gibbs free energy G = F + pL , in which the pres-sure is the (imposed) variable [16]. The system size L ,or the mean inter-particle distance is then obtained from a = ∂ ( G/N ) /∂p | T , or in terms of dimensionless variables aℓ = ∂ ( βG/N ) ∂f (cid:12)(cid:12)(cid:12)(cid:12) T . (4)Similarly, C = − a/ (cid:16) ∂a∂p (cid:12)(cid:12)(cid:12) T (cid:17) + p , and βℓC = − a ∂a∂f (cid:12)(cid:12)(cid:12) T + f . (5)Squire et al. [17] developed a formalism for a direct cal-culation of elastic parameters from the correlation func-tions of particles. In this approach stress (pressure) and elastic moduli are related to thermal averages of productsof various inter-particle forces and separations. This for-malism was extended to hard potentials in Refs. [11, 12].Since for hard potentials the forces vanish except whenthe particles touch, the results depend on various contactprobabilities. In two and three dimensions the stress andthe elastic constants are tensors, and the expressions in-volve averages over a variety of components. These re-sults simplify in 1D, and in particular, the expression forstress [Eq. (22) in Ref. [11]] can be considerably simpli-fied: If we denote the separation between two adjacentneedles by s i = x i +1 − x i , they are in contact if the ar-gument of ∆ i = δ [ s i − ℓd i,i +1 ( φ i , φ i +1 )] vanishes. Thedimensionless pressure then f becomes f = nℓ N X i h s i ∆ i i ! . (6)The first term in this expression is simply the pressureof the ideal gas, while the second term can be easily rec-ognized as the mean value of the product of the inter-particle separation and force, as appears in the virial the-orem [18]. To evaluate Eq. (6), we need the probabilitythat two particles ( i and i + 1) with specified orientations( φ i and φ i +1 ) touch each other.Similarly, the elastic coefficient C can be expressed as βCℓ = nℓ N X i h s i ∆ i i + 1 N X i h s i ∆ i i ! (7) − N X i
The partition functions of 1D models with short rangeinteractions can be found analytically using a transfermatrix method [16, 19, 20]. It convenient to considerthe isobaric ensemble with fixed external pressure (force) p , such that (the configurational part of) the partitionfunction is given by Z G = Z N Y i =1 dx i N Y i =1 dφ i e − β P Ni =1 V i − ,i − βpx N . (8)Since x N = P i s i , we can change variables and performintegrations over the separations s i between adjacent par-ticles. For the hard potential given by Eq. (2), this leadsto Z G = ( βp ) − N Z N Y i =1 (cid:16) dφ i e − βpℓd i − ,i ( φ i − ,φ i ) (cid:17) = ( βp ) − N Z N Y i =1 ( dφ i D i − ,i ( f ; φ i − , φ i )) , (9)where D i − ,i ( f ; φ i − , φ i ) ≡ e − fd i − ,i ( φ i − ,φ i ) . (10)(According to the definition of d , all D are identical, ex-cept for D , = D N,N +1 ≡ M equal segments, by setting φ k = − π/ πM k , with k = 0 , , · · · , M −
1. This replaces the function D byan M × M matrix, and the integrals in Eq. (9) are re-placed by matrix products. The partition function thenbecomes Z G = ( βp ) − N ( π/M ) N ˜ vD N − v, (11)where v is a column vector with all of its elements equalto 1. Repeated multiplications can be performed numer-ically, first multiplying D by itself, then multiplying theresulting matrix by itself, etc.. After a total of K suchiterations we arrive at D N , with N = 2 K + 1. The expo-nential dependence on K allows us to achieve very largevalues of N , in practice we used K = 20 in our simula-tions. For moderate pressures, the discretization of theangle φ has little influence on the result once M exceeds10, and we report results for M = 512. (It should benoted, that the same results can also be obtained by nu-merically finding the largest eigenvalue of D . However,in our case this alternative provides no numerical advan-tage.) From the numerical value of Z G , we then obtainthe Gibbs free energyFigure 3 depicts the scaled Gibbs free energy calcu-lated by this numerical procedure. For non-interactingneedles the partition function is Z = ( π/βp ) − N , andthe corresponding βG /N = ln( βp/π ) is indicated by thelower line in the figure. (Both curves exclude the trivialcontribution due to kinetic energy.) f= β pl -4-20246 β G / N FIG. 3: (Color online) The upper curve depicts the Gibbsfree energy per particle (made dimensionless by multiplyingby β ) as a function of the dimensionless pressure f = βpℓ .The lower curve shows, for comparison, the same quantity fornon-interacting needles; the curves begin to separate when f is larger than about one. V. SIMULATIONS AND RESULTS
A Monte Carlo procedure was used to evaluate thepressure and elastic coefficient of the system of hard nee-dles. We simulated N = 128 particles with periodicboundary conditions. Correlations between the needlesfor small and moderate densities do not persist past afew neighbors, and our choice of N caused no percepti-ble finite size effects. An elementary MC move consists ofrandomly choosing a needle, randomly deciding whetherto displace or to rotate it, and attempting to performsuch a move. The move is accepted if in the new posi-tion, or with new orientation, the needle does not overlapwith neighboring needles. The particles are sequentiallyordered, and a position change is rejected if it changesthis ordering. The attempted moves are uniformly dis-tributed over an interval, whose width is chosen to be aslarge as possible, while maintaining an acceptance ratelarger than 50%. The varying size of the interval impliesa diffusion constant for each particle that decreases withincreasing density. A single MC time unit consists of2 N attempts to move or rotate particles. The relaxationtime of the system is proportional to L , and inverselyproportional to the diffusion constant and elastic coef-ficient. (The latter increases with increasing density.)Our choice of elementary step ensured that, within theexamined range of densities, the relaxation time was ap-proximately constant and remained of order N . Thiswas verified by directly measuring several autocorrela-tion functions. For every density n the simulation timewas 5 × N . Such long times are required to ensurehigh accuracy of measured contact probabilities, as ex-plained below.The presence of the Dirac δ -function in the definitionof ∆ i necessitates delicate handling. Both Eq. (6) andEq. (7) require measuring separation s i between the adja- .2 .3 .4 .5 .6 .7 .8 .9 2 3 4 5 6 7 8 f= β pl a /l , β l C FIG. 4: (Color online) The mean distance between particles a (in units of ℓ ), as a function of dimensionless force f , for anideal gas (lower continuous line) as compared with the virialexpansion truncated at the second term (upper continuousline), and with the exact transfer matrix result (dashed line).The dashed-dotted line is the transfer matrix result for theelasticity coefficient C . Solid circles represent the relationbetween a and f from the MC simulations. Solid squaresrepresent the MC results for the elastic coefficient. cent needles at the moment of contact. Such events havezero probability, and the formulas really involve probabil-ity densities . The latter can be evaluated by examiningthe probability that the two needles are within ǫ , anddivide the result by ǫ . Of course, the number of suchnear collisions decreases with decreasing ǫ , and the sta-tistical error increases. The situation is even worse forterms of type h s i s j ∆ i ∆ j i , where two simultaneous con-tacts are supposed to appear. One may define two nearcollision events by considering intervals of size ǫ and ǫ . The opposing requirements of having ǫ i → ǫ i (to ensure statistical accuracy) can be partially recon-ciled by considering each argument of a δ -function beingin the range [ mǫ, ( m + 1) ǫ ), with m = 0 , , · · · , M . Weused M = 10, and ǫ = 0 .
002 (0.01) for high (low) par-ticle density simulations. With 11 data points for singlecontact terms, and 11 points for two-contact terms, wecould view the results as a function of one variable m ,or two variables m and m , and extrapolate the resultsto the “exact contact” limit. The accuracy and prac-ticality of such a procedure has been demonstrated inRef. [12]. The total simulation time was determined byrequirement of having sufficient number of terms in each“bin” of the statistical procedure explained above. Thetotal simulation time was dictated by the need to havea very accurate estimate of the fourth term on the r.h.s.of Eq. (7).Since the MC simulation is performed in the ensem-ble of fixed length, the density or mean inter-particledistance a are given, while the dimensionless pressure f and the dimensionless elastic modulus βℓC are calcu-lated. The full circles in Fig. 4 depict the calculated dependence of f (horizontal axis) on a (vertical axis).The error bars on f are negligible, since Eq. (6) includesonly single pair contacts, and the large statistics as wellas small ǫ ensures very high accuracy. This result iscompared with the relation between f and a obtainedfrom the Gibbs free energy via Eq. (4) by taking the nu-merical derivative of G calculated by the transfer matrixmethod. The latter is depicted by the dashed line. Ex-cellent agreement is obtained between the results fromthese two methods.The solid squares in Fig. 4 depict the MC results for thedimensionless elastic coefficient βℓC (vertical axis), as afunction of the dimensionless pressure f (horizontal axis).Since in the MC procedure f is itself a computed quan-tity, there are now also horizontal error bars, which arenegligible as explained in the previous paragraph. Theaccuracy of C , however, is much lower and depends onboth statistical errors, and systematic errors from ex-trapolation to the true contact probability densities. Wechose the values of ǫ i , and the simulation time in such away that both errors were of the same order. We esti-mate that the vertical error bars are approximately thesize of the symbol for the leftmost point, and decrease tohalf the symbol size for the rightmost point. The dashed-dotted line depicts the same relation obtained from G byusing Eq. (5) and the transfer matrix calculations. Theresults from both approaches coincide within the esti-mated errors. VI. DISCUSSION
The good agreement between the results from MC sim-ulations, based on contact probabilities, and those fromthe transfer matrix method, support the validity of theexpressions reproduced in Eqs. (6), (7) for the pressureand elastic moduli of hard potentials. While limited to1D, this is the first direct comparison between formulaederived in Ref. [11], and an exact alternative approach.We conclude by pointing out an interesting feature ofthe hard needle system: Both at small and large den-sities, the pressure and elastic coefficient are related bythe simple expression C = 2 p , while the behavior at in-termediate densities is more complicated. For low pres-sure (density) the system behaves as an ideal gas with a = ( βp ) − , and substituting in Eq. (5) immediatelyyields C = 2 p in this limit. Interestingly, as discussedby Lebowitz et al. [13], the relation between density andpressure also simplifies at very high pressure (density).In this limit the angular integrations are themselves con-strained by pressure, and a Gaussian approximation leadsto additional powers of βp in the Gibbs partition func-tion. This in turn leads to a density n = a − = 2 βp , i.e.the same functional dependence as an ideal gas but witha factor of 2. Inserting this limiting behavior into Eq. (5)again leads to C = 2 p as in the ideal gas limit.The lower solid line if Fig. 4 depicts the the dependenceof a on f for an ideal gas. The curve deviates from idealbehavior for values of f larger than about 1. At higherdensities (and pressures) we can improve upon the idealgas behavior by using a virial expansion. From the formof the interaction we compute a second virial coefficientof B = 8 /π . As indicated by the upper solid line inFig. 4, inclusion of the second virial coefficient providesa good approximation for f up to 3. Clearly there is nono simple relation between C and f in this intermediateregion.The focus of this paper was to use the model of hardneedles to validate the relation between elastic moduliand contact probabilities for the exactly solvable modelof hard needles. However, the model itself has some in-teresting features, which will be explored elsewhere [21].In particular the simplified behavior alluded above in the high density limit is related to an incipient critical point.The nature and universality of this criticality is relatedto the shapes of the hard objects (in this case, needles). Acknowledgments
This work was supported by the Israel Science Foun-dation Grant No. 193/05 (Y.K.) and by the NationalScience Foundation Grant No. DMR-08-03315 (M.K.)Part of this work was carried out at the Kavli Institutefor Theoretical Physics, with support from the NationalScience Foundation Grant No. PHY05-51164. [1] L. Tonks, Phys. Rev. , 955 (1936)[2] H. Takanishi, Proc. Phys.-Math. Soc. Japan , 60(1942). (Reprinted in Mathematical Physics in One Di-mension , ed. by E. H. Lieb and D. C. Mattis (AcademicPress, New York, 1966), p. 25.[3] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. , 1087 (1953).[4] A. P. Gast and W. B. Russel, Phys. Today , 24 (1998).[5] W. G. Hoover and F. H. Ree, J. Chem. Phys. , 3609(1968), and , 4873 (1967); B. J. Alder, W. G. Hoover,and D. A. Young, J. Chem. Phys. , 3688 (1968).[6] P. G. de Gennes and J. Prost, The Physics of LiquidCrystals , 2nd ed. (Oxford University Press, 1995).[7] L. Onsager, Ann. NY Acad. Sci. , 627 (1949).[8] D. Frenkel and R. Eppenga, Phys. Rev. Lett. , 1089(1982); R. Eppenga and D. Frenkel, Mol. Phys. , 1303(1984); M. A. Bates and D. Frenkel, Phys. Rev. E ,4824 (1998).[9] D. Frenkel and J. F. Maguire, Mol. Phys. , 503 (1983),and Phys. Rev. Lett. , 1025 (1981).[10] D. Fox, M. M. Labes and A. Weissberger, editors, Physicsand Chemistry of Organic Solid State , vol. 1 (Inter-science, New York, 1963).[11] M. Murat and Y. Kantor, Phys. Rev. E , 031124, 2006. [12] O. Farago and Y. Kantor, Phys. Rev. E , 2478, 2000[13] J. L. Lebowitz, J. K. Percus and J. Talbot, J. Stat. Phys. , 1221,1987.[14] L. D. Landau and E. M. Lifshits, Theory of Elasticity (Pergamon Press, Oxford, 1986).[15] T. H. K. Barron and M. L. Klein, Proc. Phys. Soc. ,523 (1965).[16] M. Kardar, Statistical Physics of Particles (CambridgeUniversity Press, Cambridge, 2007).[17] D. R. Squire, A. C. Holt and W. G. Hoover, Physica ,388 (1969).[18] R. K. Pathria, Statistical Mechanics , 2nd ed.(Butterworth-Heinemann, Oxford, 1996); M. Toda,R. Kubo, and N. Saito,
Statistical Physics: EquilibriumStatistical Mechanics Part I (Springer, Berlin, 1998);R. Balescu,
Equilibrium and Nonequilibrium StatisticalMechanics (Wiley, New York, 1975).[19] R. J. Baxter,
Exactly Solved Models in Statistical Me-chanics (Academic Press, London, 1982).[20] L. M. Casey and L. K. Ryunnels, J. Chem. Phys.51