Size segregation in granular media induced by phase transition %A purely thermodynamic mechanism for size segregation in granular media
M. Tarzia$^{a}$
A. Fierro$^{a,b}$
M. Nicodemi$^{a,b}$
M. Pica Ciamarra$^{a}$
and A. Coniglio$^{a,b}$
In order to study analytically the nature of the size segregation in granular mixtures, we introduce a mean field theory in the framework of a statistical mechanics approach, based on Edwards' original ideas. For simplicity we apply the theory to a lattice model for hard sphere binary mixture under gravity, and we find a new purely thermodynamic mechanism which gives rise to the size segregation phenomenon. By varying the number of small grains and the mass ratio, we find a crossover from Brasil nut to reverse Brasil nut effect, which becomes a true phase transition when the number of small grains is larger then a critical value. We suggest that this transition is induced by the effective attraction between large grains due to the presence of small ones (depletion force). Finally the theoretical results are confirmed by numerical simulations of the $3d$ system under taps. \pacs{45.70.Mg, 64.75.+g, 05.50.+q}
paper.tex
\documentclass[twocolumn,prl,showpacs]{revtex4}
\usepackage{epsfig}% Include figure files
\usepackage{graphicx}
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
%\nofiles

\def\lan{\langle}
\def\ran{\rangle}
\def\e{\mbox{e}}

\begin{document}

