text:
# The effect of entropy on macroions adsorption

## Abstract

We study macroion adsorption on planar surfaces, through a simple model. The importance of entropy in the interfacial phenomena is stressed. Our results are in qualitative agreement with available computer simulations and experimental results on charge reversal and self-assembling at interfaces.

preprint: PREPRINT

The physics of macroions adsorption at charged interfaces is subject of recent theoretical and experimental studies decher97 ; rondelez\_1998 ; tanaka\_2001 ; terao\_2002 ; gelbart00 ; tohver\_2001  because its potential technological applications. It is known that in the neighborhood of a charged surface, low-concentrated solutions of monovalent ions produce an exponentially decaying charge distribution which is well described by Poisson-Boltzmann equation. However, adsorption of multivalent ions displays important deviations from the Poisson-Boltzmann picture and there is not yet a fully accepted theory which well accounts for this situation rondelez\_1998 ; gelbart00 . In inhomogeneous charged fluids, charge reversal (CR) and charge inversion (CI) have been reported as mechanisms with which systems minimize free energy. CR is the overcompensation of the surface’s charge by a contiguous layer of oppositely charged ions from the fluid. Because of electrostatic equilibrium, next to the CR layer, a second layer of ions (with same charge sign as that of the surface) is formed, producing a CI attard\_1996 ; kjellander98 . In molecular engineering, CR and CI are basic mechanisms for self-assembling polyelectrolyte layers on a charged substrate decher97 , reverse mobility experiments lozada99 , self-assembled DNA-lipid membrane complexes raedler\_sci2 , anomalous macroions adsorption on a Lagmuir film rondelez\_1998 , and novel colloids stabilization mechanisms tohver\_2001 .

Although CR was reported since the early 80’s megen80 ; mier\_1989  and some advances in its understanding have been done, the mechanism has not been completely understood. Some authors recognized the relevance of electrostatic and ionic size correlations in CR attard\_1996 ; kjellander98 ; mier\_1989 . In particular, the unsymmetrical size effects (between small coions and counterions) has been considered by Greberg and Kjellander kjellander98 . Their analysis of CR and CI is based on mechanistic arguments and is basically correct, although incomplete: the role of the ionic concentration, entropy and charge asymmetry is not discussed, which are not minor effects tohver\_2001 , as we will see. Recently CR by macroions of an oppositely charged surface was considered by Shklovskii et al. shklovskiiL2000 . They assert in proposing that, when macroions overscreen a charged surface, these macroions are highly ordered in a bidimensional Wigner crystal (WC) structure, as a configuration of minimal free energy. However, because in the WC theory ionic correlations are only considered at the level of a two dimensional electrons liquid at zero temperature, this approach leads to incorrect conclusions. For example, it is predicted that CR decreases by increasing surface charge density (see Fig. 3 in ref. shklovskiiL2000 ) whereas molecular simulations show the opposite jimenez\_2001 ; terao\_2001 . Nevertheless, WC theory has a good agreement with non-zero Monte Carlo simulations of a two-dimensional electron liquid totsuji\_1978  and zero temperature molecular dynamics calculations of a two-dimensional layer of adsorbed ions messina00 . The role of entropy in the attraction of like-charged particles has been recognized in the past muthuJCP\_1996 . Here we will analyze the effect of electrostatic and short range correlations on macroions adsorption at interfaces. A decreasing of the system entropy always produces an increasing of macroions adsorption. Thus, it is concluded that self-assembling at interfaces, CR, CI and like-charged attraction are mainly ruled by entropy.

For the present study we consider a uniformly charged planar surface with $`\sigma _0`$ charges per unit area contiguous to a macroions solution. In our model, we include the two relevant forces of the system: the long range coulombic interaction and the hard-core interaction forces. Therefore macroions are considered to be charged hard spheres (with charge $`Q=Ze`$, diameter $`a`$ and number density $`\rho `$; being $`e`$ the electronic charge) plus an additional amount ($`\rho _c`$) of neutralizing counterions (with charge $`q=z_ce`$ and radius $`r_0=a_c/2`$). The mixture is embedded in a continuous media of dielectric constant $`ϵ`$, set equal to that of the charged surface in order to avoid image charges. By the electroneutrality condition $`\sigma _0=_0^{\mathrm{}}\rho _{el}(x)𝑑x`$ with $`\rho _{el}=_{i=1}^2q_i\rho _i(x)`$, and $`\rho _i(x)`$ being the local concentration profile of the $`i`$-th species. Our model is an extension to that considered by Kjellander and Greberg kjellander98 . The WC theory shklovskiiL2000  has some similarity with our model. However, whereas we use direct (macroion-macroion, macroion-counterion and macroion-wall) coulombic interactions, in WC screened interactions are used. More importantly, short range correlations and entropy (of fundamental importance for finite temperature calculations) are included in our theory jimenez\_2001 .

We study the system by means of integral equations using the hybrid closure hypernetted chain/mean spherical approximation (HNC/MSA). This approximation, for simpler models of charged inhomogeneous fluids, has shown to be in agreement with computer simulations where CR and other effects are studied jimenez\_2001 ; lozadaPRL\_1997 . For the present model, we found qualitative agreement with integral equations results kjellander98  and molecular simulations terao\_2002 . The HNC/MSA equations can be derived in an straightforward manner henderson92a  and are given by

$$g_i(x)=\mathrm{exp}\left\{\beta u_{\alpha i}(x)+\underset{j=1}{\overset{2}{}}\rho _jh_j(y)c_{ij}(s)𝑑v_3\right\},$$
(1)
where $`\rho _i`$ is the bulk concentration of the $`i`$-th species, $`h_j(x)g_j(x)1`$ and $`c_{ij}(s)`$ are the total and direct correlation functions, respectively; $`u_{\alpha i}(x)`$ is the ion-surface interaction potential, $`\beta =1/k_\mathrm{B}T`$ with $`k_B`$ the Boltzmann constant and $`T`$ the absolute temperature.

From HNC/MSA, the concentration profiles ($`\rho _i(x)=\rho _ig_i(x)`$) are obtained directly. These are related to the ion-surface potential of mean force, $`w_i(x)`$, by $`g_i(x)=\mathrm{exp}\{\beta w_i(x)\}`$. Eq. 1 can be written as $`w_i(x)=k_BT\mathrm{ln}[g_i(x)]=\psi _i(x)+J_i(x)`$, where $`\psi _i(x)`$ is the mean electrostatic potential and $`J_i(x)`$ are the short-range contributions. Since both, $`\psi (x)`$ and $`J_i(x)`$ are functionals of $`\rho _{el}(x)`$ and are in a non-linear equation, the charge and size correlations are non-independent. In our theory macroions adsorption is related to both: charge and ionic size correlations. In the Poisson-Boltzmann description the short range correlations are not considered in any way (i. e., $`J_i(x)=0`$). Therefore, PB does not describe correctly many of the adsorption phenomena.

In a system at constant volume V, fixed particle number N and T, its entropy is given by $`SS^{\mathrm{id}}+k\left(\frac{\mathrm{ln}Z_N}{T}\right)_{N,V}`$, where $`S^{\mathrm{id}}`$ is the ideal gas entropy and the configurational integral $`Z_N=\mathrm{}\mathrm{exp}\left\{\beta U_N\right\}𝑑𝐫_1,\mathrm{},d𝐫_N`$, being $`U_N`$ the potential energy of the N interacting particles. The exact computation of $`S`$ is not an easy task. However, from the expressions for $`s`$ and $`Z_N`$ is clear that: (i) the accessible volume my\_note  decreases as the volume fraction ($`\eta _T\frac{\pi }{6}_i\rho _ia_i^3`$) increases. Thus, the larger $`\eta _T`$ the lower the entropy. (ii) Electrostatic interactions between like-charged particles also reduce the system number of accessible configurations (hence entropy) by increasing the strength of the coulombic coupling interaction $`\xi =\beta q_i^2/ϵa_i`$. Less entropy always implies higher order, which in an inhomogeneous fluid translates in higher adsorption, as we will show. In this picture it is plausible charge reversal by monovalent ions, macroions adsorption on uncharged (or like-charged) surfaces, as well as the strong dependence of CR on $`\eta `$. These effects can not be explained only with energy arguments.

We will discuss our results in terms of $`g_i(x)`$, the local net charge per unit area $`\sigma (x)=\sigma _0+\underset{0}{\overset{x}{}}\rho _{el}(y)𝑑y`$ and the absolute maximum of this function ($`\sigma ^{}\mathrm{max}\{\sigma (x):x[0,\mathrm{})\}`$); with the charge profile given by $`\rho _{el}(x)\underset{m=1}{\overset{2}{}}ez_m\rho _mg_m(x)`$. The effective electrostatic force on an ion is $`f_i^e(x)\frac{2\pi }{\epsilon }q_i\sigma (x)`$. Therefore, $`\sigma (x)`$ not only quantifies charge adsorption but also the wall-particle effective electrical force. We have solved HNC/MSA for several values of $`Z`$ and $`\rho `$. In all cases the little ions (counterions) are monovalent, i.e., $`|z_c|=1`$. In all our calculations we have used fixed values of $`a=4`$nm, $`a_c=0.4`$nm, $`T=298`$K and $`ϵ=78.5`$. We present results for macroions next to positively charged, negatively charged and uncharged surfaces. We study the influence on macroions adsorption of $`\sigma _0`$, $`Z`$, $`\eta `$.

In Fig. 1 we show the macroions reduced concentration profile (RCP) for negative macroions ($`Z=20`$, i.e., $`Q=20e<0`$) next to a positively charged surface ($`\sigma _0=0.03`$C/m2) for three different macroions concentrations ($`\eta =0.002,0.01,0.05`$). The macroion RCP shows a very strong adsorption to the wall and a considerable amount of counterions surrounding macroions: In all figures, the counterions RCPs are not plotted for clarity. A second layer of macroions is adsorbed mediated by counterions. In the inset we show the corresponding $`\sigma (x)`$ functions: At a distance of one macroion radius a deep minimum is observed, corresponding to a very strong CR. For $`\eta =0.05`$, the maximum located around $`x=35r_0`$ shows a CI. In the inset, the maximum located at $`x=a/2`$ indicates that small positive ions are adsorbed to the wall. Hence, positive charge is accumulated next to the positively charged surface. We denominate this effect as overcharging. This effect is energetically unfavorable, thus, it must be attributed to the short range correlations since it increases with $`\eta `$. This prediction agrees qualitatively with computer simulation of Tanaka and Grosberg tanaka\_2001 . The effective wall electrical field, which is proportional to $`\sigma (x)`$, is positive before the first layer of macroions and then negative, before the second layer. Therefore, the electrical force is first attractive and then repulsive to negative ions. The behavior of the total force on an ion of species $`i`$, however, is implicit in the RCP, i.e., a $`g_i(x)`$ above (below) 1 implies that $`F_i(x)`$ is attractive (repulsive).

