Freezing of hard spheres confined in narrow cylindrical pores

Freezing of hard spheres confined in narrow cylindrical pores M. C. Gordillo, B. Martínez-Haya, and J. M. Romero-Enrique Citation: The Journal of Chemical Physics 125, 144702 (2006); doi: 10.1063/1.2358135 View online: http://dx.doi.org/10.1063/1.2358135 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/125/14?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Computational study of the melting-freezing transition in the quantum hard-sphere system for intermediate densities. II. Structural features J. Chem. Phys. 126, 164509 (2007); 10.1063/1.2718525 Computational study of the melting-freezing transition in the quantum hard-sphere system for intermediate densities. I. Thermodynamic results J. Chem. Phys. 126, 164508 (2007); 10.1063/1.2718523 Finite-size effects in the microscopic structure of a hard-sphere fluid in a narrow cylindrical pore J. Chem. Phys. 124, 154708 (2006); 10.1063/1.2187487 Voronoi neighbor statistics of hard-disks and hard-spheres J. Chem. Phys. 123, 074502 (2005); 10.1063/1.2000233 Freezing of hard spheres in confinement J. Chem. Phys. 115, 6538 (2001); 10.1063/1.1401825 Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 01 Apr 2016 08:14:08 Freezing of hard spheres confined in narrow cylindrical pores M. C. Gordilloa兲and B. Martínez-Haya Departamento de Sistemas Físicos, Químicos y Naturales, Universidad Pablo de Olavide, Carretera de Utrera km 1, 41013 Sevilla, Spain J. M. Romero-Enrique Departamento de Física Atómica, Molecular y Nuclear, Área de Física Teórica, Universidad de Savilla, Apartado de Correos 1065, 41080 Sevilla, Spain 共Received 28 July 2006; accepted 5 September 2006; published online 10 October 2006兲 Monte Carlo simulations for the equation of state and phase behavior of hard spheres confined inside very narrow hard tubes are presented. For pores whose radii are greater than 1.1 hard sphere diameters, a sudden change in the density and the microscopic structure of the fluid is neatly observed, indicating the onset of freezing. In the high-density structure the particles rearrange in such a way that groups of three particles fit in sections across the pore. © 2006 American Institute of Physics.关DOI: 10.1063/1.2358135兴 Confined systems have received lot of attention in the past years due mainly to their unusual properties compared to those of the bulk. It is well known that the behavior of any substance can be altered considerably when forced inside a cylinder pore or in between two parallel walls.1–6Extensive studies of constrained species also indicate that the final behavior of the adsorbate depends on the interaction between the particles and the confining walls.7–10 However, this is only a part of the problem since simulation results also suggest that in order to describe completely the involved phenomenology, one should take into account not only the energetics of the considered arrangement but also the steric influence of pure confinement in the properties under scrutiny.11 In systems with soft attractive and repulsive interactions, such as Lennard-Jones-type adsorbates inside a carbon nanotube or a silica pore, the energetic and the entropic influences compete to produce a net effect whose relative contributions cannot be easily separated. On the other hand, in the study of the hard sphere system considered here, energetic effects are suppressed and we are only left with the confinement/entropic influences to contend with. It has been established that a purely one-dimensional system with short-ranged interactions cannot have a phase transition at finite temperature,12–15 although some concerns have been pointed out about the generality of this result.16 Nevertheless, there is profuse information about transitionlike transformations in narrow tubes.5–10 Transformations in confined systems resembling the bulk liquid-solid phase transitions are commonly referred to as “freezing.”1In this paper, freezing represents then a behavior qualitatively similar to a solid-liquid transition 共driven by temperature changes or by purely geometrical packing effects兲, but in a context in which the existence of a real transition is not formally possible or cannot be demonstrated. Our goal is then to study the possible freezinglike transformations in narrow tubes and check the minimum cylinder radius for which the lateral interactions make their occurrence possible, obtaining at the same time the microscopic structure and the equation of state 共EOS兲of the hard spheres systems under consideration. To this latter respect, some advances have already been made,17,18 although related only to the fluid part of the phase diagram. The structure of confined fluids has also deserved considerable attention,19 and it is worthwhile to note that in the incompressible limit 共or equivalently at T=0兲spontane- ous chiral structures have been observed for hard spheres systems.20 In this work, we have performed a series of Monte Carlo 共MC兲simulations in the isothermal-isobaric ensemble 共NPT兲. The confining tubes are taken to be hard and perfectly rigid, i.e., the change in the volume of the simulation cell necessary to conform the system to the fixed pressure requirements was made only in the direction parallel to its z axis, along which periodic boundary conditions were considered. Cylinder pores with fixed constraining radii ranging R=0.75 ␴ –1.50 ␴ , with reduced pressures within P*= ␤ P ␴ 3 =0–20, were considered for the present investigation. In this context, ␴ denotes the diameter of the hard spheres absorbed in the tubes. In most cases the number of particles considered was N=800. Additional test runs were performed with N=1600 and 3200 to make sure that the behavior found was not an artifact of a too small simulation cell. The computations with the greater number of particles were concentrated on the region where the freezing was found, and similar results were obtained in all cases. Averages were performed over 106 Monte Carlo steps, run after 5⫻105equilibration ones. In the freezing region, longer simulations up to 108steps were run. Each MC step comprised an attempt to change the volume of the simulation cell and 2Ntrials to displace a hard sphere 共HS兲particle chosen at random. We tried to move every sphere in a random direction a maximum distance yielding a success ratio of 30%–40%. The maximum length change of the tube allowed at each trial of the NPT simulations was fixed to achieve an acceptance ratio within a兲Electronic mail: [email protected] THE JOURNAL OF CHEMICAL PHYSICS 125, 144702 共2006兲 0021-9606/2006/125共14兲/144702/4/$23.00 © 2006 American Institute of Physics125, 144702-1 Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 01 Apr 2016 08:14:08 20%–30%. The equation of state 共pressure versus density兲of the confined HS fluid was computed in a compression run in which the system pressure was stepwise increased from small pressure up to P*=20, followed by an expansion run back to the limit of negligible pressure. The freezing was in this way characterized by a sizeable change in density when following the equation of state in the compression and expansion runs. In addition, the transition was monitored by the sudden structural rearrangement taking place in the fluid as reflected in the distribution function along the pore axis, g共z兲, and in the structure factor in the zdirection, S共k兲, which is defined as S共k兲=N−1 冓 兺 i,j N exp关ik共zi−zj兲兴 冔 .共1兲 Here, k=2 ␲ n/L, with Lbeing the length of the simulation cell and nan integer number. The zi’s and zj’s coordinates represent the position along the pore axis of the considered pair of spheres. A solidlike structure factor would present a series of spikes at regular intervals, whereas a liquidlike one would have smoother features. The EOS’s of the confined hard spheres in tubes with sizes illustrative of the different behavior regimes observed in the present work are displayed in Fig. 1. In the figure, ␳ =N共 ␲ R2具L典兲−1 stands for the average number density. The dotted line corresponds to the Hall equation of state for the fluid and solid branches of the bulk hard sphere system,21 which is supposed to be an improvement over the traditional Carnahan and Starling expression, especially in the high density limit. A first inspection of Fig. 1reveals that the equilibrium densities corresponding to a given pressure are smaller in the confined fluids investigated in comparison with the bulk fluid, as it was already noted in Ref. 17. As the radius of the tube is increased the EOS of the confined fluid approaches that of the bulk system. However, the main conclusion one can draw from Fig. 1 and its most interesting feature is that, whereas for R⬍1.1 ␴ the pressure increases smoothly with density, for R=1.1 ␴ or greater radii, a neat discontinuity indicative of the freezing transformation is observed in the equation of state. Such feature bears a close resemblance to what happens in the bulk system,21 where we have two branches, a fluidlike one at low pressure and a solidlike one at high pressure, as can be seen in Fig. 1. Noticeably, the freezing enters at smaller densities but higher pressures in the confined fluid in comparison with the bulk. The existence of the discontinuity in the EOS is confirmed in the simulations for the larger pores, which show that with growing Rthe pressure and density of the freezing transformation converge to that of the bulk liquid-solid transition. As mentioned above, we ruled out finite size effects as being responsible for the observed transition by carrying out calculations with particles with N=800, 1600, and 3200 without relevant changes in the results. In order to characterize the freezing in more detail, we have run additional simulations for the N=800 particle system in the R=1.1 ␴ pore for closely spaced pressures around the freezing region 共in the range of P*=14–16兲and up to 108 Monte Carlo steps. We observe that within P*=14.9–15.0 the system becomes highly compressible, so that it explores a large range of densities in the same run, thereby exhibiting broad histograms for the simulation cell length covering the density gap observed in the EOS at freezing. The MC correlation simulation length for the changes in volume also exhibits a sharp increase. This fact resembles the critical slowing down which the Metropolis algorithm manifests close to a critical point. No clear bistable behavior 共characteristic in simulations corresponding to first-order transitions兲was observed. Finally, the simulation cell length fluctuations decrease again as the pressure is increased to P*=16 and larger. That seems to indicate that the freezing is associated to a sudden crossover between two different microscopic structural arrangements in a narrow pressure window, which for practical purposes would be undistinguishable from a true phase transition. Figure 2shows the structure factor for the fluidlike and solidlike states with P*=15, each on one side of the freezing transition, in the pore with R=1.1 ␴ . The higher degree of internal order of the frozen fluid has a correspondence with the presence of marked spikes in S共k兲characteristic of a solidlike structure, in comparison to the smoother and singlepeaked curve for the fluidlike state. It becomes apparent that FIG. 1. Reduced pressure vs numerical density of a system of hard spheres confined inside hard cylinders of different radii. Here, triangles indicate R = ␴ ; crosses, R=1.09 ␴ ; full circles, R=1.1 ␴ ; and squares, R=1.5 ␴ .The dotted lines indicate the corresponding EOS’s for a three-dimensional hard spheres bulk according to Ref. 21. FIG. 2. Structure factor for the case R=1.1 ␴ and ␤ P ␴ 3=15. Full line, equilibrium low density and dashed line, equilibrium high density. 144702-2 Gordillo, Martínez-Haya, and Romero-Enrique J. Chem. Phys. 125, 144702 共2006兲 Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 01 Apr 2016 08:14:08 the sudden change in density in the EOS at P*=15 is correlated with a structure transformation consistent with freezing. The behavior of the axial distribution function, g共z兲, shown in Fig. 3for the same pair of states corroborates this finding. It is observed that a qualitative change in the shape of g共z兲, reflecting a sudden internal ordering in the fluid, takes place upon freezing. The analysis of the present simulation data indicates that the threshold cylinder radius above which the onset of freezing is possible is R⬇1.10 ␴ , and that transformation is associated with a sudden organization of the system in triangular arrangements of hard spheres at roughly the same values of z in correlation with the discontinuity in the EOS. The existence of a threshold radius for freezing may be interpreted in terms of the relevance of the lateral interactions. Since those interactions are heavily determined by the size of the cylinder pore, it is conceivable that there must exist a threshold radius above which lateral interactions become of enough magnitude as to induce the appearance of the structural transformations. The fundamental question then arises of what the expected value of the threshold radius must be. The smaller cylinder pore radii for which two and three hard spheres can fit in a plane with the same axial coordinate zare, respectively, R2=1.00 ␴ and R3=共1/冑3+1/2兲 ␴ ⬇1.08 ␴ . In these two cases, a particle can “jump” over its neighbors and move along the tube, making possible the structural changes. This is opposed to the quasi-one-dimensional system associated with R⬍ ␴ in which every particle has always the same nearest neighbors. As already mentioned and shown below in more detail, our simulations indicate that the frozen branch of the confined fluid is structured as a sequence of trimers of particles. Such trimeric structure provides a more efficient packing of the HS particles in comparison to a dimeric arrangement, which may be the reason for the absence of freezing in pores with R=1.10 ␴ or slightly greater. Freezing becomes then possible only after the trimeric arrangement fits in the pore for R艌1.08 ␴ . Our observation of a threshold at R=1.10 ␴ , slightly delayed with respect to the expected 1.08 ␴ 共no transformation was found for R=1.09 ␴ 兲, may be attributed to a numerical effect linked to the requirement of some additional room for the particles to effectively change their arrangement inside the tubes in the simulations. The trimeric organization of the frozen fluid in the R =1.10 ␴ pore can already be appreciated in the axial distribution function, g共z兲. First of all, for z=0, g共0兲=0 for the fluid branch of the EOS, whereas g共0兲⫽0 on the frozen branch. This implies that freezing is associated with a structure in which neighboring particles share the same z, as opposed to the case of the fluid branch where the occurrence of such situation is negligible. A direct proof for the trimeric arrangements arises from the cumulative integration of the g共z兲, which indicates that every particle has two neighbors at approximately the same zcoordinate. This can be also inferred from the oscillations in g共z兲, which have a roughly threefold greater period in the solidlike structure in comparison to the fluidlike one. In order to characterize in greater detail the trimeric structure of the frozen system, we have considered the correlation function, d共z兲=g共z兲具cos关3共 ␪ i共0兲− ␪ j共z兲兲兴典,共2兲 where ␪ iis the angle which forms the projection of the position of the particle ion the pore section plane with an arbitrary axis contained on it. Figure 4shows that d共z兲presents peaks with alternating signs at the high-density branch, whereas the low-density branch shows appreciably weaker orientational correlations. It is found for the peaks in the frozen branch that d共z兲/g共z兲is almost 1, indicating that the particles form a virtually perfect equilateral triangle at each trimer. Furthermore, the alternating sign shows that orientations of each trimer have a ␲ /3 angle shift with respect to the two nearest-neighboring ones. In this way, consecutive trimers fit in a way similar to the hexagonal packing of the bulk hard sphere solid. We have also analyzed the long-distance behavior of the two axial correlation functions, g共z兲and d共z兲, for the highdensity phase. In both cases, the oscillation amplitudes are modulated by an approximately exponential decaying function at large particle separations. The correlation length defined by the damping factor is smaller in all the cases than half the simulation cell length, and increases as the system is FIG. 3. Axial distribution functions for the zcoordinate at R=1.1 ␴ and ␤ P ␴ 3=15. Full line, equilibrium low density and dashed line, equilibrium high density. FIG. 4. Plot of the orientational correlation function d共z兲for R=1.1 ␴ .The solid line corresponds to the solidlike structure at P*=15 and the longdashed line to the fluidlike one at the same pressure. For comparison, we show d共z兲for the toy model of trimers at P*=10 共short-dashed line兲. 144702-3 Freezing of hard spheres in nanopores J. Chem. Phys. 125, 144702 共2006兲 Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 01 Apr 2016 08:14:08 compressed. So, the solidlike structure does not possess a true long-range ordering, although its short-distance correlation features are analogous to those of a solid phase. The structural properties of the structure of the frozen fluid may be better understood by its comparison with a toy model. We have considered that the fluid forms trimers, each of them composed by three hard spheres constrained to be at the vertices of an equilateral triangle of side length 冑3共R− ␴ /2兲. The triangle fits perfectly into a pore section transversal to its axis, while it is free to rotate around it and to translate along the pore. Hard-core repulsions between nearest-neighboring trimers are considered. This model is effectively one dimensional in nature, and therefore lacks the lateral interactions of the full model presented above and solved by Monte Carlo. The one-dimensional character of this system makes it suitable for analytical study by transfermatrix formalism.18,22–24 In the isothermal-isobaric ensemble, the problem is reduced to obtaining the eigenvalues ␭and eigenfunctions ␺ of the following integral equation: 冕 − ␲ /3 ␲ /3 d ␾ ⬘K共 ␾ − ␾ ⬘兲 ␺ 共 ␾ ⬘兲=␭ ␺ 共 ␾ 兲,共3兲 where the kernel K共 ␾ − ␾ ⬘兲is defined as K= exp 冉 − ␤ P ␴ 3冑1− 冉 R ␴ 冊 2 共1 − cos共 ␾ − ␾ ⬘兲兲 冊 .共4兲 Analysis of Eq. 共3兲shows that the lowest eigenfunction 共which determines the thermodynamical properties of the system23兲is the constant function and consequently ␭0 =兰− ␲ /3 ␲ /3 d ␾ ⬘K共 ␾ − ␾ ⬘兲. Thus, there is no long-range orientational order as expected for a one-dimensional system. However, by the same transfer-matrix technique we obtained the d共z兲function of this toy model and compared it to the corresponding to the full Monte Carlo calculation in Fig. 4. There, it can be easily observed that the correlation function of the solidlike structure resembles closely that of the trimer model. Quantitative differences in the period of the oscillations in d共z兲can be attributed to the strong constrains in our toy model that forces the particles to be trimers, which intrinsically introduces order in the system. We close by noticing that within the precision of our calculations, the observed freezing transition can be regarded as undistinguishable of a true phase transition. In this, we are very close to the conclusions of Koga and Tanaka in Ref. 19 for a similar system of a Lennard-Jones adsorbate inside a carbon nanotube. We acknowledge the financial support of the “Plan Andaluz de Investigación” 共Group FQM-215兲. 1L. D. Gelb, K. E. Gubbins, R. Radhakrishnan, and M. Sliwinska- Bartkowiak, Rep. Prog. Phys. 62, 1573 共1999兲. 2K. Koga, H. Tanaka, and X. C. Zeng, Nature 共London兲408,564共2000兲. 3K. Koga and H. Tanaka, J. Chem. Phys. 112, 104711 共2005兲. 4M. C. Gordillo, J. Boronat, and J. Casulleras, Phys. Rev. Lett. 85,2348 共2000兲. 5Y. Maniwa, H. Kataura, M. Abe, S. Suzuki, and Y. Achiba, J. Phys. Soc. Jpn. 71, 2863 共2001兲. 6Y. Maniwa, H. Kataura, M. Abe, A. Udaka, S. Suzuki, Y. Achiba, H. Kira, K. Matsuda, H. Kadowaki, and Y. Okabe, Chem. Phys. Lett. 401, 534 共2005兲. 7F. R. Hung, B. Coasne, E. E. Santiso, K. E. Gubbins, F. R. Siperstein, and M. Sliwinska-Bartkowiak, J. Chem. Phys. 122, 144706 共2005兲. 8M. W. Maddox and K. E. Gubbins, J. Chem. Phys. 107, 9659 共1997兲. 9M. Sliwinska-Bartkowiak, G. Dudkiak, R. Sikorski, R. Gras, R. Radhakrishnan, and K. E. Gubbins, J. Chem. Phys. 114,950共2001兲. 10 J. Hoffmann and P. Nielaba, Phys. Rev. E 67, 036115 共2003兲. 11 H. Kanda, M. Miyahara, and K. Higashitani, Langmuir 16, 8529 共2000兲. 12 L. van Hove, Physica 共Amsterdam兲16, 137 共1950兲. 13 D. Ruelle, Statistical Mechanics: Rigorous Results 共Addison-Wesley, Reading, MA, 1989兲. 14 E. H. Lieb and D. C. Mattis, Mathematical Physics on One Dimension 共Academic, New York, 1966兲. 15 L. D. Landau and E. M. Lifshitz, Statistical Physics 共Pergamon, New York, 1980兲,Pt.1. 16 J. A. Cuesta and A. Sánchez, J. Stat. Phys. 115, 869 共2004兲. 17 K. K. Mon and J. K. Percus, J. Chem. Phys. 112, 3457 共2000兲. 18 I. E. Kamenetskyi, K. K. Mon, and J. K. Percus, J. Chem. Phys. 121, 7355 共2004兲. 19 K. Koga and H. Tanaka, J. Chem. Phys. 124, 131103 共2006兲. 20 G. T. Pickett, M. Gross, and H. Okuyama, Phys. Rev. Lett. 85, 3652 共2000兲. 21 D. Frenkel and B. Smit, Understanding Molecular Simulation 共Academic, San Diego, 1996兲. 22 J. L. Lebowitz, J. K. Percus, and J. Talbot, J. Stat. Phys. 49, 1221 共1987兲. 23 D. A. Kofke and A. J. Post, J. Chem. Phys. 98, 4853 共1993兲. 24 B. Martínez-Haya, J. M. Pastor, and J. A. Cuesta, Phys. Rev. E 59, 1957 共1999兲. 144702-4 Gordillo, Martínez-Haya, and Romero-Enrique J. Chem. Phys. 125, 144702 共2006兲 Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 01 Apr 2016 08:14:08