\title{Size segregation in granular media induced by phase transition
%A purely thermodynamic mechanism for size segregation in granular media
}
\author{M. Tarzia$^{a}$, A. Fierro$^{a,b}$, M. Nicodemi$^{a,b}$, M. Pica
Ciamarra$^{a}$, and A. Coniglio$^{a,b}$}
\affiliation{${}^a$ Dipartimento di Scienze Fisiche, Universit\`{a} degli Studi
di Napoli ``Federico II'', INFM and INFN, via Cinthia, 80126 Napoli, Italy}

\affiliation{${}^b$ INFM - Coherentia, Napoli, Italy}

\date{\today}

\begin{abstract}
In order to study analytically the nature of the size segregation in granular
mixtures, we introduce a mean field theory in the framework of a
statistical mechanics approach, based on Edwards' original ideas. For simplicity
we apply the theory to a lattice model for hard sphere binary mixture under
gravity, and we find
a new purely thermodynamic mechanism which gives rise to
the size segregation phenomenon.
By varying the number of small grains and the mass ratio, we find a crossover 
from Brasil nut to reverse Brasil nut effect, which becomes a true phase
transition when the number of small grains is larger then a critical value.
We suggest that this transition is induced by the effective attraction between
large grains due to the presence of small ones (depletion force).
Finally the theoretical results are confirmed by numerical simulations of
the $3d$ system under taps.
\pacs{45.70.Mg, 64.75.+g, 05.50.+q}
\end{abstract}

\maketitle

In the last decades a great attention has been devoted to the
study of the problem of vertically shaken granular mixtures under
gravity. It was observed that such systems depending on the
control parameters can mix or segregate their components
spontaneously according to a mechanism which is still largely
unclear, although of deep practical and conceptual relevance
\cite{rev_segr,nagel}. Rosato et al. \cite{rosato} showed that
large grains surrounded by a sea of smaller ones rise to the top
when subjected to vertical shaking. This is the well known
``Brazil Nut Effect'' (BNE), while the opposite one (i.e. large
grains on the bottom and small ones on the top) is known as
``Reverse Brazil Nut Effect'' (RBNE) \cite{luding,breu}. The BNE
was originally explained \cite{rosato} in terms of a geometric
effect, based on  ``percolation'' arguments.
%In a vibrated bed smaller grains are
%more likely to find a void through which they can
%percolate down to the bottom \cite{rosato,bridgewater}.
Along with geometry, dynamical effects,
such as convection \cite{knight} or  inertia \cite{shinbrot},
were shown to play a crucial role.
The hydrodynamic equations for binary mixtures \cite{jenkins} were also applied,
and in this framework the regimes of BNE and RBNE \cite{herrmann} were found
to depend on the presence of inelastic dissipation.
%in order to investigate this phenomenon, and in this framework
%the size segregation was explained \cite{herrmann} in terms of competition
%between buoyancy and geometric forces.

Recent results finally outlined that segregation processes may involve
global mechanisms, such as condensation \cite{luding} or
phase separation \cite{mullin}.
This suggested a change of perspective on the issue and
the idea to formulate a statistical mechanics description of
these phenomena \cite{luding,both,seg}.
To this aim, in Ref.s~\cite{luding, both} the interplay between
size and mass was studied and a phase diagram for BNE/RBNE transition
based on the competition between percolation and condensation was proposed.
%Recently, an
%long range
%effective attraction \cite{sanders}
%between two large grains in a $2d$ vertically vibrated granular bed of smaller
%ones was also shown to play a role.

In the present paper, applying a statistical mechanics approach for
non-thermal systems \cite{Edwards, e1}, a new purely thermodynamic
mechanism is found responsible for size segregation. We follow a
statistical mechanics approach for granular materials under taps,
%(where the energy is pumped into the system in pulses, and the
%measures are performed when the system is at rest)
which was recently developed \cite{nfc} on the basis of Edwards'
original ideas \cite{Edwards, e1}. This approach postulates that
time averages coincide with suitable ensemble averages over the
mechanically stable states (i.e. those where the system is found
at rest). Such hypothesis was shown \cite{nfc} to hold with good
approximation in a lattice model for granular mixtures under taps.
In order to study size segregation, in this framework we apply a mean field 
theory to such model.

%In particular in Ref.~\cite{seg}
%a phase diagram for BNE/RBNE transition was
%obtained in the Bethe approximation.
%the partition function of the model considered in \cite{nfc}
%was solved in the Bethe approximation.
%A phase diagram for BNE/RBNE transition was given as a
%function of the configurational temperatures and the mass ratio
%of the two species.
%The calculations were limited to the
%high shaking amplitude region, while the low shaking amplitude
%region, where a real phase transition is present, is not
%considered at all.

We find that, by fixing the number of grains of the two species, it is
possible to obtain a smooth crossover from BNE to RBNE by
changing the mass ratio of the grains. However if the number of
small grains increases above a critical value, a sharp
transition occurs resulting in a phase separation between the two
species: Vertical size segregation is induced by the presence of gravity
which forces the heavier phase to move downward. At the critical
point, critical fluctuations are found typical of a second order
phase transition. We suggest that  this demixing
phase separation is due to the presence of the depletion force.
This is an effective attraction between large grains, obtained
by tracing over all the degrees of freedom of small grains. The
strength of the interaction is proportional to the number of small
grains. The depletion force, typically present in thermal
binary mixtures \cite{Evans}, in fact may be also present in $3d$
granular mixtures under gravity \cite{depletion}. The
theoretical results are confirmed by numerical simulations of the
$3d$ system.

%We find that by varying suitably the control parameters a genuine
%phase separation between the two species occurs; Size segregation
%is induced by the presence of gravity which forces the heavier
%phase to move downward. We suggest that the depletion force (an
%effective attraction between large particles obtained by tracing
%over all the degrees of freedom of small particles) typically
%present in thermal binary mixtures \cite{Evans}, might be also
%present in $3d$ granular mixture under gravity \cite{depletion},
%and induce such a demixing phase separation. Finally the
%theoretical results are compared with numerical simulations of the
%$3d$ system.
{\em The model~}-- The model \cite{nfc,seg} is a hard sphere binary mixture 
under gravity
made up of two species, $1$ (small) and $2$ (large), respectively with
diameter $a_0=1$ and $\sqrt{2} a_0$, and masses $m_1$ and $m_2$.
In order to simplify the calculations we divide the space in
cubic cells of linear size $a_0$ whose vertices can be or not
occupied by grains. Each site of the lattice is labeled by the two indices
$i=1,\dots,L^2$ and $z=1,\dots,L_z$, where $L$ and $L_z$ are respectively the
horizontal and vertical dimensions of the lattice.
A microstate of the model is completely identified by the occupancy variables
$\{n_i^z\}$, associated to each site of the lattice:
$n_i^z= 0, 1, 2$
respectively if the site $(i,z)$ is empty, filled by a small grain or by a
large one.
In terms of these variables the model Hamiltonian is
\begin{equation}
{\cal H}={\cal H}_{HC}+ m_1gH_1 + m_2gH_2,
\label{H1}
\end{equation}
where $H_1=\sum_{i,z} z \delta_{n_i^z 1}$,
$H_2=\sum_{i,z} z \delta_{n_i^z 2} $ are respectively the heights of the two
species,
and ${\cal H}_{HC}$ is the hard core potential, preventing two nearest
neighbor sites to be both occupied if at least one contains a large grain.

In Ref.s~\cite{nfc} it was shown that the model, Eq.~(\ref{H1}), under taps is
described with good approximation by a statistical mechanics approach based
on Edwards' originally ideas. In particular it was found that
the weight of a microscopic state, $\{n_i^z\}$, is given by
\begin{equation}
P(\{n_i^z\})= \Pi (\{n_i^z\})\cdot e^{-{\cal H}_{HC}(\{n_i^z\})-
%\beta_1m_1gH_1(\{n_i^z\}) -\beta_2m_2gH_2(\{n_i^z\})},
\beta_1m_1gH_1 -\beta_2m_2gH_2},
\end{equation}
where $\Pi (\{n_i^z\})$ is a projector over the
``mechanically stable'' states, namely $\Pi (\{n_i^z\})=1$ if $\{n_i^z\}$
is mechanically stable, and $\Pi (\{n_i^z\}) =0$ otherwise,
$\beta_1$ and $\beta_2$ are the variables conjugated respectively to the
gravitational energies of the two species, and
the two ``configurational temperatures'' $T_1=\beta_1^{-1}$ and
$T_2=\beta_2^{-1}$ are increasing functions of the tap amplitude.
The system partition function is given by ${\cal Z}=\sum_{\{n_i^z\}}
P(\{n_i^z\})$
%\begin{equation}
%\nonumber
%{\cal Z} =  \sum_{\{n_i^z\}} \Pi \cdot \e^{-\left[{\cal H}_{HC}  + \beta_1
%m_1gH_1 + \beta_2 m_2gH_2 \right]},
%\label{Zed}
%\end{equation}
where the sum is over all the microstates but, due to the projector,
only the mechanically stable ones are taken into account.
%In numerical simulations of the model, Eq. (\ref{H1}), under taps we
%found \cite{prep} that the configurational temperatures,
%$T^{conf}_1 \equiv \beta_1^{-1}$ and $T^{conf}_2\equiv \beta_2^{-1}$,
%increase as function of the tap amplitude, and coincide each other in the
%limit of low tap amplitudes. Thus for simplicity in the following
%we put $\beta_{conf} \equiv \beta_1= \beta_2$.
In the following calculations we adopt a simple definition of ``mechanical
stability'': a grain is considered stable if it has a grain underneath
\cite{nota_proj}.

{\em Mean field approximation~}-- 
In the present paper the partition function, ${\cal Z}$,
is evaluated  in mean field approximation:
we consider a generalization of Bethe-Peierls method for anisotropic systems
(due to gravity) already used in previous papers \cite{seg,epl,jamming}.
In particular, we consider a Bethe lattice made up by
$L_z$ horizontal layers (i.e., $z\in\{1,...,L_z\}$) whose edges can be
occupied by grains. Each layer is a random graph of given
connectivity, $k-1=4$. Each site in layer $z$ is also
connected to its homologous site in $z-1$ and $z+1$
(the total connectivity is thus $k+1=6$).
We use the grancanonical ensemble by introducing two chemical potentials,
$\mu_1$ and $\mu_2$, relative to species $1$ and $2$.

In the following we give the main ideas of the calculations and
refer to a longer paper \cite{jamming} for the details (see also
\cite{MP}). The Bethe-Peierls recursion equations, which allow to calculate
the partition function, are written in terms of
the local ``cavity fields'' defined by:
$e^{s_{n}^{(i,z)}} = Z_{n,s}^{(i,z)}/Z_{0,s}^{(i,z)}$,
$e^{u_{n}^{(i,z)}} = Z_{n,u}^{(i,z)}/Z_{0,u}^{(i,z)}$,
$e^{v_{n}^{(i,z)}} = \overline{Z}_{n,u}^{(i,z)}/Z_{0,u}^{(i,z)}$,
$e^{d_{n}^{(i,z)}} = Z_{n,d}^{(i,z)}/Z_{0,d}^{(i,z)}$,
$e^{c_{n}^{(i,z)}} = \overline{Z}_{n,d}^{(i,z)}/Z_{0,d}^{(i,z)}$
(with $n=1,2$).
Here $Z_{0,s}^{(i,z)}$ and $Z_{n,s}^{(i,z)}$ are the partition functions
of the ``side'' branch restricted respectively to configurations
where the site $(i,z)$ is empty or filled by a particle of specie $n$;
analogously, $Z_{n,u}^{(i,z)}$, $Z_{0,u}^{(i,z)}$ and
$\overline{Z}_{n,u}^{(i,z)}$ (resp. $Z_{n,d}^{(i,z)}$, $Z_{0,d}^{(i,z)}$
and $\overline{Z}_{n,d}^{(i,z)}$) are the partition functions of the
``up'' (resp. ``down'') branch restricted to configurations where
the site $(i,z)$ is filled by a grain of species $n$, empty with the upper
(resp. lower) site empty and empty with the upper (resp. lower) site
filled by a grain of species $n$.

The calculations are limited to the ``fluid'' phase corresponding
to a solution of Bethe-Peierls equations
where local fields in each layer are site independent \cite{nota_cryst}.
Such a solution, characterized by horizontal translational invariance, is
given by the fixed points of the following equations:
\begin{eqnarray} \label{ricorrenza}
&&\!\!\!\!\!\!\!\!\!\!\!\!e^{s_{n}^{(z)}}  =
%e^{\beta_n(\mu_n-m_ngz)}
\phi_n^{(z)}
\left[
%\frac{A_n^{(z)}}{S^{(z)}}
{\cal K}_n^{(z)}
\right]^{k-2}
\frac{P_n^{(z+1)} R_n^{(z-1)}}{Q^{(z)}},
\\
\nonumber
&&\!\!\!\!\!\!\!\!\!\!\!\!e^{u_{n}^{(z)}} =
\phi_n^{(z)}
%e^{\beta_n(\mu_n-m_ngz)}
\left[
%\frac{A_n^{(z)}}{S^{(z)}}
{\cal K}_n^{(z)}
\right]^{k-1} P_n^{(z+1)},
%\\
%\nonumber
\,\,\,\,\,\,\,\,
e^{d_{n}^{(z)}}  = e^{u_{n}^{(z)}} \frac{R_n^{(z-1)}}{P_n^{(z-1)}},
%\phi_n^{(z)}
%e^{\beta_n(\mu_n-m_ngz)}
%\left[
%\frac{A_n^{(z)}}{S^{(z)}}
%{\cal K}_n^{(z)}
%\right]^{k-1} R_n^{(z-1)},
\\
\nonumber
&&\!\!\!\!\!\!\!\!\!\!\!\!e^{v_{n}^{(z)}}=e^{u_{n}^{(z+1)}},
\qquad \,\,\,\,\,\,\,
e^{c_{n}^{(z)}}  =  e^{d_{n}^{(z-1)}} \left[1 + R_2^{(z-1)}
\right]^{-1},
\end{eqnarray}
where $n=1,2$, and $\phi_n^{(z)} = e^{\beta_n(\mu_n-m_ngz)}$.
The functions appearing in Eq.~(\ref{ricorrenza}) have a simple expression
in terms of the local fields \cite{nota_func}.
From the local fields we calculate the free energy \cite{nota_free}, $F$, and
all the quantities of interest, such as
density profiles of small and large grains, $\sigma_1 (z)$ and $\sigma_2 (z)$,
the grain number per unit surface, $N_{n} = \sum_z \sigma_n (z)$, and their
average height,
$h_n=\langle{z_n}\rangle=\sum_z z \sigma_n(z) / \sum_z \sigma_n(z)$
(with $n=1,2$).
\begin{figure}[ht]
%\vspace{-0.05cm}
\begin{center}
\epsfig{figure=deh.ps,scale=0.3,angle=270}
\end{center}
\vspace{-0.3cm}
\caption{(color online). {\bf Main frame} 
$\frac{\Delta h}{h}$ as a function of $\delta$ in
the mean field approximation for fixed value of $N_2=~0.685$, and $N_1=~0.5$ 
(blue curve), $N_1=~0.9$ (green curve), $N_1=~N_{1c}=~1.37$ (red curve) 
and $N_1=~2.4$ 
(black curve) ($N_1$ increases along the arrow). 
{\bf Inset:} The susceptibility, $\chi$, as a function of $\delta$ for fixed 
value of $N_2=~0.685$, and
$N_1=~0.5$ (blue curve), $N_1=~0.9$ (red curve) and $N_1=~N_{1c}=~1.37$ 
(black curve)  ($N_1$ increases from bottom to top).}
\label{fig:deh}
%\vspace{-0.3cm}
\end{figure}

In the following instead of using the chemical potential variables, $\mu_1$
and $\mu_2$, we use the conjugate variables, $N_1$ and $N_2$, respectively
the number of the small and large grains per unit surface. For simplicity we
fix the configurational temperatures, $T_1=~T_2=~1$, 
and the number of large grains, $N_2$, and vary $N_1$.
In order to study the BNE/RBNE transition we also vary the masses of the 
two species, $m_1$ and $m_2$ (keeping constant $m_1+m_2$), and measure
the height difference, $ \frac{\Delta h}{h} \equiv \frac{h_1 - h_2}{h_1 + h_2}$:
This quantity, that is thus a function of $N_1$ and $\delta \equiv \frac{2 m_1 
- m_2}{2 m_1 + m_2}$, is positive when small grains are in average above 
large grains (RBNE) and negative otherwise (BNE).

\begin{figure}[ht]
%\vspace{-0.05cm}
\begin{center}
\epsfig{figure=maxsusc.ps,scale=0.3,angle=270}
\end{center}
\vspace{-0.3cm}
\caption{(color online).
The coexistence curve (black curve on the left) and the maximum of the 
susceptibility (blue curve on the right) as a function of $N_1^{-1}$.} 
\label{fig:susc} 
%\vspace{-0.1cm}
\end{figure}
In the Main Frame of Fig.~\ref{fig:deh},
$ \frac{\Delta h}{h}(\delta,N_1)$ is plotted for fixed $N_1$ as a function
of $\delta$: The different curves correspond to increasing values of $N_1$
(along the arrow). 
We find that for low enough values of $N_1$,
$\frac{\Delta h}{h}$
is a continuous monotonous decreasing function of $\delta$, due to entropic
and buoyancy effect.
By increasing $N_1$ above some critical value ($N_{1c} \simeq
1.37$) $\frac{\Delta h}{h}$ instead exhibits first order transition with a
finite jump from a positive value to a negative one at 
$\delta=\delta_{coex}(N_1)$, where a BNE phase and a RBNE one
coexist with the same free energy. As shown in Fig.~\ref{fig:susc},
the jump becomes smaller and smaller  approaching 
$N_{1c}$ from above, and goes to zero continuously at $N_{1c}$ (the critical 
point is $(\delta_c, N_{1c})$, where $\delta_c\equiv \lim_{N_1\rightarrow 
N_{1c}} \delta_{coex}(N_1)$).

Along the coexistence curve near the critical point in both phases
the order parameter $\frac{\Delta h}{h}\simeq (N_{1c}^{-1}-N_1^{-1})^{\beta}$,
where, solving the mean field equations numerically, $\beta$ is found
consistent with the mean field value $\beta=~0.5$.

As in any second order phase transition, the fluctuations of the order 
parameter diverge at the critical point.
In the Inset of Fig.~\ref{fig:deh}, 
the susceptibility, $\chi(\delta,N_1)\equiv N_{tot}(\langle\Delta h
^2\rangle-\langle\Delta h\rangle^2)$ (where $N_{tot}$ is the total number of
grains) is plotted as a function
of $\delta$ (the curves, from bottom to top, correspond to increasing values 
of $N_1$).  As $N_1$ approaches $N_{1c}$,
the susceptibility diverges at $\delta = \delta_c$.
For $N_1<N_{1c}$, $\chi(\delta,N_1)$ displays a maximum (for the value 
of $\delta$ where $\Delta h$ is zero), which diverges as a power at $N_{1c}$,
$(N_1^{-1}-N_{1c}^{-1})^{-\gamma}$ with an exponent $\gamma$ found, solving
numerically the mean field equations, consistent with
the mean field value $\gamma=~1.0$ (see Fig.~\ref{fig:susc}).

In the present paper we have chosen as free parameters  $N_1$ and
$\delta$: $\delta$ plays the role of an external field, and $N_1$ plays the 
role of an inverse temperature, since it controls the amplitude of the
depletion force \cite{Evans}. To confirm this idea, we have studied the case of
mixture of grains with the same 
radii ($R_2/R_1 = 1$), where the depletion force is completely 
absent. We find that the critical point disappears, and a smooth crossover from
positive to negative values of $\Delta h$, varying the control parameters $N_1$
and $\delta$. 

{\em Monte Carlo simulations~}-- 
The theoretical results obtained in mean field approximation 
are confirmed by Monte Carlo simulations of the model,
Eq.~(\ref{H1}), under taps \cite{nfc,nota_tap} on a cubic lattice.
In Fig.~\ref{fig:bne}, the height difference, $\frac{\Delta h}{h}$, and the 
susceptibility, $\chi$, are plotted as function of 
$\delta$ for fixed $N_2$ and increasing values of $N_1$. 

%On the bottom the density profiles obtained for a fixed value of $\delta$ in
%the mean field approximation and
%in the simulations are plotted: A surprising agreement is found showing
%that the mean field theory gives not only qualitative
%predictions, but also quantitative.
%Two different pairs ($N_1$, $N_2$) are considered:
%In the first case (figure
%on the left) the mean field approximation predicts a continuous crossover from
%the BNE to the RBNE by decreasing $\delta$, instead in the second case
%(figure on the right) a discontinuous jump.
%($N_1=1$, $N_2=0.5$), and ($N_1=0.2$, $N_2=0.15$).
%In the simulations  we have varied $m_1$ and $m_2$ in the following
%way. Starting from positive values of $\delta$, we decrease $m_1$ and increase
%$m_2$ at a fixed rate $\dot{m}$; when $\delta$ reaches a value close to
%$-1$ we invert the rate sign. In the figure on the left
%$\dot{m}=~2\cdot 10^{-6}~taps^{-1}$,
%and in the figure on the right the orange dots and
%the red ones correspond respectively to $\dot{m}=10^{-5}~taps^{-1}$ and
%$2\cdot10^{-6}~taps^{-1}$. In the first case the mean field approximation
%predicts a continuous crossover from the BNE to the RBNE by decreasing
%$\delta$; in the simulations a similar behavior is found and the curve of
%$\Delta H/H$ as
%a function of $\delta$ is reversible. In the second case the mean field
%approximation predicts a discontinuous jump from the BNE to the RBNE
%by decreasing $\delta$;
%in the simulations
%In the simulations $m_1$ and $m_2$ are varied at a fixed rate $\dot{m}$:
%In the first case (figure on the left) we find that the equilibrium state
%is easily reached,
%the curve of $\Delta H/H$ as a function of $\delta$ is
%reversible, and a continuous crossover from BNE to RBNE is observed, analogous to
%that predicted by the mean field approximation. Differently in the second
%case (figure on the right) the equilibrium is not reached also at rate lower
%than ones here shown: The data present an evident hysteresis, which becomes
%smaller at lower rates.
%, tending to the equilibrium limit where a jump is
%expected close to $\delta=~0$.

\begin{figure}[ht]
%\vspace{-0.1cm}
\begin{center}
\epsfig{figure=deh_MC.ps,scale=0.3,angle=270}
\end{center}
\vspace{-0.3cm}
\caption{(color online). {\bf Main frame} 
$\frac{\Delta h}{h}$ as a function of $\delta$ in the Monte
Carlo simulations on the cubic lattice for fixed value of $N_2=~1.2$, and
$N_1=~0.6$ (blue curve and filled circles), $N_1=~0.8$ (green curve and empty 
circles), $N_1=~1.0$ (red curve and stars) and $N_1=~1.2$ (black curve and 
triangles). {\bf Inset} The susceptibility, $\chi$, as a function of $\delta$ 
for fixed value of $N_2=~1.2$, and $N_1=~0.6$ (blue curve and filled circles),
$N_1=~1.0$ (red curve and stars) and $N_1=~1.2$ (black curve and
triangles)}
\label{fig:bne}
%\vspace{-0.3cm}
\end{figure}


%\begin{figure}[ht]
%\vspace{-0.1cm}
%\begin{center}
%\epsfig{figure=fluct.eps,scale=0.4,angle=0}
%\end{center}
%\vspace{-0.6cm}
%\caption{(color online). The susceptibility, 
%$\chi$, as a function of $\delta$ in
%the Monte Carlo simulations on the cubic lattice for fixed value of
%$N_2=~1.2$, and $N_1=~0.6$ (blue curve and filled circles),  
%$N_1=~1.0$ (red curve and stars) and $N_1=~1.2$ (black curve and 
%triangles)}
%\label{fig:fluct}
%\vspace{-0.3cm}
%\end{figure}

As found in the mean field approximation, at fixed value of $N_1$, 
$\frac{\Delta h}{h}$ as a function of $\delta$ displays a steeper and steeper 
slope (and the susceptibility, $\chi$, a larger and larger maximum) 
as $N_1$ is increased.
A power law divergence of the maximum of the susceptibility is found at
$N_{1c}\simeq 2.2$, with an exponent $\gamma\simeq 1.6$. 
As usually in mean field
approximation, the critical point ($N_{1c}^{-1}$) is over-estimated.  

As $N_1$ is further increased, $\frac{\Delta h}{h}$ 
obtained by varying $\delta$ at a given rate, becomes irreversible, and a 
strong hysteresis appears, signaling the presence of a first order transition.


{\em Conclusions}-- We focus on the problem of the vertical size
segregation in binary granular mixtures. In the framework of a
statistical mechanics approach to granular media, based on
Edwards' original ideas, we apply a mean field theory to a simple
hard sphere lattice model under gravity, and explain BNE and RBNE
with purely thermodynamic argument. 

%Due to depletion
%force caused by small grains,
%to an effective attraction between large grains, analogous to the depletion force
%in usual thermal systems,
%we find a thermodynamic phase transition associated to the
%transition from  BNE  to RBNE, resulting in a strong size
%segregation. The theoretical results are confirmed by the
%numerical simulations of the $3d$ system.
%In principle this phase transition could be investigated
%experimentally, and its effect may be revealed above the critical point by
%the growth of the fluctuations of the height difference. 
%, and below the critical point, by the
%presence of a strong hysteresis.

We find that by varying the control parameters a transition from BNE to RBNE
can occur as a smooth crossover or as a sharp transition, depending whether
one is above or below the critical point.
The presence of such a critical point is manifested by the divergence
of the height fluctuations of the two species.
Therefore we suggest that it might be important to measure in
experiments these fluctuations in the vicinity of the transition
from BNE to RBNE. The larger the fluctuations the closer to the critical
point are the control parameters.


We would like to thank Hans Herrmann
for many interesting discussions and suggestions.
Work supported by EU Network Number  MRTN-CT-2003-504712, MIUR-PRIN 2004,
MIUR-FIRB 2001, CRdC-AMRA, INFM-PCI.

\begin{thebibliography}{40}
\bibitem{rev_segr}
J. Bridgwater, Chem. Eng. Sci. {\bf 50}, 4081 (1995).
J.M. Ottino and D.V. Khakhar,
Ann. Rev. Fluid Mech. {\bf 32}, 55 (2000).
T. Shinbrot and F.J. Muzzio, Physics Today {\bf 53}, 25 (2000).
Rosato et al., Chem. Eng. Sci. {\bf 57}, 265 (2002).
\bibitem{nagel}  M.E. Mobius,
X. Cheng, G.S. Karczmar, S.R. Nagel, H.M. Jaeger,
{\em cond-mat}/0405328.
\bibitem{rosato} A. Rosato,
K.J. Strandburg, F. Prinz,  R.H. Swendsen,
Phys. Rev. Lett. {\bf 58}, 1038 (1987).
\bibitem{luding}  D.C. Hong,
P.V. Quinn, S. Luding,
Phys. Rev. Lett. {\bf 86}, 3423 (2001).
\bibitem{breu} A. P. J. Breu,
H.-M. Ensner, C. A. Kruelle, and I. Rehberg,
Phys. Rev. Lett. {\bf 90}, 014302 (2003)

%\bibitem{bridgewater} J. Bridgwater,
%Powder Technol. {\bf 15}, 215 (1976).
%J.C. Williams,
%Powder Technol. {\bf 15}, 245 (1976).
\bibitem{knight}  J.B. Knight,
H.M. Jaeger, S.R. Nagel,
Phys. Rev. Lett. {\bf 70}, 3728 (1993).

\bibitem{shinbrot}  T. Shinbrot and F.J.  Muzzio,
Phys. Rev. Lett. {\bf 81}, 4365 (1998).
\bibitem{jenkins} J.T. Jenkins and D.K. Yoon,
Phys. Rev. Lett. {\bf 88}, 194301 (2002).
\bibitem{herrmann} L. Trujillo,
M. Alam and H. Herrmann,
Europhys. Lett. {\bf 64}, 190 (2003).
\bibitem{mullin}
P.M. Reis and T. Mullin, Phys. Rev. Lett. {\bf 89}, 244301 (2002).
P. M. Reis,
G. Ehrhardt, A. Stephenson and T. Mullin
Europhys. Lett. {\bf 66}, 357 (2004).
\bibitem{both}
J.A. Both and D.C. Hong, Phys. Rev. Lett. {\bf 88}, 124301 (2002).
\bibitem{seg}
M. Tarzia,
A. Fierro, M. Nicodemi, and A. Coniglio,
Phys. Rev. Lett. {\bf 93},
198002 (2004).

\bibitem{Edwards}  S.F. Edwards and R.B.S. Oakeshott,
Physica A {\bf 157}, 1080 (1989).
A. Mehta and S.F. Edwards, Physica A {\bf 157}, 1091 (1989).

\bibitem{e1}  M. Nicodemi, Phys. Rev. Lett. {\bf 82}, 3734 (1999).
A. Barrat,
J. Kurchan, V. Loreto and M. Sellitto,
Phys. Rev. Lett. {\bf 85}, 5034 (2000).
J.J. Brey,
A. Prados, B. S\'{a}nchez-Rey,
Physica A {\bf 275}, 310 (2000).
A. Coniglio and M. Nicodemi, Physica A {\bf 296}, 451 (2001).
D. S. Dean and A. Lef\`{e}vre, Phys. Rev. Lett. {\bf 86}, 5639 (2001).
H. A. Makse and J. Kurchan, Nature {\bf 415}, 614 (2002).
A. Fierro,
M. Nicodemi and A. Coniglio,
Europhys. Lett. {\bf 59}, 642 (2002).
G. De Smedt,
C. Godreche, J.M. Luck,
Eur. Phys. J. B {\bf 32}, 215-225 (2003).
\bibitem{nfc} M. Nicodemi,
A. Fierro, A. Coniglio,
Europhys. Lett. {\bf 60}, 684 (2002). A. Fierro,
M. Nicodemi and A. Coniglio,
Phys. Rev. E {\bf 66}, 061301 (2002). A. Fierro,
M. Nicodemi, A. Coniglio,
Jour. Phys.: Cond. Matt. {\bf 15}, S1095 (2003).

\bibitem{Evans} M. Dijkstra,
R. van Roij and R. Evans,
Phys. Rev. E {\bf 59}, 5744 (1999). M. Fasolo and P. Sollich, Phys. Rev. Lett.
{\bf 91}, 068301 (2003).

\bibitem{depletion} M. Bose,
P.R. Nott and V. Kumaran,
Europhys. Lett. {\bf 68},
508 (2004). An effective attraction
between two large grains
was also found in a $2d$ vertically vibrated granular bed of smaller ones
in Ref.  D.A. Sanders,
M.R. Swift, R.M. Bowley, and P.J. King,
Phys. Rev. Lett. {\bf 93}, 208002 (2004).

\bibitem{nota_proj} For a given grain configuration, $\{n_i^z\}$,
$\Pi (\{n_i^z\})$ has a tractable expression:
$\Pi (\{n_i^z\}) =\lim_{K\rightarrow\infty}\exp\left\{-K{\cal H}_{Edw}\right\}$
where ${\cal H}_{Edw}=
\sum_{i,z} \left[ \delta_{n_i^z 2}\delta_{n_i^{z-1} 0}\delta_{n_i^{z-2} 0}
+ \delta_{n_i^z 1}\delta_{n_i^{z-1} 0}
\left(1 - \delta_{n_i^{z-2} 2}\right) \right]$.

\bibitem{epl} M. Tarzia,
A. de Candia, A. Fierro, M. Nicodemi and
A. Coniglio,
Europhys. Lett. {\bf 66}, 531 (2004).
\bibitem{jamming} A. Fierro,
M. Nicodemi, M. Tarzia, A. de Candia and A. Coniglio,
{\em cond-mat}/0412120.
\bibitem{MP} M. M\'{e}zard and G. Parisi, Eur. Phys. J. B {\bf 20}, 217 (2001).
G. Biroli and M. M\'{e}zard, Phys. Rev. Lett. {\bf 88}, 025501 (2002).

\bibitem{nota_cryst} Along with the fluid solution, also a crystalline solution
characterized by the breakdown of the horizontal translational invariance,
can be found. Here we consider the case where the crystallization is avoided.
%and limit the calculation to the fluid phase.

\bibitem{nota_func}
${\cal K}_n^{(z)} = A_n^{(z)} / S^{(z)}$ where
$A_n^{(z)} = 1+\delta_{n,1}e^{s_{1}^{(z)}}$ and
$S^{(z)}=A_1^{(z)}+e^{s_{2}^{(z)}}$. Moreover we have that
$P_n^{(z)} = 1+ \delta_{n,1} e^{u_{1}^{(z)}} + \delta_{n,2} e^{v_{1}^{(z)}}+
e^{v_{2}^{(z)}}$,
$R_1^{(z)} = \delta_{n,1} e^{d_{1}^{(z)}} + \delta_{n,2}e^{c_{1}^{(z)}}
+ e^{c_{2}^{(z)}}$,
$Q^{(z)} = 1+e^{d_{1}^{(z-1)}} %\left(
[1+ e^{u_{2}^{(z+1)}}]%\right)
+e^{d_{2}^{(z-1)}} %\left(
[1+e^{u_{1}^{(z+1)}} +  e^{u_{2}^{(z+1)}}]
%\right)
+ R_2^{(z)}$.

\bibitem{nota_free}
The free energy, $F$, is given by
$F=\sum_{z=0}^H \Delta F_s^{(z)}-\frac{(k-1)}{2}\sum_{z=0}^H
\Delta F_{l,2}^{(z)}-\sum_{z=0}^{H-1}\Delta F_{l,1}^{(z)}$,
where,
$e^{- \Delta F_{s}^{(z)}}=
%\left
(S^{(z)}
%\right
)^{k-1}Q^{(z)}+\sum_{n=1,2}
\phi_n^{(z)}
%e^{\beta(\mu_n-m_ngz)}
%\left
( A_n^{(z)}
%\right
)^{k-1}P_n^{(z+1)}R_n^{(z-1)}$,
$e^{-\Delta F_{l,1}^{(z)}}=1+\sum_{n=1,2}e^{c_n^{(z)}}
%\left
(1+\delta_{n,2}
e^{u_1^{(z+1)}}+e^{u_2^{(z+1)}}
%\right
)+\sum_{n=1,2}e^{d_n^{(z)}}P_n^{(z+1)}$,
$e^{-\Delta F_{l,2}^{(z)}}=1+2 e^{s_1^{(z)}}+2 e^{s_2^{(z)}}+e^{2 s_1^{(z)}}$.

%\bibitem{nota} $\Delta h / h$ is
%positive when the small grains is in average above than large grains (RBNE)
%and negative otherwise (BNE).
%The value $\delta_c$, defined by $\Delta h / h(\delta_c)\equiv 0$, does
%not in general coincide with zero, depending on what the local densities of the
%two species are.

%\bibitem{nota1} The susceptibility is defined by 
%$\chi \equiv \langle \left(z_1-z_2  \right)^2 \rangle
%- \langle \left(  z_1-z_2 \right) \rangle^2$. In the mean field approximation 
%this quantity is be computed using the
%fluctuation-dissipation theorem: $\chi \simeq \left(
%\frac{\partial}{\partial (m_1-m_2)} \right)^2 \ln {\cal Z}$.

\bibitem{nota_tap} Each
tap consists in raising the bath temperature to a finite value
and, after a lapse of time quenching the bath temperature back to zero.
The measurements are performed when the shake is off and the system is at rest.

\end{thebibliography}
\end{document}


cond-mat0505724.extracted_bib
\begin{thebibliography}{40}
\bibitem{rev_segr}
J. Bridgwater, Chem. Eng. Sci. {\bf 50}, 4081 (1995).
J.M. Ottino and D.V. Khakhar,
Ann. Rev. Fluid Mech. {\bf 32}, 55 (2000).
T. Shinbrot and F.J. Muzzio, Physics Today {\bf 53}, 25 (2000).
Rosato et al., Chem. Eng. Sci. {\bf 57}, 265 (2002).
\bibitem{nagel}  M.E. Mobius,
X. Cheng, G.S. Karczmar, S.R. Nagel, H.M. Jaeger,
{\em cond-mat}/0405328.
\bibitem{rosato} A. Rosato,
K.J. Strandburg, F. Prinz,  R.H. Swendsen,
Phys. Rev. Lett. {\bf 58}, 1038 (1987).
\bibitem{luding}  D.C. Hong,
P.V. Quinn, S. Luding,
Phys. Rev. Lett. {\bf 86}, 3423 (2001).
\bibitem{breu} A. P. J. Breu,
H.-M. Ensner, C. A. Kruelle, and I. Rehberg,
Phys. Rev. Lett. {\bf 90}, 014302 (2003)

%\bibitem{bridgewater} J. Bridgwater,
%Powder Technol. {\bf 15}, 215 (1976).
%J.C. Williams,
%Powder Technol. {\bf 15}, 245 (1976).
\bibitem{knight}  J.B. Knight,
H.M. Jaeger, S.R. Nagel,
Phys. Rev. Lett. {\bf 70}, 3728 (1993).

\bibitem{shinbrot}  T. Shinbrot and F.J.  Muzzio,
Phys. Rev. Lett. {\bf 81}, 4365 (1998).
\bibitem{jenkins} J.T. Jenkins and D.K. Yoon,
Phys. Rev. Lett. {\bf 88}, 194301 (2002).
\bibitem{herrmann} L. Trujillo,
M. Alam and H. Herrmann,
Europhys. Lett. {\bf 64}, 190 (2003).
\bibitem{mullin}
P.M. Reis and T. Mullin, Phys. Rev. Lett. {\bf 89}, 244301 (2002).
P. M. Reis,
G. Ehrhardt, A. Stephenson and T. Mullin
Europhys. Lett. {\bf 66}, 357 (2004).
\bibitem{both}
J.A. Both and D.C. Hong, Phys. Rev. Lett. {\bf 88}, 124301 (2002).
\bibitem{seg}
M. Tarzia,
A. Fierro, M. Nicodemi, and A. Coniglio,
Phys. Rev. Lett. {\bf 93},
198002 (2004).

\bibitem{Edwards}  S.F. Edwards and R.B.S. Oakeshott,
Physica A {\bf 157}, 1080 (1989).
A. Mehta and S.F. Edwards, Physica A {\bf 157}, 1091 (1989).

\bibitem{e1}  M. Nicodemi, Phys. Rev. Lett. {\bf 82}, 3734 (1999).
A. Barrat,
J. Kurchan, V. Loreto and M. Sellitto,
Phys. Rev. Lett. {\bf 85}, 5034 (2000).
J.J. Brey,
A. Prados, B. S\'{a}nchez-Rey,
Physica A {\bf 275}, 310 (2000).
A. Coniglio and M. Nicodemi, Physica A {\bf 296}, 451 (2001).
D. S. Dean and A. Lef\`{e}vre, Phys. Rev. Lett. {\bf 86}, 5639 (2001).
H. A. Makse and J. Kurchan, Nature {\bf 415}, 614 (2002).
A. Fierro,
M. Nicodemi and A. Coniglio,
Europhys. Lett. {\bf 59}, 642 (2002).
G. De Smedt,
C. Godreche, J.M. Luck,
Eur. Phys. J. B {\bf 32}, 215-225 (2003).
\bibitem{nfc} M. Nicodemi,
A. Fierro, A. Coniglio,
Europhys. Lett. {\bf 60}, 684 (2002). A. Fierro,
M. Nicodemi and A. Coniglio,
Phys. Rev. E {\bf 66}, 061301 (2002). A. Fierro,
M. Nicodemi, A. Coniglio,
Jour. Phys.: Cond. Matt. {\bf 15}, S1095 (2003).

\bibitem{Evans} M. Dijkstra,
R. van Roij and R. Evans,
Phys. Rev. E {\bf 59}, 5744 (1999). M. Fasolo and P. Sollich, Phys. Rev. Lett.
{\bf 91}, 068301 (2003).

\bibitem{depletion} M. Bose,
P.R. Nott and V. Kumaran,
Europhys. Lett. {\bf 68},
508 (2004). An effective attraction
between two large grains
was also found in a $2d$ vertically vibrated granular bed of smaller ones
in Ref.  D.A. Sanders,
M.R. Swift, R.M. Bowley, and P.J. King,
Phys. Rev. Lett. {\bf 93}, 208002 (2004).

\bibitem{nota_proj} For a given grain configuration, $\{n_i^z\}$,
$\Pi (\{n_i^z\})$ has a tractable expression:
$\Pi (\{n_i^z\}) =\lim_{K\rightarrow\infty}\exp\left\{-K{\cal H}_{Edw}\right\}$
where ${\cal H}_{Edw}=
\sum_{i,z} \left[ \delta_{n_i^z 2}\delta_{n_i^{z-1} 0}\delta_{n_i^{z-2} 0}
+ \delta_{n_i^z 1}\delta_{n_i^{z-1} 0}
\left(1 - \delta_{n_i^{z-2} 2}\right) \right]$.

\bibitem{epl} M. Tarzia,
A. de Candia, A. Fierro, M. Nicodemi and
A. Coniglio,
Europhys. Lett. {\bf 66}, 531 (2004).
\bibitem{jamming} A. Fierro,
M. Nicodemi, M. Tarzia, A. de Candia and A. Coniglio,
{\em cond-mat}/0412120.
\bibitem{MP} M. M\'{e}zard and G. Parisi, Eur. Phys. J. B {\bf 20}, 217 (2001).
G. Biroli and M. M\'{e}zard, Phys. Rev. Lett. {\bf 88}, 025501 (2002).

\bibitem{nota_cryst} Along with the fluid solution, also a crystalline solution
characterized by the breakdown of the horizontal translational invariance,
can be found. Here we consider the case where the crystallization is avoided.
%and limit the calculation to the fluid phase.

\bibitem{nota_func}
${\cal K}_n^{(z)} = A_n^{(z)} / S^{(z)}$ where
$A_n^{(z)} = 1+\delta_{n,1}e^{s_{1}^{(z)}}$ and
$S^{(z)}=A_1^{(z)}+e^{s_{2}^{(z)}}$. Moreover we have that
$P_n^{(z)} = 1+ \delta_{n,1} e^{u_{1}^{(z)}} + \delta_{n,2} e^{v_{1}^{(z)}}+
e^{v_{2}^{(z)}}$,
$R_1^{(z)} = \delta_{n,1} e^{d_{1}^{(z)}} + \delta_{n,2}e^{c_{1}^{(z)}}
+ e^{c_{2}^{(z)}}$,
$Q^{(z)} = 1+e^{d_{1}^{(z-1)}} %\left(
[1+ e^{u_{2}^{(z+1)}}]%\right)
+e^{d_{2}^{(z-1)}} %\left(
[1+e^{u_{1}^{(z+1)}} +  e^{u_{2}^{(z+1)}}]
%\right)
+ R_2^{(z)}$.

\bibitem{nota_free}
The free energy, $F$, is given by
$F=\sum_{z=0}^H \Delta F_s^{(z)}-\frac{(k-1)}{2}\sum_{z=0}^H
\Delta F_{l,2}^{(z)}-\sum_{z=0}^{H-1}\Delta F_{l,1}^{(z)}$,
where,
$e^{- \Delta F_{s}^{(z)}}=
%\left
(S^{(z)}
%\right
)^{k-1}Q^{(z)}+\sum_{n=1,2}
\phi_n^{(z)}
%e^{\beta(\mu_n-m_ngz)}
%\left
( A_n^{(z)}
%\right
)^{k-1}P_n^{(z+1)}R_n^{(z-1)}$,
$e^{-\Delta F_{l,1}^{(z)}}=1+\sum_{n=1,2}e^{c_n^{(z)}}
%\left
(1+\delta_{n,2}
e^{u_1^{(z+1)}}+e^{u_2^{(z+1)}}
%\right
)+\sum_{n=1,2}e^{d_n^{(z)}}P_n^{(z+1)}$,
$e^{-\Delta F_{l,2}^{(z)}}=1+2 e^{s_1^{(z)}}+2 e^{s_2^{(z)}}+e^{2 s_1^{(z)}}$.

%\bibitem{nota} $\Delta h / h$ is
%positive when the small grains is in average above than large grains (RBNE)
%and negative otherwise (BNE).
%The value $\delta_c$, defined by $\Delta h / h(\delta_c)\equiv 0$, does
%not in general coincide with zero, depending on what the local densities of the
%two species are.

%\bibitem{nota1} The susceptibility is defined by 
%$\chi \equiv \langle \left(z_1-z_2  \right)^2 \rangle
%- \langle \left(  z_1-z_2 \right) \rangle^2$. In the mean field approximation 
%this quantity is be computed using the
%fluctuation-dissipation theorem: $\chi \simeq \left(
%\frac{\partial}{\partial (m_1-m_2)} \right)^2 \ln {\cal Z}$.

\bibitem{nota_tap} Each
tap consists in raising the bath temperature to a finite value
and, after a lapse of time quenching the bath temperature back to zero.
The measurements are performed when the shake is off and the system is at rest.

\end{thebibliography}