Some of the physics of macroions adsorption is not explicit when considering macroions adsorbed on oppositely charged surfaces. In Fig. 2 we show the macroions RCP for negative macroions ($`Z=20`$) next to a negatively charged surface ($`\sigma _0=0.03`$C/m2). Counterions adsorb on the charged surface (not shown) and, for sufficiently high values of $`\eta _T`$, a layer of macroions is adsorbed to the wall, after a relatively large void of macroions. Hence a surface CR is induced by little ions (according to the $`\sigma (x)`$ functions in the inset) and they, in turn, induce macroions adsorption and the formation of the macroions layer. For sufficiently low values of $`\eta =\frac{\pi }{6}\rho a^3`$ there is no macroions adsorption. The qualitative agreement of our Fig. 2 results for $`\eta =0.01`$, with Monte Carlo data reported in ref. terao\_2002 , is encouraging. From Figs. 1 and 2 it is seen that higher $`\eta _T`$ produces higher adsorption. This prediction is in agreement with computer simulations tanaka\_2001 ; terao\_2001 ; jimenez\_2001 .

In Fig. 3 the RCP are shown for negative macroions next to an uncharged surface. For this macroions volume fraction ($`\eta =0.05`$) we observe colloids adsorption even for low charged colloids \[as low as $`Z=1`$ (not shown)\]. Since the wall is uncharged, the strong macroions adsorption and wall overcharging is due to the low system entropy, which in this case is due to both: high $`\eta `$ and high $`\xi `$. Increasing $`Z`$ increases ($`\xi `$), and thus decreases the entropy and increases macroions adsorption and inducing a surface charge. This induced charge is in agreement with the nanoparticles halo reported by Tohver tohver\_2001 , which seems to be a new colloidal stabilization mechanism. Our results show that there is not a minimum surface or macroions charge for the occurrence of CR. Therefore, the assembled nanostructures at interfaces are stabilized by entropy-induced attractive forces, rather than the surface-macroions electrostatic correlations as it is proposed in the the WC theory. These findings are in agreement with experiments on polyelectrolyte adsorption where it is found that CR of a charged substrate “is more a property of the polymer than a property of the surface” decher97 , and with previous theoretical results, where the importance of particle’s excluded volume (entropy) is recognized as an attractive force mechanism between like charged particles muthuJCP\_1996 . Please notice that CR in Fig. 2 through Fig. 3, is of the same order of magnitude whether the substrate charge is positive, negative or zero.

In Fig. 4, $`\sigma ^{}/\sigma _0`$ is plotted as a function of $`Z`$ for $`\eta =0.002,0.01,0.05`$, for both $`\sigma _0=0.03`$C/m2and $`\sigma _0=0.03`$C/m2(in the inset). For $`\sigma _0=0.03`$C/m2, macroions and surface are both negatively charged, and CR is carried out by the small counterions (see discussion of Fig. 2). Nevertheless $`\sigma ^{}/\sigma _0`$ strongly depends on $`\eta `$, being negligible for low $`\eta `$. For $`\sigma _0<0`$, $`\sigma ^{}/\sigma _0`$ monotonically increases as $`Z`$ with a positive concavity ($`\frac{1}{\sigma _0}\frac{d^2\sigma ^{}}{d^2Z}>0`$) whereas for $`\sigma _0>0`$ $`\frac{1}{\sigma _0}\frac{d^2\sigma ^{}}{d^2Z}<0`$. These results imply that oppositely charged macroions are much more efficiently adsorbed into a charged surface. Our results for oppositely-charged macroions resemble those shown in Fig. 3 of ref. shklovskiiL2000  (WC theory) where $`\sigma ^{}/\sigma _0`$ increases in the same way with $`Z`$ and $`\frac{1}{\sigma _0}\frac{d^2\sigma ^{}}{d^2Z}>0`$. However, the role of entropy is not properly considered in the WC theory shklovskiiL2000 : in WC there is a minimum in $`\sigma _0`$ and $`Z`$ must be large to have overcharging which is in opposition to computer simulations jimenez\_2001  and our theory. When entropy is considered we observe the following: i) $`\sigma ^{}/\sigma _0`$ depends more strongly on $`\eta `$ than on $`Z`$ my\_note2 . ii) For a sufficiently high value of $`\eta `$ our theory predicts $`\sigma ^{}>0`$, even for $`Z=1`$, independently of the value of $`\sigma _0`$ (even for $`\sigma _0=0`$, see Fig. 4). iii) Although not shown explicitly, we found that $`\sigma ^{}/\sigma _0`$ (i.e., CR) increases by increasing $`|\sigma _0|`$ (positive or negative, not shown), which is in agreement with recent simulations results jimenez\_2001 ; terao\_2001  and, as we pointed out above, in opposition to the WC theory.

From Figs. 2 to 4 it is clear that macroions adsorption strongly depends on $`\eta _T`$. This effect is far more important for CR than surface or ionic charge. In an experiment of macroion adsorption on a Langmuir film of amphiphilic molecules Rondelez rondelez\_1998  finds that macroions adsorption is increased (hence CR increases) as the lipid layer charge density decreases, by increasing the surface area. This is consistent with our results since an increase of the area per molecule decreases the surface charge density but not the macroions adsorption which is strongly driven by macroions excluded volume. According to this the CR should also increase by increasing the macroions size and/or concentration. In another important experiment Rädler et al. find that the effective charge on DNA-cationic liposome complexes depends on the DNA-to-cationic lipid mixing ratio, i.e., for a certain DNA-to-DNA distance, $`d_{DNA}`$=$`d_{iso}`$ the complex is neutral. However, if the DNA solution concentration is increased $`d_{DNA}`$ becomes lower than $`d_{iso}`$, which implies higher DNA adsorption on the membranes, i.e., CR on the membrane. This effect can not be explained by energy arguments, but according to our findings an increase in macroion concentration produces higher adsorption.

In summary, a model for macroions solution next a charged surface has been considered. We found that entropy plays a fundamental role in macroions adsorption phenomena and, in consequence, for assembly of organized nanostructures at interfaces. Some of our predictions seem to be in agreement with previous theoretical studies of polyelectrolyte solutions muthuJCP\_1996 , computer simulations jimenez\_2001 ; tanaka\_2001 ; terao\_2002  and experimental results decher97 ; rondelez\_1998 ; tohver\_2001 ; raedler\_sci2 .

We gratefully acknowledge the financial support of INDUSTRIAS NEGROMEX.


format: text

text:
# Specific heats of dilute neon inside long single-walled carbon nanotube

## 1 Introduction

Since the discovery of carbon nanotubes , the peculiar electronic and mechanical properties of these structures have attracted much attention . Experiments have also revealed that they can also be used to store hydrogen  and other gases . Many physicists expected that gases in nanotubes or nanotube bundles may display novel 1-dimensional (1D) behavior as a consequence of the remarkable aspect ratio of the length of tubes to their radius. The group led by Cole  and other researchers  have theoretically studied properties of gases in nanotubes or nanotube bundles. One of the most fantastic properties they found is specific heat of dilute gas inside single-walled carbon nanotubes (SWNT’s) as a functions of temperature: With temperature increasing it shows the thermal behavior changing from 1D, 2D, and 3D. However, there are still two question arisen: Do 1D and 2D behaviors always exist for dilute gas inside SWNT only if the temperature is low enough? Is classical statistical mechanics (CSM) sufficient to deal with this problem?

In Ref. , we has given a brief answer to the first question. Here we will exhibit the full calculations and continue to discuss the second question. In recent papers, S̆iber et al treat the specific heat of dilute He atoms adsorbed in interstitial channels and grooves of carbon nanotube bundles . Here we partially follow their aprroach.

This paper is organized as follows: In Sec.2, we give the expressions of the coordinates of carbon atoms in a unit cell of a SWNT. In Sec.7, we analytically calculate the potential of Ne in a carbon nanotube. For example, we calculate the potential of Ne in (20,20) tube. In Sec.4, we list the related formula of thermodynamics in statistical mechanics to calculate the specific heat. In Sec.5, we calculate the specific heats of Ne inside (20,20) at different temperatures without considering $`\theta ,z`$ effects ($`\theta ,z`$ are defined in Sec.2). In Sec.6, we calculate the specific heats of Ne inside (20,20) at low temperature. In Sec.7, we discuss why we do not use CSM to calculate the specific heat. In this section, we present a coarse method to estimate the characteristic temperature, below which the CSM can not be used. In Sec.8, we discuss the reliability of our results and give brief conclusions.

## 2 the coordinates of carbon atoms in a unit cell of a SWNT

A SWNT without two caps can be constructed by wrapping up a single sheet of graphite such that two equivalent sites of hexagonal lattice coincide . To describe the SWNT, some characteristic vectors require introducing. As shown in Fig.1, the chiral vector $`𝐂_h`$, which defines the relative location of two sites, is specified by a pair of integers $`(n_1,n_2)`$ which is called the index of the SWNT and relates $`𝐂_h`$ to two unit vectors $`𝐚_1`$ and $`𝐚_2`$ of graphite ($`𝐂_h=n_1𝐚_1+n_2𝐚_2`$). The chiral angle $`\theta _0`$ defines the angle between $`𝐚_1`$ and $`𝐂_h`$. For $`(n_1,n_2)`$ nanotube, $`\theta _0=\mathrm{arccos}[\frac{2n_1+n_2}{2\sqrt{n_1^2+n_2^2+n_1n_2}}]`$. The translational vector $`𝐓`$ corresponds to the first lattice point of 2D graphitic sheet through which the line normal to the chiral vector $`𝐂_h`$ passes. The unit cell of the SWNT is the rectangle defined by vectors $`𝐂_h`$ and $`𝐓`$, while vectors $`𝐚_1`$ and $`𝐚_2`$ define the area of the unit cell of 2D graphite. The number $`N`$ of hexagons per unit cell of SWNT is obtained as a function of $`n_1`$ and $`n_2`$ as $`N=2(n_1^2+n_2^2+n_1n_2)/d_R`$, where $`d_R`$ is the greatest common divisor of ($`2n_2+n_1`$) and ($`2n_1+n_2`$). There are $`2N`$ carbon atoms in each unit cell of SWNT because every hexagon contains two atoms. To denote the $`2N`$ atoms, we use a symmetry vector $`𝐑`$ to generate coordinates of carbon atoms in the nanotube and is defined as the site vector having the smallest component in the direction of $`𝐂_h`$. From a geometric standpoint, vector $`𝐑`$ consists of a rotation around the nanotube axis by an angle $`\mathrm{\Psi }=2\pi /N`$ combined with a translation $`\tau `$ in the direction of $`𝐓`$; therefore, $`𝐑`$ can be denoted by $`𝐑=(\mathrm{\Psi }|\tau )`$. Using the symmetry vector $`𝐑`$, we can divide the $`2N`$ carbon atoms in the unit cell of SWNT into two classes: one includes $`N`$ atoms whose site vectors satisfy

$$𝐀_l=l𝐑[l𝐑𝐓/𝐓^2]𝐓(l=0,1,2,\mathrm{},N1),$$
(1)
another includes the remainder $`N`$ atoms whose site vectors satisfy

$`𝐁_l=l𝐑+𝐁_0[(l𝐑+𝐁_0)𝐓/𝐓^2]𝐓`$
$`[(l𝐑+𝐁_0)𝐂_h/𝐂_h^2]𝐂_h(l=0,1,\mathrm{},N1),`$ (2)
where $`𝐁_0\left(\mathrm{\Psi }_0|\tau _0\right)=\left(\frac{2\pi a_{cc}\mathrm{cos}(\theta _0\frac{\pi }{6})}{|𝐂_h|}|a_{cc}\mathrm{cos}(\theta _0\frac{\pi }{6})\right)`$ represents one of the nearest neighbor atoms to A0.

If we introduce a cylindrical coordinate system $`(r,\theta ,z)`$ whose $`z`$-axis is the tube axis. Its $`r\theta `$-plane is perpendicular to $`z`$-axis and contains atom $`A_0`$ in the nanotube. $`r`$ the distance from some point to $`z`$-axis, and $`\theta `$ the angle rotating around $`z`$-axis from an axis which is vertical to $`z`$-axis and passes through atom $`A_0`$ in the tube to the point. In this coordinate system, we can express Eqs.(1) and (2) as

$$𝐀_l=(\rho ,l\mathrm{\Psi },l\tau [l\tau /T]T)(l=0,1,2,\mathrm{},N1),$$
(3)
and

$`𝐁_l`$ $`=`$ $`(\rho ,l\mathrm{\Psi }+\mathrm{\Psi }_02\pi [{\displaystyle \frac{l\mathrm{\Psi }+\mathrm{\Psi }_0}{2\pi }}],l\tau +\tau _0[{\displaystyle \frac{l\tau +\tau _0}{T}}]T)`$ (4)
$`(l=0,1,2,\mathrm{},N1),`$
where $`\rho =\frac{|𝐂_h|}{2\pi }`$. In Eqs.(1)-(4), the symbol $`[\mathrm{}]`$ denotes the largest integer smaller than $`\mathrm{}`$, e.g., $`[5.3]=5`$.

## 3 the potential of ne inside carbon nanotubes

To obtain the potential of Ne inside the nanotube, we firstly consider another simple system shown in Fig. 2: Many atoms distribute regularly in a line form an infinite atom chain and an atom Q is out of the chain. The interval between neighbor atoms in the chain is $`T`$, and the site of atom Q relative to atom 0 can be represented by numbers $`c_1`$ and $`c_2`$. We take the Lennard-Jones potential $`U(R_j)=4ϵ[(\sigma /R_j)^{12}(\sigma /R_j)^6]`$ between atom Q and atom $`j`$ in the chain, where $`R_j`$ is the distance between Q and atom $`j`$, and $`ϵ=\sqrt{ϵ_cϵ_{ne}}`$, $`\sigma =(\sigma _c+\sigma _{ne})/2`$ with $`ϵ_{ne}=35.6`$ K, $`\sigma _{ne}=2.75`$ Å, $`ϵ_c=28`$ K and $`\sigma _c=3.4`$ Å. We calculate the potential between atom Q and the chain as

$$U_{QC}=4ϵ[\sigma ^{12}U_6(c_1,c_2)\sigma ^6U_3(c_1,c_2)],$$
(5)
where $`U_k(c_1,c_2)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\frac{1}{[(c_1+nT)^2+c_2^2]^k}(k=1,2,\mathrm{})`$ which can be calculated through the following recursion :

$$\begin{array}{c}U_1(c_1,c_2)=\frac{\pi \mathrm{sinh}(2\pi c_2/T)}{c_2T[\mathrm{cosh}(2\pi c_2/T)\mathrm{cos}(2\pi c_1/T)]},\hfill \\ U_{k+1}(c_1,c_2)=1/(2kc_2)U_k/c_2.\hfill \end{array}$$
(6)
The $`(20,20)`$ tube, for example, with infinite length can be regarded as $`2N=80`$ chains. Thus the potential of any point Q inside the tube can be calculated as

$$U(r,\theta ,z)=\underset{i=1}{\overset{2N}{}}U_{QC},$$
(7)
where $`(r,\theta ,z)`$ is coordinates of Q in the cylindrical coordinate system. As an approximation, we neglect the potential varying with $`z`$ and $`\theta `$ because we find that it is much smaller than the potential varying with $`r`$ through calculations, and fit the potential with $`U(r)=4\epsilon [(\frac{\stackrel{~}{\sigma }}{\rho r})^{10}(\frac{\stackrel{~}{\sigma }}{\rho r})^4]`$, where $`\rho =13.56`$ Å is the radius of the tube, $`\epsilon =390`$ K, and $`\stackrel{~}{\sigma }=2.63`$ Å (see also Fig.3). Moreover, we simplify it as

$$U(r)=\{\begin{array}{c}4\epsilon [(\frac{\stackrel{~}{\sigma }}{\rho r})^{10}(\frac{\stackrel{~}{\sigma }}{\rho r})^4],r<\rho \stackrel{~}{\sigma },\hfill \\ \mathrm{},r>\rho \stackrel{~}{\sigma }.\hfill \end{array}$$
(8)
## 4  The basic formulae of thermodynamics in statistical mechanics

When we deal with the thermal behavior of dilute Ne in nanotube, we will begin with the free energy $`=𝒯\mathrm{ln}𝒵`$, where $`𝒯`$ is temperature and $`𝒵=_ne^{E_n/𝒯}=\mathrm{tr}\mathrm{exp}(/𝒯)`$ is the partition function. We have let the Boltzmann factor be 1 . $`E_n`$ is the eigenvalues of the Schrödinger equation: $`\psi _n=E_n\psi _n`$. From the free energy, we can derive out the specific heat $`c_v=𝒯(^2/𝒯^2)_v=\frac{E_n^2E_n^2}{𝒯^2}`$, where $`E_n^2=_nE_n^2e^{E_n/𝒯}/𝒵`$ and $`E_n=_nE_ne^{E_n/𝒯}/𝒵`$.

We call above statistical mechanics ordinary statistics (OS) and list the corresponding formula in CSM: The free energy $`_{cl}=𝒯\mathrm{ln}𝒵_{cl}`$, where $`𝒵_{cl}=^{}e^{E(p,q)/𝒯}𝑑\mathrm{\Gamma }`$ is the partition function and here the prime means we integrate only over the regions of phase space which correspond to physically different states of particles. The specific heat $`c_{vcl}=𝒯(^2_{cl}/𝒯^2)_v=\frac{E^2E^2}{𝒯^2}`$, where $`E^2=^{}E^2e^{E(p,q)/𝒯}𝑑\mathrm{\Gamma }/𝒵_{cl}`$ and $`E=^{}Ee^{E(p,q)/𝒯}𝑑\mathrm{\Gamma }/𝒵_{cl}`$.

Under the sufficient high temperature, $``$ converge to the $`_{cl}`$. In other words, the applicable domain of $``$ is larger than that of $`_{cl}`$. Therefore, the conclusion derived from $``$ is much more reliable than that derived from $`_{cl}`$. We will use OS in the following discussions if we do not make special statement.

## 5  Specific Heats of Ne gas in an approximate potential

Because we consider the dilute Ne, we can neglect the interaction between Ne atoms and write the single particle Schrödinger equation  as $`H\psi =E\psi `$, where $`H=\frac{\mathrm{}^2}{2\mu }^2+U(r)`$ and $`\psi =\varphi e^{i(m\theta +\kappa z)}`$. It follows that

$$\begin{array}{c}E=\frac{\mathrm{}^2\kappa ^2}{2\mu }+E_m(\kappa ,m=0,\pm 1,\pm 2,\mathrm{}),\hfill \\ H(r)\varphi =E_m\varphi ,\hfill \\ H(r)=\frac{\mathrm{}^2}{2\mu }(\frac{d^2}{dr^2}+\frac{1}{r}\frac{d}{dr}\frac{m^2}{r^2})+U(r),\hfill \end{array}$$
(9)
where $`m`$ is angular quantum number and $`E_m`$ is the corresponding energy. Setting $`r=(\rho \stackrel{~}{\sigma })\xi `$, $`\epsilon _0=\frac{\mathrm{}^2}{2\mu \rho ^2}`$ and $`\eta =\stackrel{~}{\sigma }/\rho `$, Eqs.(8) and (9) are transformed into

$$u(\xi )=\{\begin{array}{c}4\epsilon [(\frac{\eta }{1(1\eta )\xi })^{10}(\frac{\eta }{1(1\eta )\xi })^4],\xi <1,\hfill \\ \mathrm{},\xi >1.\hfill \end{array}$$
(10)
and

$$\begin{array}{c}\phi (\xi )=E_m\phi (\xi ),\hfill \\ =\frac{\epsilon _0}{(1\eta )^2}(\frac{d^2}{d\xi ^2}+\frac{1}{\xi }\frac{d}{d\xi }\frac{m^2}{\xi ^2})+u(\xi ).\hfill \end{array}$$
(11)
If we let $`|\phi =_na_n|\chi _n`$, we will obtain the secular equation

$$det(_{jn}E_m𝒮_{jn})=0,$$
(12)
where $`_{jn}=_0^1\chi _j(\xi )\chi _n\xi 𝑑\xi `$, and $`𝒮_{jn}=_0^1\chi _j\chi _n\xi 𝑑\xi `$. If we let $`\chi _n=J_{|m|}(\nu _n\xi )`$, where $`J_{|m|}(\xi )`$ is the m-th order Bessel function of the first class and $`\nu _n`$ is the n-th zero point of Bessel function , we can calculated $`E_{mn},(m=0,\pm 1,\pm 2,\mathrm{};n=1,2,3,\mathrm{})`$ from Eq.(12).

If there are $`𝒩`$ Ne atoms inside the tube, we have the free energy $`=𝒩𝒯\mathrm{ln}𝒵_1`$, where $`𝒵_1=_{mn}e^{E_{mn}/𝒯}_{\mathrm{}}^{\mathrm{}}e^{\frac{\mathrm{}^2\kappa ^2}{2\mu 𝒯}}𝑑\kappa `$. We can easily obtain the specific heat per atom is

$$c_v=\frac{𝒯^2}{𝒩𝒯^2}=\frac{1}{2}+\frac{E_{nm}^2E_{nm}^2}{𝒯^2},$$
(13)
where $`E_{nm}=\frac{_{mn}E_{mn}e^{E_{mn}/𝒯}}{_{mn}e^{E_{mn}/𝒯}}`$ and $`E_{nm}^2=\frac{_{mn}E_{mn}^2e^{E_{mn}/𝒯}}{_{mn}e^{E_{mn}/𝒯}}`$.

In Fig.4, the symbols “$`\mathrm{}`$” reflect $`c_v`$ varying with the temperature $`𝒯`$, which implies Ne atoms inside (20,20) tube behave as 3D gas at high temperature (specific heat approaches to 3/2) and 2D gas at low temperature (specific heat is 1). Therefore, we can naturally assume that all atoms in the valley of potential $`U(r)`$ at low temperature, i.e. lie on the shell $`S^{}`$ with radius $`\varrho =\rho [1(5/2)^{1/6}\eta ]`$.

At high temperature, the Ne atoms maybe evaporate from the nanotube, but our model cannot recover this effect because our potential is infinite near the surface of carbon nanotute.

## 6  Specific Heats at low temperature

Now we consider the thermal property of Ne inside the nanotube at low temperature in detail. We assume that all atoms lie on the shell $`S^{}`$ with radius $`\varrho =\rho [1(5/2)^{1/6}\eta ]`$ at temperature lower than some critical temperature $`T_{lc}`$ which depends on the energy difference between the tow lowest energy of radial excitations. From Eqs.(5)-(7) we can easily calculate the potential $`U_s(v,z)`$ on $`S^{}`$, where $`v=\varrho \theta `$. The Hamiltonian of single particle can be expressed as $`H^{}=\frac{\mathrm{}^2}{2\mu }(\frac{^2}{v^2}+\frac{^2}{z^2})+U_s(v,z).`$

In fact, $`U_s`$ has periodic structure. If we denote $`𝜶_1=(2\pi \varrho /N,\tau )`$, $`𝜶_2=(0,T)`$ and $`𝜸_l=l_1𝜶_1+l_2𝜶_2`$, we have $`U_s(𝐫+𝜸_l)=U_s(𝐫)`$, where $`l_1,l_2`$ and $`𝐫=(v,z)`$. On the one hand, we have the Bloch’s theorem :

$$\begin{array}{c}H^{}\mathrm{\Phi }(𝜿,𝐫)=E_𝜿\mathrm{\Phi }(𝜿,𝐫),\hfill \\ \mathrm{\Phi }(𝜿,𝐫+𝜸_l)=e^{i𝜿𝜸_l}\mathrm{\Phi }(𝜿,𝐫),\hfill \end{array}$$
(14)
which suggests that $`\mathrm{\Phi }(𝜿,𝐫)=_ja(𝜿+𝐆_j)e^{i(𝜿+𝐆_j)𝐫}`$, where $`𝐆_j=j_1𝜷_1+j_2𝜷_2`$ with $`j_1,j_2`$, $`𝜷_1=(N/\varrho ,0)`$ and $`𝜷_2=(\tau N/(T\varrho ),2\pi /T)`$. From Eq.(30), we obtain the secular equation

$$det(_{lj}E_𝜿\delta _{lj})=0,$$
(15)
where $`_{lj}=\frac{\mathrm{}^2}{2\mu }(𝜿+𝐆_j)^2\delta _{lj}+𝒰_{lj}`$, $`𝒰_{lj}=\frac{1}{\mathrm{\Omega }_0}_{\mathrm{\Omega }_0}e^{i(𝐆_l𝐆_j)𝐫}U_s(𝐫)𝑑𝐫`$ and $`\mathrm{\Omega }_0=|𝜶_1\times 𝜶_2|`$.

On the other hand, Periodic boundary condition along the circumference of the shell $`S^{}`$ suggests that we just need to consider the first Brillouin zone which consists of $`𝜿=(m_v/\varrho ,\kappa _z)`$ where $`m_v,0m_v<N`$ and $`\kappa _z,0\kappa _z<2\pi /T`$.

From Eq.(15), we can calculate the energy $`E_{m_v,\kappa _z}`$, and then the free energy $`=𝒩𝒯\mathrm{ln}𝒵_1`$ , where $`𝒵_1=_{m_v}_0^{2\pi /T}e^{E_{m_v,\kappa _z}/𝒯}𝑑\kappa _z`$. Moreover, the specific heat per atom is

$$c_v=\frac{𝒯^2}{𝒩𝒯^2}=\frac{E^2E^{}^2}{𝒯^2},$$
(16)
where $`E^{}=\underset{m_v=0}{\overset{N1}{}}_0^{2\pi /T}E_{m_v,\kappa _z}e^{E_{m_v,\kappa _z}/𝒯}𝑑\kappa _z/𝒵`$ and $`E^2=\underset{m_v=0}{\overset{N1}{}}_0^{2\pi /T}E_{m_v,\kappa _z}^2e^{E_{m_v,\kappa _z}/𝒯}𝑑\kappa _z/𝒵`$.

In Fig.5, the symbols “$`\mathrm{}`$” reflect $`c_v`$ varying with the temperature $`𝒯`$, which implies Ne atoms inside (20,20) tube behave as the lattice gas  at ultra low temperature (specific heat is 0) and 2D gas at low temperature (specific heat approaches 1). There is no 1D gas inside (20,20) tube, which is quite different from our usual notion.

## 7 why not use CSM?

At the sufficient high temperature, $``$ converge to the $`_{cl}`$. If a carbon nanotube, for example (20, 20) tube, just a geometric tube with diameter $`d2`$ nm, the characteristic temperature $`𝒯_c=h^2/(2\pi md^2)=0.04`$ K. We just need to use the CSM if $`𝒯>>𝒯_c`$. However, here the carbon nanotube provides potential to Ne. In the potential, it is difficult matter to calculate the characteristic temperature is not a simple matter. For example, consider some atoms in a harmonic potential $`U(r)=\frac{1}{2}m\omega ^2r^2`$. We can easily calculate the specific heat per atom $`c_v=(3\omega ^2/𝒯^2)\mathrm{exp}(\omega /𝒯)/[1\mathrm{exp}(\omega /𝒯)]^2`$ . Here we have set $`\mathrm{}=1`$. We know $`c_v=3`$ if $`𝒯>>\omega `$ and this is the classical case. Thus we can set $`𝒯_c=\omega `$ (In the whole form $`𝒯_c=\mathrm{}\omega /k_B`$). Above discussion suggests that $`𝒯_c`$ depends on the potential. We propose a coarse method to obtain $`𝒯_c`$ which is the least root of the following equation:

$$\{\begin{array}{c}U_m+𝒯_c=U(R_1)=U(R_2)\hfill \\ R_2R_1=h/\sqrt{2\pi m𝒯_c}\hfill \end{array}.$$
(17)
The physical meanings of parameters in Eq.(17) are shown in the Fig.6.

For the potential $`U(r)=\frac{1}{2}m\omega ^2r^2`$ in above example, $`U_m=0,R_1=R_2,R_2R_1=2R_2`$, using Eq.(17) we arrive at $`𝒯_c=\frac{\sqrt{\pi }}{2}\omega \omega `$ (the whole form is $`𝒯_c=\frac{\sqrt{\pi }\mathrm{}\omega }{2k_B}\mathrm{}\omega /k_B`$). For (20,20) nanotube, the approximate potential is $`U(r)=4\epsilon [(\frac{\stackrel{~}{\sigma }}{\rho r})^{10}(\frac{\stackrel{~}{\sigma }}{\rho r})^4]`$, with $`\rho =13.56`$ Å, $`\epsilon =390`$ K, and $`\stackrel{~}{\sigma }=2.63`$ Å. We can numerical solve Eq.(17) and give $`𝒯_c60`$ K. If $`𝒯>>𝒯_c`$ (e.g. $`𝒯>3𝒯_c=180`$ K), the classical case is expected and CSM is valid. But CSM can not be used if $`𝒯𝒯_c`$ (e.g. $`𝒯<2𝒯_c=120`$ K), especially $`𝒯<𝒯_c=60`$ K.

In Fig.4 we give results obtained from OS and CSM. The information shown in Fig.4 agrees with our coarse estimate to $`𝒯_c`$ very well. Therefore, using the CMS is insufficient. It is necessary to consider the quantum mechanics. In other words, quantum mechanics must be used when studying the behaviors of gases in the nanometer scale.

## 8 discussion and conclusion

We notice that the potential $`U`$ is a function of $`(r,\theta ,z)`$ and $`U(r,\theta +2\pi j_1/N,z+j_1\tau +j_2T)=U(r,\theta ,z)`$, where $`N,\tau ,T`$ are parameters of a carbon nanotube, and $`j_1,j_2`$ are two integers. Under this symmetry, we must use a generalized Bloch’s theorem (see appendix) to find the eigenvalues of Schrödinger equation. We can also obtain the secular equation of energies. Although we have 16 Pentium IV computers, but the data and computing time exceed our computers’ ability if we calculate the specific heat with temperature varying from 0 to 400 K. But, fortunately, we find that the variation of potential with $`\theta ,z`$ is less than 30 K, much smaller than that with $`r`$. In Sec.5, we neglect the effect of the variation of potential with $`\theta ,z`$. Indeed, this effect can be neglected under high temperature in terms of intuition. At high temperature, the values of $`c_v`$ approach to that calculated from CSM. We show them in Fig.7, in which the triangles are values neglected the effect of $`\theta ,z`$ and the line is the result calculated from CSM considering the effect of $`\theta ,z`$. It suggests our intuition is reliable.

We know the effect of $`\theta ,z`$ is important under the low temperature. In Sec.6, we consider this effect in detail. Because the shape of $`U`$ is very sharp and deep with r, the particles are confined in the bottom of $`U(r)`$ at low temperature. Thus we can consider the effect of $`\theta ,z`$ but neglect the effect of $`r`$ at low temperature.

Therefore, we believe that we have grasped the main factors under at high and low temperature and our qualitative conclusions are reliable.

The interactions (include the contribution from the polarization of the pi-electrons) between Ne and C atoms are intermolecular forces, van der Waals type interactions. As the theoretical study, the simplest form of van der Waals force is the Lenard-Jones potential. We believe that the eventual results based on the Lenard-Jones potential are qualitatively reasonable. Maybe, many researchers  use it because of the same reasons. Indeed, the more detailed study needs to consider the Gay-Berne potential  or other complex potentials. We are considering to overcome this difficult problem by using numerical simulations.

It is necessary to notice that we just consider the dilute Ne and the interactions between Ne atoms are naturally neglected. Therefore we do not describe the phase transition of Ne in SWNT’s. If we intend to deal with the phase transition, we must include the interactions between Ne atoms and use the methods of quantum field theory, such as the renormalization group approach , etc.

In conclusion, we calculate the specific heat of Ne in (20,20) SWNT at low and high temperature and find the dimensional crossover of thermal behavior in this system. Especially, the dilute Ne gas exhibits as 2D behavior at low temperature or 2D lattice behavior at ultra low temperature. The simple physical picture is that all atoms are confined on a shell blow a critical temperature.

We are grateful to kind discussion of Prof. H. W. Peng, Dr. R. An, M. Li, Y. Zhang, and L. R. Dai.

## Appendix: Generalized Bloch’s theorem

In Sec.5, we calculate the specific heat with neglecting the effect of $`\theta ,z`$; and in Sec.6, in fact, we neglect the effect of $`r`$. In other words, we do not obtain the specific heats in a consistent way from low to high temperature. Indeed, the potential $`U`$ is a function of $`(r,\theta ,z)`$ and satisfies $`U(r,\theta +2\pi j_1/N,z+j_1\tau +j_2T)=U(r,\theta ,z)`$, where $`j_1,j_2`$ are two integers. If a potential satisfies this condition, it is called the spiral symmetric. Under this spiral symmetry, we must use a generalized Bloch’s theorem to find the eigenvalues of Schrödinger equation and then calculate the specific heat in a consistent way .

Above all, we go over traditional translational symmetry and Bloch’s theorem .

The hamiltonian of a system with the translational symmetry is expressed as

$$H=\frac{\mathrm{}^2}{2\mu }^2+V(𝐫);V(𝐫+𝐑_j)=V(𝐫),$$
(18)
where $`𝐑_j=n_{j1}𝐚_1+n_{j2}𝐚_2+n_{j3}𝐚_3`$ and $`n_{j1},n_{j2},n_{j3}`$ are integers.

Define translational operators $`𝒥(𝐑_j)`$, which act on a function $`f(𝐫)`$ as:

$$𝒥(𝐑_j)f(𝐫)=f(𝐫+𝐑_j).$$
(19)
It follows that

$$𝒥(𝐑_j)𝒥(𝐑_l)=𝒥(𝐑_l+𝐑_j)=𝒥(𝐑_j+𝐑_l)=𝒥(𝐑_l)𝒥(𝐑_j),$$
(20)
and

$`𝒥(𝐑_j)Hf(𝐫)={\displaystyle \frac{\mathrm{}^2}{2\mu }}𝒥(𝐑_j)^2f(𝐫)+𝒥(𝐑_j)V(𝐫)f(𝐫)`$
$`={\displaystyle \frac{\mathrm{}^2}{2\mu }}^2𝒥(𝐑_j)f(𝐫)+V(𝐫+𝐑_j)f(𝐫+𝐑_j)`$
$`=[{\displaystyle \frac{\mathrm{}^2}{2\mu }}^2+V(𝐫)]f(𝐫+𝐑_j)=H𝒥(𝐑_j)f(𝐫).`$ (21)
Therefore $`\{𝒥(𝐑_j),H\}`$ is the set of conserved quantities. In this case, an eigenfunction of the Hamiltonian must be an eigenfunction of the translational operators, i.e., $`𝒥(𝐑_j)\psi (𝐫)=\psi (𝐫+𝐑_j)=\lambda (𝐑_j)\psi (𝐫)`$ if $`H\psi (𝐫)=E\psi (𝐫)`$. Furthermore, the electron density must be periodic, i.e., $`|\psi (𝐫+𝐑_j)|^2=|\psi (𝐫)|^2`$. It follows that

$$|\lambda (𝐑_j)|^2=1.$$
(22)
But from Eq.(20), we know

$$\lambda (𝐑_j)\lambda (𝐑_l)=\lambda (𝐑_j+𝐑_l).$$
(23)
The solution of Eq.(23) under the constraint Eq.(22) is $`\lambda (𝐑_j)=e^{i𝜿𝐑_j}`$. Thus we have Bloch’s theorem:

$$\psi (𝐫+𝐑_j)=e^{i𝜿𝐑_j}\psi (𝐫).$$
(24)
Next, we will turn to the spiral symmetry and generalize the Bloch’s theorem. The hamiltonian with the spiral symmetry in cylindrical coordinates is expressed:

$$H=\frac{\mathrm{}^2}{2\mu }(\frac{^2}{r^2}+\frac{}{rr}+\frac{^2}{r^2\theta ^2}+\frac{^2}{z^2})+V(r,\theta ,z),$$
(25)
$$V(r,\theta ,z)=V(r,\theta +j_1\vartheta ,z+j_1\tau +j_2T),$$
(26)
where $`\vartheta =2\pi /N,N\tau =MT`$, $`N,M`$, $`j_1,j_2`$. Define vectors $`𝐫=(\theta ,z)`$, $`𝐑_j=(j_1\vartheta ,j_1\tau +j_2T)`$ in the space $`[0,2\pi )\times `$, and operators $`𝒥(𝐑_j)`$ which act on a function $`f(r;𝐫)`$ as $`𝒥(𝐑_j)f(r;𝐫)=f(r;𝐫+𝐑_j)`$. It follows that $`𝒥(𝐑_j)𝒥(𝐑_l)=𝒥(𝐑_l+𝐑_j)=𝒥(𝐑_j+𝐑_l)=𝒥(𝐑_l)𝒥(𝐑_j)`$ and $`𝒥(𝐑_j)V(r;𝐫)=V(r;𝐫)𝒥(𝐑_j)`$. Otherwise,

$$𝒥(𝐑_j)(\frac{^2}{r^2\theta ^2}+\frac{^2}{z^2})f(r;𝐫)=(\frac{^2}{r^2\theta ^2}+\frac{^2}{z^2})𝒥(𝐑_j)f(r;𝐫).$$
(27)
Thus

$$𝒥(𝐑_j)H=H𝒥(𝐑_j).$$
(28)
We can obtain generalized Bloch’s theorem analogizing the method to obtain the traditional Bloch’s theorem:

$$\psi (r;𝐫+𝐑_j)=e^{i𝜿𝐑_j}\psi (r;𝐫),$$
(29)
where $`𝜿=(l,\kappa )`$, $`l=0,1,\mathrm{},N1,0\kappa <2\pi /T`$.

If we set $`𝜶_1=(\vartheta ,\tau )`$, $`𝜶_2=(0,T)`$, $`𝜷_1=(N,0)`$, $`𝜷_2=(\tau N/T,2\pi /T)`$, $`𝐆_𝐣=j_1𝜷_1+j_2𝜷_2=(j_1Nj_2N\tau /T,2\pi j_2/T)`$, then we can expand the wave function as the superposition of planar waves:

$$\psi (r;𝐫)=\underset{𝐣}{}C_𝐣\phi (r)e^{i(𝜿+𝐆_𝐣)𝐫}.$$
(30)
In cylindrical coordinates, Eqs.(29) and (30) are expressed

$`\psi (r,\theta +{\displaystyle \frac{2\pi j_1}{N}},z+j_1\tau +j_2T)=e^{i[\frac{2\pi lj_1}{N}+\kappa (j_1\tau +j_2T)]}\psi (r,\theta ,z),`$ (31)
$`\psi _{l\kappa }(r,\theta ,z)={\displaystyle \underset{j_1j_2}{}}C_{j_1j_2}\phi _{l\kappa }(r)e^{i[(l+j_1Nj_2N\tau /T)\theta +(\kappa +2\pi j_2/T)z]}.`$ (32)
If Expanding $`\phi _{l\kappa }(r)`$ with normalized orthogonal basis $`\{\chi _{nl\kappa }(r),n\}`$ we obtain the secular equation from Eq.(32) and Shrödinger equation

$$det[(_{m_1m_2j_1j_2nn^{}}E_{nl\kappa }\delta _{j_1m_1}\delta _{j_2m_2}\delta _{nn^{}})]=0,$$
(33)
where $`_{j_1j_2m_1m_2nn^{}}=𝒯_{j_1j_2nn^{}}\delta _{j_1m_1}\delta _{j_2m_2}+𝒰_{j_1j_2m_1m_2nn^{}}`$, $`𝒯_{j_1j_2nn^{}}=\frac{\mathrm{}^2}{2\mu }_0^\rho \chi _{nl\kappa }(r)[\frac{d^2}{dr^2}+\frac{d}{rdr}\frac{(l+j_1N\frac{j_2N\tau }{T})^2}{r^2}(\kappa +\frac{2\pi j_2}{T})^2]\chi _{n^{}l\kappa }(r)r𝑑r`$, and
$`𝒰_{j_1j_2m_1m_2nn^{}}=\frac{N}{2\pi T}_0^{2\pi }𝑑\theta _0^T𝑑z_0^\rho r𝑑r\chi _{nl\kappa }(r)U(r,\theta ,z)\chi _{n^{}l\kappa }(r)e^{i[(j_1m_1)\frac{(j_2m_2)\tau }{T}]N\theta +\frac{i2\pi (j_2m_2)z}{T}}`$.

## References


text:
# Studies of Mass and Size Effects in Three-Dimensional Vibrofluidized Granular Mixtures

## I Introduction

Granular systems exhibit many properties that are different from systems composed of elastic particles. For example, driven granular systems display standing and traveling waves Douday et al. (1989); Pak and Behringer (1993); Melo et al. (1995), oscillons Umbanhowar et al. (1996), heaping, and convection Evesque and Rajchenbach (1989); Ehrichs et al. (1995). In addition, granular mixtures show size segregation Rosato et al. (1987) and steady-state kinetic energies that are not equal for each component in the mixture Losert et al. (1999). This departure from equipartition is not unexpected, but it is one of the most striking differences between elastic and inelastic systems. Understanding the properties of mixtures is particularly important for granular systems since, unlike molecular systems, they are never completely monodisperse.

Theoretical studies of granular systems have focused on two distinct classes. One consists of systems that are not driven or heated. The initial energy decays over time as a result of inelastic collisions. During this “cooling” process there is a period during which the density is homogeneous. Several workers have presented kinetic theories McNamara and Luding (1998); Garzo and Dufty (1999); Dahl et al. (2002) and mean-field theories based on Maxwell models Marconi and Puglisi (2002a) to describe the properties of mixed granular systems in this homogenous cooling state. In the other class of systems, an energy source, such as a vibrating wall, is present. This leads to a non-equilibrium steady state that has been studied by several researchers McNamara and Luding (1998); Marconi and Puglisi (2002b); Barrat and Trizac (2002a, b). In both classes the components are predicted to have different kinetic energies, or granular temperatures, that depend on the mass, size, and restitution coefficient of the grains Dahl et al. (2002); Barrat and Trizac (2002a).

Recently, two-dimensional Feitosa and Menon (2001) and three-dimensional systems Losert et al. (1999); Wildman and Parker (2002) of driven, granular mixtures have been studied experimentally. Losert *et al.*, who first reported the difference in granular temperature Losert et al. (1999), observed that the velocity distributions deviated from those observed in systems with elastic collisions. Feitosa and Menon studied density distributions and granular temperature profiles in two-dimensional systems with and without gravity Feitosa and Menon (2001). Wildman and Parker have studied the convection patterns, density distributions, and temperature profiles in three-dimensional systems Wildman and Parker (2002). These studies determined that the heavier particles are at a higher granular temperature than the lighter particles. In both two and three-dimensional systems, the ratio of the temperatures varies as the relative proportion of the heavy and light particles is changed. The temperature ratio, however, is independent of the inelasticity of the particles Feitosa and Menon (2001); Wildman and Parker (2002).

While the experimental techniques employed in the studies cited above have provided many useful insights into granular behavior, they cannot easily isolate the effects of particle mass, size, and inelasticity. Theoretical and computational methods are useful in this respect. Molecular dynamics simulations of granular mixtures can accurately reproduce the phenomena observed in experiment Talbot and Viot (2002); Paolotti et al. , while providing information on the effects of the individual properties mentioned above. For example, Paolotti *et al.* Paolotti et al.  and Barrat and Trizac Barrat and Trizac (2001) investigated the effects of rotation, mass ratio, and relative density in two-dimensional vibrated systems. Mixtures have also been studied under two-dimensional shear flow conditions Clelland and Hrenya (2002). Simulations of the homogeneous cooling state in two-dimensional systems are consistent with experiment and theory Douday et al. (1989); Barrat and Trizac (2001); Clelland and Hrenya (2002).

It is important to stress that the conclusions drawn from simulations of two-dimensional systems cannot necessarily be extended to three-dimensions. In particular, the system boundaries have a much larger influence in three-dimensions, as recently demonstrated by Talbot and Viot Talbot and Viot (2002).

This paper presents a three-dimensional, event-driven, molecular dynamics simulation of a mixture of inelastic hard-spheres. The simulation methodology is discussed in Section II followed by a comparison to the available experimental results, Section III.1. Finally, the effects of isolated changes in the mass ratio (Section III.2) and size ratio (Section III.3) on the energy distribution and component temperatures are examined.

## II Model and Simulation

The three-dimensional system Talbot and Viot (2002) (Figure 1) contains a mixture of inelastic hard spheres in an infinitely tall cylinder of radius $`R`$ under the influence of gravity. The mixture is composed of $`n_1`$ spheres of mass $`m_1`$ and diameter $`d_1`$ and $`n_2`$ spheres of mass $`m_2`$ and diameter $`d_2`$. Energy is injected into the system by means of the base of the cylinder, which vibrates in a symmetric saw-tooth waveform of amplitude $`A`$ and frequency $`\nu `$.

Three kinds of collisions occur in the system: particle-particle, particle-wall and particle-base. The post-collisional velocities ($`𝐯_{\alpha ,𝐢}^{}`$ and $`𝐯_{\beta ,𝐣}^{}`$) resulting from a collision between a particle from component $`\alpha `$ and a particle from component $`\beta `$ with masses $`m_\alpha `$ and $`m_\beta `$ and initial velocities $`𝐯_{\alpha ,𝐢}`$ and $`𝐯_{\beta ,𝐣}`$, respectively, are given by

$$𝐯_{\alpha ,𝐢}^{}=𝐯_{\alpha ,𝐢}\frac{m_\beta }{m_\alpha +m_\beta }(1+c)[(𝐯_{\alpha ,𝐢}𝐯_{\beta ,𝐣})\widehat{𝐧}]\widehat{𝐧}$$
(1a)
$$𝐯_{\beta ,𝐣}^{}=𝐯_{\beta ,𝐣}+\frac{m_\alpha }{m_\alpha +m_\beta }(1+c)[(𝐯_{\alpha ,𝐢}𝐯_{\beta ,𝐣})\widehat{𝐧}]\widehat{𝐧}$$
(1b)
The unit vector between the centers of the colliding particles is $`\widehat{𝐧}`$, and $`c`$ is the appropriate restitution coefficient. While Equations 1a and 1b conserve momentum, they imply an energy loss of

$$\mathrm{\Delta }E=\frac{1}{2}\mu (1c^2)[(𝐯_{\alpha ,𝐢}𝐯_{\beta ,𝐣})\widehat{𝐧}]^2.$$
(2)
where $`\mu =m_\alpha m_\beta /\left(m_\alpha +m_\beta \right)`$ is the reduced mass of the particles involved in the collision.

In a binary mixture it is generally necessary to specify three restitution coefficients for particle-particle collisions: $`c_{11}`$ and $`c_{22}`$ for intra-component collisions and $`c_{12}(=c_{21})`$ for inter-component collisions. Several authors have reported on random Barrat and Trizac (2002b); Barrat et al. (2001) and velocity-dependent Brilliantov et al. (1996); Goldman et al. (1998) restitution coefficients. Luding and McNamara have proposed a contact time model that also leads to a variable restitution coefficient Luding and McNamara (1998). For simplicity, we have chosen to use a constant value for the restitution coefficients. We further assume that $`c_{11}=c_{22}=c_{12}=c`$.

Particle-wall collisions are governed by

$$𝐯_{\alpha ,𝐢}^{}=𝐯_{\alpha ,𝐢}(1+c_{\alpha ,w})(𝐯_{\alpha ,𝐢}\widehat{𝐫})\widehat{𝐫}$$
(3)
where $`c_{\alpha ,w}`$ is the appropriate restitution coefficient for component $`\alpha `$, and $`\widehat{𝐫}`$ is the radial unit vector. We assume that the restitution coefficient for collision with the wall is constant for both species and $`c_{1,w}=c_{2,w}=c_w`$. Particle-base collisions are governed by

$$𝐯_{\alpha ,𝐢}^{}=𝐯_{\alpha ,𝐢}(1+c_{\alpha ,b})[(𝐯_{\alpha ,𝐢}𝐯_𝐰)\widehat{𝐤}]\widehat{𝐤}$$
(4)
where $`c_{\alpha ,b}`$ is the appropriate restitution coefficient for component $`\alpha `$, $`𝐯_𝐰`$ is the velocity of the base at the instant of collision, and $`\widehat{𝐤}`$ is the unit vector in the z-direction. The restitution coefficients for collisions with the base are also assumed to be constant and equal for both species (i.e. $`c_{1,b}=c_{2,b}=c_b`$).

A phenomenon similar to inelastic collapse can be observed in these simulations. For certain ranges of the velocity a given particle will collide repeatedly with the side wall. As its energy is dissipated, the particle approaches the wall ever more closely. This is accompanied by an increase in collision frequency that eventually “freezes” the simulation. To prevent this phenomenon from occurring, a small impulse is imparted to the particle toward the center of the cylinder once its radial velocity falls below a certain value. This method has been used previously Talbot and Viot (2002), and the threshold value was set such that the injected energy does not discernibly influence the simulation output. It is also possible for a particle to come to rest on the base for a time corresponding to $`1/2`$ a cycle of the vibration. To avoid this possibility, the sign of the z-component of the velocity for the particle is inverted when the velocity of the colliding particle is found to match the velocity of the base. This condition was found to occur once every $`12.5\times 10^6`$ collisions for the parameter values used in this paper. Thus, this method causes little perturbation in the simulation output.

We calculated a number of properties from the particle positions and velocities generated by the simulation. The packing fraction $`\eta _\alpha `$ for component $`\alpha `$ is defined as

$$\eta _\alpha =\frac{n_\alpha v_\alpha }{V}$$
(5)
where $`n_\alpha `$ is the number of particles of component $`\alpha `$ in the volume element $`V`$, and $`v_\alpha =\pi d_\alpha ^3/6`$ is the volume of a particle of this component. Another property that we calculated is the kinetic energy or granular temperature, $`T_\alpha `$, of each component using the following equations.

$`T_{\alpha ,x}=m_\alpha v_{\alpha ,x}^2`$
$`T_{\alpha ,y}=m_\alpha v_{\alpha ,y}^2`$ (6a)
$`T_{\alpha ,z}=m_\alpha v_{\alpha ,z}^2`$
$$T_\alpha =\frac{(T_{\alpha ,x}+T_{\alpha ,y}+T_{\alpha ,z})}{3}=\frac{m_\alpha v_\alpha ^2}{3}$$
(6b)
The “partial” temperatures in the x-,y-, and z-directions are $`T_{\alpha ,x}`$, $`T_{\alpha ,y}`$, $`T_{\alpha ,z}`$, and the angular brackets denote a time average over all particles of component $`\alpha `$.

Our first objective was to model the experimental system studied by Wildman and Parker Wildman and Parker (2002). We chose simulation parameters that correspond to those of the experiment, i.e., a cylinder of diameter 145mm that is shaken at 50 Hz with an amplitude of 1.74 mm. The acceleration due to gravity is taken as $`g=9.81m/s^2`$. The restitution coefficient for particle-particle collisions is $`c=0.91`$, for particle-wall collisions is $`c_w=0.68`$, and for particle-base collisions is $`c_b=0.88`$.

We performed simulations in which we varied the relative proportions of large and small particles while maintaining enough particles to cover the base of the cylinder with a monolayer for comparison to the experimental work of Wildman and Parker Wildman and Parker (2002). We then performed additional simulations to examine the effect of varying the mass ratio $`m_2/m_1`$ and the size ratio $`d_2/d_1`$ independently. In the discussion of these three studies, component 2 will always refer to the smaller and/or lighter component.

The particles of each component were randomly placed in the cylinder with random velocities. We then equilibrated all systems for approximately 5000 collisions per particle. Data were collected over $`2.4\times 10^4`$ collisions per particle at intervals of approximately 10 collisions per particle. These data were then averaged for each component to obtain a representation of the system at steady state.

## III Results and Discussion

### III.1 Systems with varying composition

For the simulations discussed in this section, the size ratio of the two components was set at $`d_2/d_1=0.8`$ and the mass ratio was set at $`m_2/m_1=0.512`$. Systems with the following compositions were then studied: $`n_1=525`$, $`n_2=270`$; $`n_1=350`$, $`n_2=540`$; and $`n_1=175`$, $`n_2=810`$.

#### III.1.1 Velocity Field

The velocity fields of components 1 and 2 for a system with $`n_1=525`$ and $`n_2=270`$ particles are shown in Figures 2a and 2b, respectively. It can be seen that both components circulate in a pattern that rises in the center of the cylinder and falls at the walls. This convection pattern has been observed previously in both experiment and simulation Wildman and Parker (2002); Talbot and Viot (2002); Wildman et al. (2001a). The patterns shown here are very similar to those reported by Wildman and Parker Wildman and Parker (2002) with the center of the convection for both species present at a radius of approximately 50 mm and a height of approximately 40 mm.

#### III.1.2 Packing Fraction

The packing fraction, as a function of radius and height for each component, is shown in Figure 3. It can be seen that there is a density gradient in both the radial and vertical directions for both components, just as is observed in monodisperse systems Talbot and Viot (2002); Wildman et al. (2001a). The data from the simulation are qualitatively similar to those observed in experiment except that the simulation shows a higher concentration at the bottom of the cylinder near the wall for both components. It should be noted that while the maximum density occurs at the same point for both components, component 2 obtains a greater height than component 1. This trend is opposite that observed for size segregation in weakly tapped systems where the larger particles rise above the smaller ones Rosato et al. (1987).

The radially averaged packing fractions as a function of height are shown in Figure 4. The three curves correspond to the three relative fractions of component 1 and component 2. In all cases, the packing fraction of both components increases steeply at small heights, reaches a maximum and decays at large heights. The packing fraction of component 1 decreases and it increases for component 2 as the relative amount of each component is changed. The changes in the relative fractions only affects the magnitudes of the packing fractions. There is no noticeable variation in the details of the packing fraction profiles (i.e. the position of the maximum, the rate of decay, etc.) as the relative amounts of the components are changed.

#### III.1.3 Temperature

The granular temperature of each component is also studied as a function of the relative proportions. Figure 5 shows contour plots of the granular temperature of component 1 in the x- (Figure 5a), y- (Figure 5b), and z- (Figure 5c) directions for the system presented in Figure 2. The symmetry evident in the x- and y-directions is produced by the unbiased introduction of energy into these directions by particle-particle collisions. These two partial temperatures decay rapidly in both the radial and vertical dimensions from a maximum near the center of the cylinder, close to the vibrating base. The z-direction, however, is different because of the bias introduced by the vibrating base. This partial temperature decays in the vertical dimension with very little radial dependence at small heights. At larger heights, the center of the cylinder is slightly warmer than the surrounding area, as would be expected given the toroidal flow profile presented in Figure 2.

Figure 6 shows the height dependence of the radially averaged, partial granular temperatures of each component. It can be seen that the height profiles of the temperature are similar for both components at all relative proportions. The temperatures in the x- and y-directions show the formation of a maximum for both components as height increases. These maxima occur very close to the height that corresponds to the maximum in packing fraction. The increase in the partial temperatures can then be attributed to an increase of particle-particle collisions that inject energy into the x- and y-directions, increasing the corresponding temperatures. Figure 6 also shows that the temperature in the z-direction is larger than that in the other two directions. The minimum observed as height increases has been predicted by Brey *et al.* for systems in which the particles do not interact with a top barrier, but are under the influence of gravity Brey et al. (2001). Wildman *et al.* observed the minimum the experimental systems Wildman et al. (2001a, b), and Ramirez and Soto presented a hydrodynamic theory that addresses this phenomenon Ramirez and Soto (2002).

Figures 6a and 6b also show that the temperatures of the two components change as the relative proportion of the two components changes. The temperature in the three spatial directions decreases as the relative proportion of the larger particles is decreased. As the fraction of component 1 decreases, the number of collisions with the smaller particles increases, causing the temperature of the larger particles to decrease. Also, the temperature of the entire system decreases because component 2 does not gain as much kinetic energy from the base since it is lighter than component 1. The decrease in the temperature, however, does not change the height at which the extrema in $`T_x`$, $`T_y`$ and $`T_z`$ are observed. Only the magnitude of the measured temperature appears to change in these systems (Figure 6a and 6b).

While the trends observed in the temperature of the two components are similar, there are differences between the two components. The temperatures of component 1 (Figure 6a) are greater than those of component 2 (Figure 6b) in all the systems, an effect particularly pronounced in the z-direction. The temperatures in the x- and y-directions differ only slightly. It is difficult to determine whether the differences in temperature between the two components are dominated by the differences in size or mass from these data. Thus, we conducted further simulations to determine the individual effects of mass and size.

### III.2 Systems with varying mass ratio

To determine the effects of particle mass, we simulated systems with mass ratios $`m_2/m_1=`$ 0.01, 0.125, 0.25, 0.5, and 1.0 at constant size, $`d_2/d_1=1.0`$, and relative fraction $`n_2/n_1=1.0`$ and with a total $`n_1+n_2=1050`$ particles in the system.

#### III.2.1 Energy Distribution

We expect changes in the mass ratio to affect the exchange of energy, and hence the temperature, of each component. We therefore computed the average change in energy, $`\mathrm{\Delta }E_\alpha _\beta `$, of a particle of component $`\alpha `$ resulting from collisions with particles of component $`\beta `$. We examined the energy exchanges resulting from particle-particle, particle-wall and particle-base collisions. The data we collected from the simulations for each kind of collision are presented in Figure 7.

As expected, both components lose energy on collision with the wall, while they experience a net energy gain on collision with the base (7a). Both components also lose energy from intra-component, particle-particle collisions (squares in 7b). The inter-component, particle-particle collisions (circles in 7b), however, show trends that are not intuitively obvious. Component 1 shows a loss of energy for all mass ratios, while component 2 shows that there may be a loss or gain of energy depending on the mass ratio. We obtained a theoretical estimate for this quantity by assuming each component has a Maxwell-Boltzmann velocity distribution but with a temperature specific to the component. The details of the derivation are presented in the Appendix. The average energy loss for a particle of component $`\alpha `$ resulting from collisions with particles of component $`\beta `$ and average component temperatures (kinetic energies) of $`T_\alpha `$ and $`T_\beta `$ is

$$\mathrm{\Delta }E_\alpha _\beta =k_B\mu \left(1+c\right)\left(\left(1+c\right)\frac{T_\beta m_\alpha +T_\alpha m_\beta }{m_\alpha \left(m_\alpha +m_\beta \right)}2\frac{T_\alpha }{m_\alpha }\right)$$
(7)
where $`k_B`$ is the Boltzmann constant. The total average energy loss per collision between particles of components $`\alpha `$ and $`\beta `$ with temperatures of $`T_\alpha `$ and $`T_\beta `$ is

$$\mathrm{\Delta }E_\alpha _\beta +\mathrm{\Delta }E_\beta _\alpha =k_B\left(1c^2\right)\frac{T_\beta m_\alpha +T_\alpha m_\beta }{m_\alpha +m_\beta }$$
(8)
Equation 8 indicates that for inelastic collisions, $`c<1`$, the total energy of the colliding pair always decreases, even if the temperatures of the two components are different. For equal temperatures, Equation 7 shows that this is also true for the individual energies of the components. If, however, the temperatures are different, it is possible for the energy of the lighter component to *increase* on average due to collisions with the heavier component.

We calculated values of $`\mathrm{\Delta }E_\alpha _\beta `$ for inter-component and intra-component, particle-particle collisions for each component using Equation 7. The masses of the particles and the restitution coefficient are set by the input parameters, but the temperatures of each component are not known *a priori*. Therefore, the temperature was obtained from the simulation output for each mass ratio. The results from Equation 7 are presented in Figure 7b along with the results obtained from the simulation. It can be seen that the average energy changes calculated for the two kinds of particle-particle collisions compare favorably with those obtained from the simulation. There is better agreement for intra-component collisions because there is no temperature or mass difference between the colliding particles. The predicted energy changes for inter-component collisions show the same trend that is observed in the simulation results, but the magnitude of the change is incorrect. Specifically, the equation overestimates the change in energy resulting from inter-component, particle-particle collisions. This error probably arises from the inhomogeneities in the particle density and temperature (see Figures 3 and 5), which are not accounted for in the Maxwell-Boltzmann distribution.

More generally, we note that as the mass ratio decreases, the energy change associated with any collision also decreases. This variation in the energy change will affect the bulk properties observed for these systems, such as the packing fraction and the temperature discussed below.

#### III.2.2 Packing Fraction

The packing fraction of each component as a function of height is presented in Figure 8. The general behavior is the same for both components, and is similar to that already discussed in Section III.1.2. The packing fraction of component 1 varies little with mass ratio changes. In particular, the maximum density is at the same height for all the systems examined. At greater heights, however, component 1 condenses and then expands as the mass of component 2 decreases. This phenomenon is more clearly visible in the insert of Figure 8a, which shows a linear relationship between $`\mathrm{ln}(\eta )`$ and the height. At large altitudes, the slope of the lines in the insert increase as component 1 condenses and decrease as component 1 expands in the system. This behavior corresponds to the increased energy loss that is observed for inter-component particle-particle collisions (see Figure 7).

The packing fraction profiles of component 2 undergo a much more significant change as the mass ratio is decreased. Specifically, Figure 8b indicates a steady depletion of this component from around the maximum, with a compensating increase at large altitudes, as the particles expand into the upper reaches of the cylinder. For mass ratios of 0.25 and below, two local maxima are present. These are most distinct for the systems in which component 2 has a net gain in energy due to collisions with particles of component 1 (mass ratios of 0.125 and 0.01). Thus, the two maxima are formed as the particles of component 2 try to separate from component 1. The increase in the energy forces the particles of component 2 toward the base and toward higher altitudes. This is what is observed in the packing fractions shown in Figure 8b with one maximum very close to the base, and one maximum that increases in altitude as the mass ratio decreases. The plots of $`\mathrm{ln}(\eta )`$, shown in the insert, display a steady decrease in the slope as the mass ratio decreases.

#### III.2.3 Temperature

The effect of mass ratio on the temperature in the z-direction for components 1 and 2 is shown in Figures 9a and 9b, respectively. The minimum in the temperature is obvious for both particles. It is also easily seen that the temperature in the z direction of component 1 goes through a minimum as the mass ratio of the two components decreases. The temperatures in the x- and y-directions (not shown) also follow the same trend. This indicates that total temperature for component 1 goes through a minimum as the mass ratio is decreased. The changes in the temperature coincide with the changes observed in the packing fraction (Figure 8a) and the energy changes for the different kinds of collisions (Figure 7 ). This implies that the changes in the velocities of the particles of component 1 affect the temperature, just as expected from Equations 6. As shown by Brey *et al.* Brey et al. (2001) and Warr *et al.* Warr et al. (1995), a relationship exists between the temperature and the packing fraction in a single component system. Extending their results to a multi-component system, we obtain

$$\frac{d\mathrm{ln}(\eta _\alpha )}{dz}\frac{m_\alpha g}{k_BT_\alpha }$$
(9)
This equation holds for high altitudes and restitution coefficients close to 1. Thus, a limiting temperature can be calculated for each system using the data presented in the insert of Figure 8a. The results, plotted as horizontal lines in Figure 9b, correspond well with the asymptotic temperatures. This indicates that the decaying edge of the packing fraction is a good indicator of the temperature at those altitudes.

We observed trends for component 2 that are very different from those just discussed for component 1. Figure 9b shows that the temperature decreases as the mass ratio decreases. The decrease is expected since temperature is directly related to the mass of the particle (equations (6)). Thus, the lighter particles will have a lower temperature than the heavier particles. Component 2 exhibits a minimum as the height increases, just as in the case of component 1. However, the minimum becomes shallower as the temperature decreases. Figure 9b also shows the temperature obtained from the packing fraction using Equation 9. Again, we find that there is good agreement between the temperature calculated by Equation 9 and that calculated from Equations 6 for large heights.

### III.3 Effect of varying size ratio

Finally, we studied the effects of particle size by simulating systems with size ratios of $`d_2/d_1=`$ 1.0, 0.8, 0.5, and 0.1. The mass ratio was held constant at $`m_2/m_1=1.0`$, and the relative fraction was held constant at $`n_2/n_1=1.0`$, with $`n_1+n_2=1050`$. The changes observed in the energy distributions, the packing fraction and the partial temperatures are presented and discussed below.

#### III.3.1 Energy Distribution

Changes in particle size results in changes in the mean free path Brey et al. (2001) and the pair correlation function at contact Garzo and Dufty (1999); Barrat and Trizac (2002a, b). These changes affect the particle velocities by changing the number of collisions that a particle experiences in a given amount of time. Thus, changes in the size ratio are expected to result in changes in the distribution of the energy, just as observed for the mass ratio.

Figures 10a and 10b show the average energy changes that particles of components 1 and 2 experience as a result of collision. Collisions with the wall cause a loss of energy for both components while collision with base increase the energy of the particles for all size ratios (Figure 10a). The collisions between particles, however show different trends than seen in the mass ratio study. As seen in Figure 10b, the inter-component collisions do not result in an energy loss for component 1 for all size ratios. The smaller the size ratio, the greater the amount of energy injected per collision into component 1 from collisions with component 2. All particle-particle collisions decrease the energy of component 2. The energy loss per collision for for both inter-component and intra-component, particle-particle collisions increase as the size ratio decreases. It is interesting to note that the energy loss due to both inter-component and intra-component collisions are the same for component 2 at the smallest size ratio. The energy changes for all but the smallest size ratio agree well with those predicted by Equation 7. The smallest size ratio shows a large deviation between the energy loss predicted by the equations and that determined from the simulation. One reason for the discrepancy is the assumption of Maxwell-Boltzmann velocities used in determining the equations. The differences in the distribution of energy within each component and between the two components that we observe here will affect the packing fraction and temperature for these systems.

#### III.3.2 Packing Fraction

Figure 11 shows the effect of the size ratio on the packing fraction of the system. We can see that component 1 reaches higher altitudes as the diameter of component 2 decreases without any noticeable shift in the position of the maximum in packing fraction. This implies that the particles of component 1 expand through the system as the size ratio decreases. There are two possible causes of the increase in the tail of the packing fraction at large heights. First, the particles of component 1 are able to retain more energy because of a decrease in the collision rate between particles as the size ratio decreases. The decrease in the collision rates is the result of a decrease in the total excluded volume of the system as the size of component 2 is reduced (deduced from Figures 11a and 11b). The change in the total excluded volume of the system reduces the amount of energy lost to particle-particle collisions. The other cause of the increased altitude is the energy gain that comes from collisions with component 2. Since energy can be gained at positions above the base, the particles will be able to travel to higher altitudes in the system.

Figure 11b shows the effect of the size ratio on component 2. We can see an overall reduction in the packing fraction of component 2 as the particle size is reduced. This occurs because the number of particles of component 2 is held fixed as the size is decreased. We also see that the particles of component 2 are able to reach greater altitudes in the system as the size ratio decreases (see the insert in Figure 11b). The increase in the altitude is the result of the decrease in the number of collisions discussed above.

#### III.3.3 Temperature

Figure 12a shows the z-component of the granular temperature as a function of height for component 1 for the four size ratios. Figure 12b shows the same for component 2 of the mixture. A minimum is again observed in the temperatures of each system for each component. The other general features of the temperature profiles in the x- and y-directions, while not shown here, are the same as seen and discussed in connection with Figure 6. It can be easily seen in Figure 12 that the temperature of both the large and small particles increase as the size ratio decreases. This increase in temperature results from a decrease in the energy lost due to particle-particle collisions. The reduction in the collision rate as the size ratio decreases can be observed in Table 1. It is interesting that the energy loss due to intra-component collisions increases for component 2 as the collision rate for intra-component collisions decreases (Figure 10b). In addition, there is a discrepancy between the theory and the simulation results observed in Figure 10 at small size ratios that should be noted. The deviations are the result of the systems being dominated by collisions with the wall and not particle-particle collisions (see Table 1), as assumed in the theory (see the Appendix).

We used Equation 9 and the high altitude tails of the packing fractions shown in the inserts in Figures 11a and 11b to calculate a temperature. The temperatures calculated in this manner are also plotted in Figure 12. The temperature calculated from Equation 9 coincides with the temperature in the z-direction at large heights, just as in the case of the mass ratio. The agreement between the two temperatures is very good for both components for each size ratio, except the smallest. This may indicate that the packing fraction is not a very sensitive measure of the temperature when the difference in size between the components is large. The difference in the two temperatures may also be indicative of the changes in the dominant processes in the energy distribution. For example, the smallest size ratio is not dominated by particle-particle collisions like the other systems. The edge effects associated with the walls of the cylinder become more important, but they are not considered in the theory used to obtain Equation 9.

## IV Conclusions

We have shown that, like the single component system Talbot and Viot (2002), the simulation results reproduce the phenomena observed in experimental studies. Specifically, the experimentally observed Wildman and Parker (2002) flow pattern, radial dependence of the packing fraction and temperature of the two components are qualitatively reproduced by the simulation.

There are a number of possible reasons for the quantitative differences between the simulation and experiments. For example, for each component the simulated particles are identical spheres, while in the experimental system there is a small distribution of shape, size, and mass. The simulation also assumes that the particles are frictionless with constant restitution coefficients, which is not the case for ballotini glass spheres used in the experiments. The saw-tooth wave form of the vibrating base is an idealization of the sinusoidal vibration of the experiment. It is thought, however, that this does not have a large influence on the system behavior McNamara and Barrat (1997). In any case, we feel confident that our model captures the key physical aspects of the experimental system and can therefore be used to study the influence of various system parameters.

We also studied the effects of mass and size ratio in this paper. Generally, we observed that changes in either ratio do not result in any segregation of the particles. The lighter particles attain greater heights than the heavier particles in the mass ratio studies. For both components the rate of decay of the packing fraction at large heights is a good indication of the granular temperature in that region Brey et al. (2001); Warr et al. (1995).

As the mass ratio decreases, the overall temperature of the system decreases. This is consistent with the lower amount of energy gained by the lighter component from the vibrating base and an overall lowering of the efficiency of energy transfer as the mass ratio decreases. The packing fraction and temperature of the individual components, however, have a non-trivial dependence on the mass ratio. For the heavier component, both these quantities exhibit a minimum as the mass ratio decreases. We also observed that the energy changes due to inter-component collisions exhibit extrema for each component as the mass ratio decreases. Specifically, the heavier component shows a minimum, whereas the lighter component shows a maximum in the energy change. The lighter component actually gains energy from collisions with the other component around a mass ratio of $`m_2/m_1=0.5`$. The energy changes due to intra-component collisions for both components, however, are negative for all mass ratios. The energy loss due to intra-particle collisions increases for the heavy component as the mass ratio decreases, while it decreases for the lighter component. For comparison, we developed an approximate theory to calculate the energy changes for particle-particle collisions by assuming that the particle velocities follow a Maxwell-Boltzmann distribution with a component specific temperature. While this assumption is not strictly correct given the inhomogeneities in density and temperature, the theory is qualitatively accurate for the inter-component collisions and in near quantitative agreement with the simulation for the intra-component collisions.

As the size ratio decreases at constant mass ratio, the overall temperature of the system increases and the total packing fraction decreases. The particle-particle and particle-boundary collision rates decrease and increase, respectively, as the size ratio decreased. At the same time, the larger particles begin to gain energy from collisions with the smaller component. The approximate theory of the energy changes is again able to reproduce qualitatively the observed trends. However, the agreement between the theory and simulation results worsens as the size ratio decreases.

Wildman and Parker Wildman and Parker (2002) observed a decrease in the temperature as they decreased the ratio of the number of large to small ballotini spheres. Our results show that this effect is dominated by the difference in mass of the two components, and not the difference in size.

###### Acknowledgements.

We thank Ricky Wildman for useful discussions and the National Science Foundation (CHE-9814236) for financial support.

*

## Appendix A Energy dissipation equations

The energy change of particle $`i`$ of component $`\alpha `$ resulting from collision with particle $`j`$ of component $`\beta `$ is

$`2m_\alpha \mathrm{\Delta }E_{\alpha ,i}=\mu ^2\left(1+c\right)^2\left(𝐯_{\mathrm{𝐢𝐣}}\widehat{𝐧}\right)^2m_\alpha \mu \left(1+c\right)`$
$`\times \left[\left(𝐯_{\alpha ,𝐢}\widehat{𝐧}\right)^2\left(𝐯_{\beta ,𝐣}\widehat{𝐧}\right)^2+\left(𝐯_{\mathrm{𝐢𝐣}}\widehat{𝐧}\right)^2\right]`$ (10)
where $`𝐯_{\mathrm{𝐢𝐣}}=𝐯_{\alpha ,𝐢}𝐯_{\beta ,𝐣}`$ and $`\mu =m_\alpha m_\beta /(m_\alpha +m_\beta )`$. The total energy loss due to the collision is simply $`\mathrm{\Delta }E_{total}=\mathrm{\Delta }E_{\alpha ,i}+\mathrm{\Delta }E_{\beta ,j}`$, or

$$\mathrm{\Delta }E_{total}=\frac{1}{2}\mu \left(1c^2\right)\left(𝐯_{\mathrm{𝐢𝐣}}\widehat{𝐧}\right)^2$$
(11)
These equations are exact for each collision that occurs in the system. In order to determine the average values, however, some assumptions must be made. Specifically, we assume that the velocities of each component are described by a Maxwell-Boltzmann distribution that is homogeneous, isotropic and characterized by a component-specific temperature $`T_\alpha `$:

$$f\left(𝐯_{\alpha ,𝐢}\right)d𝐯_{\alpha ,𝐢}=\left(\frac{2\pi m_\alpha }{k_BT_\alpha }\right)^{3/2}\mathrm{exp}\left(\frac{m_\alpha v_{\alpha ,i}^2}{2k_BT_\alpha }\right)d𝐯_{\alpha ,𝐢}$$
(12)
In order to average over the velocities $`𝐯_{\alpha ,𝐢}`$, $`𝐯_{\beta ,𝐣}`$, and $`𝐯_{\mathrm{𝐢𝐣}}`$ appearing in Equations A and 11, we introduce center-of-mass, $`𝐯_𝐜`$, and relative $`𝐯_𝐫`$ velocities:

$`𝐯_{\alpha ,𝐢}={\displaystyle \frac{m_\beta /T_\beta }{m_\alpha /T_\alpha +m_\beta /T_\beta }}𝐯_𝐫+𝐯_𝐜`$ (13)
$`𝐯_{\beta ,𝐣}={\displaystyle \frac{m_\alpha /T_\alpha }{m_\alpha /T_\alpha +m_\beta /T_\beta }}𝐯_𝐫+𝐯_𝐜`$ (14)
$`𝐯_{\mathrm{𝐢𝐣}}=𝐯_{\alpha ,𝐢}𝐯_{\beta ,𝐣}=𝐯_𝐫`$ (15)
In this coordinate system the energy change for particle $`i`$ becomes

$`2m_\alpha \mathrm{\Delta }E_{\alpha ,i}=`$ $`\mu ^2\left(1+c\right)^2\left(𝐯_𝐫\widehat{𝐧}\right)^2`$
$`2T_\alpha \mu ^{}\mu \left(1+c\right)\left(𝐯_𝐫\widehat{𝐧}\right)^2`$
$`2m_\alpha \mu \left(1+c\right)\left(𝐯_𝐫\widehat{𝐧}\right)\left(𝐯_𝐜\widehat{𝐧}\right)`$
where

$$\mu ^{}=\frac{(m_\alpha /T_\alpha )(m_\beta /T_\beta )}{m_\alpha /T_\alpha +m_\beta /T_\beta }$$
(17)
The total energy loss (shown in Equation 11) becomes

$$\mathrm{\Delta }E_{total}=\frac{1}{2}\mu \left(1c^2\right)\left(𝐯_𝐫\widehat{𝐧}\right)^2$$
(18)
in the new coordinate system.

It is now necessary to average Equations A and 18 over the fraction of collisions between a particle of component $`\alpha `$ with a particle of component $`\beta `$ with a relative velocity between $`v_r`$ and $`v_r+dv_r`$. Straightforward modification of the standard kinetic theory result McQuarrie (1976) gives:

$$p\left(v_r\right)dv_r=\frac{1}{2}\left(\frac{\mu ^{}}{k_B}\right)^2v_r^3\mathrm{exp}\left(\frac{\mu ^{}v_r^2}{2k_B}\right)dv_r$$
(19)
We then compute

$$\left(𝐯_𝐫\widehat{𝐧}\right)^2=_0^{\mathrm{}}𝑑v_rp(v_r)_0^\sigma 𝑑bh(b)\left(𝐯_𝐫\widehat{𝐧}\right)^2$$
(20)
where $`h(b)db=(2b/\sigma ^2)db`$ is the probability that the impact parameter lies between $`b`$ and $`b+db`$. Substituting $`\left(𝐯_𝐫\widehat{𝐧}\right)^2=v_r^2(1(b/r)^2)`$ and evaluating the integrals leads to

$$\left(𝐯_𝐫\widehat{𝐧}\right)^2=\frac{2k_B}{\mu ^{}}$$
(21)
Using this result and the fact that $`\left(𝐯_𝐫\widehat{𝐧}\right)\left(𝐯_𝐜\widehat{𝐧}\right)=0`$ gives Equations 7 and 8.