\section{Introduction} 
Granular materials are constituted by macroscopic 
rigid objects that are too big to undergo Brownian 
motion, and that dissipate energy when they collide 
\cite{duran_2000,jaeger_RMP_1996}. Therefore, a granular 
assembly in the presence of gravity tends rapidly 
to a static configurations, where all the material 
simply forms a pile on the bottom of whichever container 
holds it. However, using a continuous injection of 
energy, it is possible to maintain a granulate in a 
non-static configuration. The energy is usually injected 
through vertically shaking the system, using a sinusoidal 
displacement. In such a case, the motion's amplitude times 
the square of the angular frequency gives the maximal 
acceleration of the exciting system, and this value 
divided by the gravitational acceleration is called 
its non-dimensional acceleration $\Gamma$. Due to the 
permanent flow of energy in this kind of granular aggregate,
one may observe the appearance 
of self-organization, allowing for some very 
interesting phenomena such as oscillons, and a wide 
variety of beautiful patterns \cite{bizon_PRL_1998}. 
But it is also possible, with
this permanent provision of energy, to obtain 
statistically stationary configurations containing 
just a few homogeneous phases. Under these conditions, 
granular materials may resemble gases, liquids 
\cite{roeller_PRL_2011} or crystals 
in equilibrium; and in fact these systems have been 
used as models for atomic, molecular, and colloidal 
ensembles in thermal equilibrium 
\cite{jaeger_RMP_1996,tatal_PRL_2000}.  

Thus, condensed phases resembling 
liquids or solids, both crystalline and glassy, have 
been reported in granular materials under vertical 
shaking excitation. Their stationary structures have 
been described \cite{knight_PRE_1998,eshuis_Phys_Flu_2007}, 
and the jammed-like system relaxation dynamics that 
they display has been reported \cite{reis_PRL_2007}. 
In particular, two-dimensional crystalline configurations with 
hexagonal order \cite{olafsen_PRL_1998,olafsen_PRE_1999} 
(and, with the appropriate confinement, 
square order \cite{melby_JPCM_05,vega_reyes_2008}), have been found.

These crystalline clusters have been 
obtained in the laboratory 
either by increasing the packing fraction, that is, 
the total volume (area in 2D systems) occupied by 
the particles, divided by the 
volume (area) of the container; or by reducing the 
amplitude of the driving. 
In the first case the process is basically equal to 
the crystallization of hard disks in equilibrium 
\cite{olafsen_PRL_2005,reis_PRL_2006}: as density is 
increased the aggregate displays a transition from 
a fluid to an hexagonal crystalline phase, presumably 
going through a narrow hexatic phase \cite{olafsen_PRL_2005,reis_PRL_2007}, 
a behavior similar in all ways to the well-known 
entropicaly-driven crystallization of monodisperse 
hard disks \cite{alder_PR_1962,engel_PRE_2013}.
But surprisingly, the formation of hexagonal crystals 
has been also observed, both experimentally and numerically,
for values of $\Gamma \approx 1$, 
\cite{olafsen_PRL_1998,nie_EPL_2000,perez-angel_PRE_2011}, 
at very low packing fractions ($\phi = 0.42$ in the 
experiment vs.\ $\phi \approx 0.71$ for equilibrium 
crystallization in a gas of hard disks).
This crystallization gives an extreme example of the anomalous 
behavior of vibrated dissipative monolayers, in both 
their morphology and dynamics, but the differences with  
a two-dimensional equilibrium gas are not exclusive to 
condensate phases. One particular issue of interest in 
these systems is the evolution of density fluctuations, 
including the creation of large dense clusters
\cite{perera-burgos_PRE_2010}, 
a situation that resembles the density fluctuations in a 
gas approaching a gas-liquid transition. 

The existence of effective 
potentials related to the physical parameters of the
shaken granulate, especially to its restitution coefficient, 
has been proposed in \cite{bordallo-favela_EPJE_2009},
in the 
spirit of this analogy to the gas-liquid transition, and 
also to the gas-crystal transition mentioned before. 
It has been shown in that work how the structure of 
a given granulate, described by its Pair Distribution Function 
(PDF) $g(r)$ can be related to an attractive potential, 
and how the naive expectation of more dissipation implying 
a more attractive potential does get realized. 
This idea of introducing effective potentials in a
granular medium ---dissipative and out of equilibrium 
as it is---, is an extension of the previous generalizations of 
interparticle potentials: for instance, the well known
Dejaurmin, Landau, Vervey and Overbeek (DLVO) potential
\cite{derjaguin_APCURSS_1941,verwey_1948,hiemenz_1997} 
for colloids with electrostatic interactions, and the
depletion interactions due to entropic effects,
modeled originally by Asakura and Oosawa 
\cite{asakura_JCP_1954,asakura_JPS_1958}.

In extending these ideas to vibrated granular matter
one simply wants to know 
whether it is possible to represent 
the structure of a stationary state as if it were 
an equilibrium state of a well defined thermodynamic 
system. 
Besides the monodisperse quasi-2D granular gas mentioned before 
\cite{bordallo-favela_EPJE_2009},
that we continue exploring here, 
this approach has also been tried in the case of
a two-species 
granular gas with horizontal shaking, whose dynamics has been 
interpreted in terms of a depletion interaction 
system \cite{ciamarra_PRL_2006}.
The interaction of particles trough a potential well 
implies a conservative system in which the energy associated 
to the relative position of two particles is always balanced 
by their total kinetic energy, producing a characteristic 
structure, that is, characteristic spatial 
distribution and correlations.
What we propose here is that a set of particles interacting trough 
dissipative collisions which uniformly receives energy from 
an oscillating support could, in general, behave effectively 
as a conservative system, in the sense of statistically 
reproducing the characteristic spatial distribution of 
particles produced by the last. 


In this paper we present a computational study 
(molecular dynamics and Monte Carlo simulations)  
of a quasi-2D gas of spheres that collide inelastically,
excited via vertical shaking. We intend to show that 
the homogeneous steady state of 
this system may be described by means of an analogy 
with equilibrium liquid theory. This is done by means 
of performing molecular dynamics simulations, getting 
the corresponding pair distribution functions, and 
obtaining from them  ---using the Ornstein-Zernike 
equation with the Percus-Yevick closure--- the associated 
effective potentials. Furthermore, in order to test 
the correctness of this analogy, we performed Monte Carlo 
simulations of hard discs interacting through the 
effective potentials obtained by the aforementioned 
molecular dynamic simulations and the simple 
liquid theory, and compare their respective pair 
correlation functions as a check of the correctness of the
approach. 


\section{Simulation scheme} 

We simulate a system quite similar to the usual 
experimental set-ups that have been used in many 
studies of granular gases \cite{olafsen_PRL_1998,olafsen_PRE_1999}. 
A monodisperse granulate composed of $N = 2000$ hard 
spheres of diameter $\sigma = 0.5$ cm is confined by
two horizontal flat interfaces (``floor'' and ``ceiling'') 
perpendicular to the direction of gravity, 
separated by a distance $h = 3$, $\sigma = 1.5$ cm, making sure in all cases 
that in the stationary state the grains never reach the ceiling. 
As an exception, for some control simulations conducted without 
dissipation, the vertical spacing was set at $h = 1.6$, $\sigma = 0.8$ cm.
The two-dimensional packing fraction is defined by the total projected 
area covered by particles divided by the simulation cell area, that is,
$\phi = N \pi \sigma^2/ (4 L^2)$. The horizontal cell is a square 
whose side is adjusted depending on the desired packing fraction,
and where periodic boundary conditions are used in the 
horizontal directions. A sketch of the system is presented in 
top and lateral views in Fig.~(\ref{esquema_sim}).

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=8.0cm,angle=0]{esquema_sim_22a.eps} \\
a) Top view. \\ 
 \\
\includegraphics[width=8.0cm,angle=0]{esquema_sim_22b.eps} \\
b) Lateral view. \\ 
\end{tabular}
\end{center}

\caption{Scheme of the simulated system in a) Top view, and b) Lateral view. 
2000 spheres are vertically shaken by a vertical sinusoidal motion of 
the containing cell.} 
\label{esquema_sim} 
\end{figure}

The assembly is sinusoidally vibrated in the vertical direction 
(parallel to gravity), 
defining the position of the two horizontal plates as 
$z_b = A \sin(2 \pi f t)$ for the floor and $z_t = z_b + h$
for the ceiling;
here $A$ is the amplitude and $f$ is the frequency
of the shaking. The simulations are conducted with 
7 different values of the restitution coefficient,
given by $\epsilon = $ 0.40, 0.50, 0.60, 0.70, 0.80, 0.90 and 1.0, 
this last one used for the non-dissipative test.
We also use four different values of the packing fraction, given by
$\phi = $ 0.20, 0.30, 0.40 and 0.50, and at least 5 values of the 
non-dimensional vertical acceleration $\Gamma$, defined,
as mentioned before, 
by $\Gamma = (2 \pi f)^2 A / g$, and given by 
$\Gamma = $ 1.05, 1.15, 1.25, 1.35 and 1.45 for each case. 
For the study of resonances, to be discussed 
latter, we also used a much finer sweep over $\Gamma$.
The gravitational acceleration is
$g = 9.81$ m/s$^2$, the frequency is kept fixed at 
$f = 70$ Hz, and the tangential friction coefficient is 
fixed at $\mu = $ 0.25 (again, with the exception of the control 
simulation that has $\mu = 0$). Notice that it is 
possible to scale out some magnitudes, but we prefer 
to use real physical quantities; it should be understood 
that the reported results can be scaled to 
other sets of parameters, making sure that the process
is such that $g$ keeps its numerical value.
The bounds in the ranges of parameters explored are set
by the need to keep some of the basic characteristics
of the system; in particular, larger values of $\Gamma$
allow sometimes a grain to jump over another, and larger
values of $\phi$ begin to introduce crystallization.

The simulations are conducted using an Event-Driven 
Molecular Dynamics (EDMD) algorithm 
\cite{poschel_2005,rapaport_2004}. 
In this algorithm, one calculates analytically 
the trajectory in space of all spheres, and finds 
the times for all possible future collisions. 
These collision times are sorted, and the 
system is then ballistically advanced to the
time of the first collision. This collision is then
realized, preserving
linear momentum and angular momentum around the 
point of contact, but not energy. 
Immediately after, the list of future collisions
is updated, then sorted, and the system is again advanced 
ballistically to the next collision. Dissipation arises
due to two effects: as seen in the center of mass, 
the particles separate with normal velocities which are 
only a fraction $\epsilon$ of the incident normal 
velocities, making the collision viscoelastic; 
also, the tangential velocities are reduced by
a tangential restitution coefficient $\epsilon_s$,
which depends on the friction coefficient $\mu$,
keeping the conservative limit $\epsilon_s = 1$ for
$\mu = 0$.
Both dissipative processes have
also an effect on the rotations of the particles.
The full equations used to solve the collision are 
given in \cite{perez_2008_Pramana, schwager_2008_EPJE}.
A couple of extra points should be mentioned: 
first: collision times between spheres can be calculated 
exactly even in the presence of gravity, in 
the same way that they can 
be solved without gravity, just solving a 
quadratic equation.
Second: solving for the time of collision 
between a particle falling in a parabola and the
floor (or the ceiling) moving sinusoidally requires a 
numerical scheme. Here we have implemented
root-finding routine specifically designed to calculate 
this type of intersection.

It is understood that we always use $\Gamma > 1$, to avoid 
situations where a grain no longer does finite 
length bounces against 
either other grain or a wall. Such a situation  
forces an EDMD code to perform an infinite number of 
collisions in a finite time, that is, to undergo inelastic 
collapse \cite{mcnamara_PRE_1996}.
Some remedies can be applied, such as to take the restitution 
coefficient $\epsilon$ equal to 1 for very small relative 
velocities \cite{mcnamara_PRE_1998}, or to ignore collisions that are less than some 
small $\Delta t$ apart, and even freezing the dynamics
of very dense clusters \cite{gonzalez_EPJST_2009}.
But any way one looks at it, ED is 
unable to deal with non-bouncing grains.
This issue is avoided for $\Gamma > 1$.

\section{Results: Pair Distribution Functions and Effective Potentials}
For our purposes, the main results from the simulations
mentioned above are their PDFs, which are extracted from
a sampling of configurations using well known methods. These
PDFs are then assumed to originate in an equilibrium 2D fluid,
and to be determined, therefore, by an interparticle potential, 
along with temperature and density.
However, the relationship between these potentials and the corresponding PDFs
is not easily resolved. We will use here the
homogeneous Ornstein-Zernike (OZ) equation
$$
h(r_{12}) = c(r_{12}) + \rho \int 
\mathrm{d}^3{\bf r}_3 \, c(r_{13}) \, h(r_{32})
$$
where $r_{ij} = |{\bf r}_i - {\bf r}_j|$, $\rho$ is the number particle density,
and $h(r) \equiv g(r) - 1$ and $c(r)$ are the \emph{total correlation function\/} and
the \emph{direct correlation function\/}, respectively. 
In practice, $c(r)$ is actually 
defined by the OZ equation itself \cite{fisher_JMP_1964}. 
The relation between these two correlation functions can be
written in a much simpler form using Fourier transform, where one finds
\begin{equation}
\hat h(k) = \frac{\hat c(k)}{1 - \rho\, \hat c(k)}.
\label{hdek}
\end{equation}
Still, in the absence of an independent relationship, this equation
is of little use.
For hard spheres ---and, in general, short 
range potentials, in either dense or diluted si\-tua\-tions---, a good
closure relationship is given by the 
Percus-Yevick (PY) approximation \cite{textbooks}
\begin{equation}
c(r) = [h(r) + 1][1 - \exp(u(r)/k_B T) ].
\label{PYapp}
\end{equation}
With this approximation one can then generate a value for the quotient
$u(r)/k_B T$, once the $g(r)$ and $\rho$ are given. This requires: the evaluation 
of $\hat h(k)$, which in 2D involves the zeroth order Hankel transform of $g(r) - 1$,
solving (\ref{hdek}) for $\hat c(k)$, doing the inverse Hankel transform to get
$c(r)$, and finally solving (\ref{PYapp}) to get $u/k_B T$. 

\subsection{Effective Potentials for the Conservative Case}

Let us start showing the behavior of the 
calculated Effective Potentials for the 
conservative case, where we have set 
$\epsilon = 1$ and $\mu = 0$.  This is a 
situation that clearly falls outside of 
experimental possibilities, and we are 
reporting it here only as a validation 
of the method, that is, of the use of the 
OZ-PY scheme to extract effective potentials 
from the simulated PDFs. In 
Figs.~(\ref{gder-conservative}) 
and~(\ref{bu-conservative}) we have an 
example of the PDFs and effective potentials 
for the non-dissipative granulate, with 
the parameters given in the figure caption. 
It is clear that the OZ-PY procedure 
manages to extract potentials that are
very close to those of a hard disk fluid, 
even if we have to take into account that 
our system is not really perfectly 2D, 
something that is noticeable in: (a) the 
behavior of the PDFs for $r$ slightly smaller 
than $\sigma$, where a non-zero  probability 
is found, (b) the mild hump between $\sigma$ 
and $2 \sigma$ found in the potential,
and (c) a very narrow potential well that 
appears at contact. These variations with 
respect to the perfect hard disk potential 
are related to the apparent penetrability 
of spheres, that is, to the fact that the 
2D projection of this 3D system shows 
distances between centers smaller than one 
diameter, obviously without any real 
interpenetration involved. This apparent 
overlapping could lead to the very narrow 
potential at contact, while the mild hump 
may be an artifact produced by confinement 
of the spheres in a narrow quasi-2D space.
Still, it should be noticed that the numerical 
values for the variations of the potential are 
quite small, even more if compared to the 
non-conservative examples to be covered next.
It should also be noticed that the potential 
shows some short wavelength oscillations, 
an unavoidable by product of the truncation in 
the back and forth Fourier transforms required 
for the solution. 

\begin{figure}
\includegraphics[width=8.0cm,angle=270]{gder-conservativo.ps} 
\caption{Pair Distribution Function for the 
conservative case, $\epsilon = 1$, $\mu = 0$. 
Here we show the $\phi = 0.2$ case, for all 
the $\Gamma$ values considered ($\Gamma =$ 
1.05, 1.15, 1.25, 1.35 and 1.45). The different 
lines superpose almost perfectly, implying that
the vibration strength is irrelevant in the 
confined conservative case. PDFs for larger 
values of $\phi$ behave in a completely similar 
manner, just showing some more structure close 
to contact.} 
\label{gder-conservative} 
\end{figure} 

\begin{figure}
\includegraphics[width=8.0cm,angle=270]{bu-conservativo.ps} 
\caption{Effective potentials corresponding 
to the PDFs given in the previous figure. 
The small oscillations induced by truncation 
in the Fourier transforms appear clearly, 
together with a narrow and not too deep well 
at contact, plus a mild potential barrier 
between $r/\sigma = 1$ and $r/\sigma = 2$. 
As mentioned in the text, this narrow well 
and the mild potential barrier are most 
probably artifacts due to the imperfect 
two-dimensional character of the granulate.}
\label{bu-conservative} 
\end{figure} 

\subsection{Effective Potentials as a Function 
of Energy Dissipation, and Forcing Resonances}

Going now into the dissipative case, let us begin
by showing some examples of the behavior of the 
effective potentials found for different sets of 
parameters. In Fig.~(\ref{graf_extra}) we show the 
effective potentials obtained for a constant packing 
fraction $\phi = 0.40$ for four different 
non-dimensional accelerations $\Gamma =$ 1.05, 1.15, 
1.35 and 1.45. 
(Some of the PDFs from where these potentials were 
extracted are shown latter, in Fig.~(\ref{varias})).
As expected, and in coincidence with the experimental 
results reported in \cite{bordallo-favela_EPJE_2009},
there is an attractive well in the effective potential, 
and a tendency to get more attractive potentials as
the restitution coefficient is lowered. This effect can 
be understood as a result of the smaller distances 
between particles that are found after mutual inelastic 
collisions, as compared with those distances expected 
for elastic ones, something that occurs simply because 
the smaller separation velocities found after 
the contacts. In this sense we find that the PDFs for
steady states of dissipative systems register, in general, 
an increase of the spatial correlation at contact for 
decreasing values of the restitution coefficient, and 
thus, a tendency of the effective potentials to display 
deeper wells for more dissipative dynamics. 

There is, however, one general exceptions to this rule,
within the range of parameters we have explored. An
example is given in Fig.~(\ref{graf_extra}-b), where 
the deepest well actually corresponds to 
$\epsilon = 0.80$, which is not the most dissipative 
value studied, and other example is given in 
Fig.~(\ref{graf_extra}-c), where the deepest well 
corresponds now to the most elastic case $\epsilon = 0.90$. 
Similar anomalies appear for other combinations of 
parameters, although "normal" behavior always 
happens for very low or very high $\Gamma$ values.  

\begin{figure}
\includegraphics[width=8.0cm,angle=270]{graf_extra_2.ps} 
\caption{Effective potentials corresponding to the
PDFs obtained for 2D packing fraction $\phi = 0.40$, at 
four different values of the restitution coefficient, 
at five different values of forcing $\Gamma$.
(a): $\Gamma = 1.05$; 
(b): $\Gamma = 1.15$; 
(c): $\Gamma = 1.25$;
(d): $\Gamma = 1.35$;
(e): $\Gamma = 1.45$.} 
\label{graf_extra} 
\end{figure} 

It is clear then that there is no simple scaling
of the effective potentials as $\epsilon$ changes:
a larger restitution coefficient may be associated 
to a deeper potential well, or the potentials for 
different values of $\epsilon$ may cross each other, 
although keeping the general trend of more attractive 
potentials for smaller $\epsilon$. To understand 
this peculiarities, let us concentrate in the effect 
that different values of the forcing $\Gamma$ 
have in the dynamics: we find that the clustering 
which is mainly responsible for the apparent 
increasing attraction between particles can be 
enhanced for a narrow range of this parameter. 
This is due to resonances that appear at narrow 
windows in forcing, whose 
location depend on the restitution coefficient 
but are independent of the packing fraction. 
The resonance arises from a simplification in 
the motion of the granulate, in the following 
sense: ignoring for the moment being the 
particle-particle collisions, the motion of these 
particles is just a collection of parabolic flights 
above the floor. In terms of vertical displacement, 
these parabolas may be of just 
one or a few types (modes), that is, of just one 
or a few heights, and this is what constitutes a
resonance. On the contrary, the granulate may 
display a form of motion composed of many different 
modes (parable with different heights).
Due to the dissipative-forced character of the motion, 
a mode corresponds not just to the presence of 
one or a few parabolic parabolic heights, 
but also involves synchronization, and so in the
granulate we will have many particles displaying
the same vertical motion. Parabolas of different
heights will of course be unsynchronized. 

The evidence for these two forms of motion ---simple
an complex, synchronous and asynchronous---, now 
incorporating the effects of collisions, is given 
in Fig.~(\ref{prob_z}), which shows the probability 
density for the height of the particles above the 
floor, $p(z)$, vs.\ height $z$. In this figure we can 
see now a single and well defined peak that appears 
for $\Gamma = 1.14$, decreases in height and 
finally splits as $\Gamma$ grows, going in this 
way from a mode with a one-height parabola to a 
mode with two different parabolas. This unity in 
the height of the vertical motion is relevant to the 
horizontal behavior of the system because, among 
particle-particle collisions, only those that 
happen between particles with different heights 
introduce horizontal scattering, and therefore 
reduce plane correlations. On the contrary, for 
parabolic motions of the same height, collisions 
are perfectly lateral, and always reduce the in-plane
kinetic energy. 

\begin{figure}
\includegraphics[width=8.0cm,angle=270]{z_rs_e070_p040.ps} 
\caption{Probability density for the height of the 
particles above the floor for 2D packing fraction 
$\epsilon = 0.70$ and $\phi = 0.40$, at several different 
values of $\Gamma$.
$\Gamma = 1.100$;  
$\Gamma = 1.120$; 
$\Gamma = 1.140$;
$\Gamma = 1.150$;
$\Gamma = 1.165$;  
$\Gamma = 1.168$; 
$\Gamma = 1.170$;
$\Gamma = 1.180$;
$\Gamma = 1.200$;  
$\Gamma = 1.220$; 
$\Gamma = 1.250$.} 
\label{prob_z} 
\end{figure} 

Thus, it is necessary to consider not only 
the absolute value of the exciting non-dimensional 
acceleration $\Gamma$, but also whether or not there 
are coincidences (resonances) in the vertical motion 
of the particles. These resonant phenomena have 
already been seen in the case of oscillons 
\cite{hill_AMM_00}. In Fig.~(\ref{snapshots_p030}) 
we show the snapshots of the simulated gas for 
$\Gamma = 1.15$, $\phi = 0.40$ and four different 
restitution coefficients $\epsilon$, corresponding 
to the Fig.~(\ref{graf_extra}-c) in which anomalous 
behavior is detected, in particular for the value 
$\epsilon = 0.80$.

Another way of showing the effects of these resonances 
is given in Fig.~(\ref{gMaxvsGamma1}), where 
correlations at contact, defined by the maximum value 
of the PDF, are plotted as a function of $\Gamma$,
for a fine tuned exploration of this parameter for 
$\epsilon = 0.80$ and $\phi = $ 0.20, 0.30 and 0.40. 
There one can see that the overall response of the 
maximum value of $g(r)$ is to decrease with growing 
$\Gamma$ (growing shaking amplitude), since large 
shaking intensity or ``granular temperature'' tends 
to destroy correlations. If the maxima of the PDFs 
are plotted as a 
function of $\Gamma$, avoiding resonant values, they 
show a roughly exponential decay with $\Gamma$.
Seeing that the correlation becomes very large for
$\Gamma \to 1$, we try a fit
$g_{Max}(r) = g_0 + A \exp \left((\Gamma - 1)/(\Gamma_0 - 1) \right)$
For the $\phi = 0.4$, $\epsilon = 0.8$ case we get
$g_0 = 3.2$, $A = 8.8$ and $\Gamma_0 = 1.041$. 
This fit can be seen in the inset 
of Fig.~(\ref{gMaxvsGamma1}). But, as already
described, this monotonic behavior is interrupted 
by a narrow peak centered at $\Gamma = 1.22$ 
showing a resonant ``cooling down'' of the 
horizontal degrees of freedom of the system. 
In Fig.~(\ref{bumin_e080}) effective
potential depths for $\epsilon = 0.80$ 
and $\phi = $ 0.20, 0.30 and 0.40, are plotted. 
The monotonic decrement of the maximum interaction 
energy is reversed by a resonant peaking 
attraction in all cases.

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=6.0cm,angle=270]{e060_g115_p040.ps} & \includegraphics[width=6.0cm,angle=270]{e070_g115_p040.ps} \\
a) $\epsilon = 0.60$, $\Gamma = 1.15$ and $\phi = 0.40$ & b) $\epsilon = 0.70$, $\Gamma = 1.15$ and $\phi = 0.40$ \\ 
\includegraphics[width=6.0cm,angle=270]{e080_g115_p040.ps} & \includegraphics[width=6.0cm,angle=270]{e090_g115_p040.ps} \\
c) $\epsilon = 0.80$, $\Gamma = 1.15$ and $\phi = 0.40$ & d) $\epsilon = 0.90$, $\Gamma = 1.15$ and $\phi = 0.40$ \\ 
\end{tabular}
\end{center}
\caption{Snapshots of the simulated granular gas 
corresponding to the same conditions used in the 
previous Figure, showing the clustering tendency 
of the gas as the restitution coefficient decreases, 
cooling down the steady state of the granular gas.}
\label{snapshots_p030}
\end{figure}

\begin{figure}
\includegraphics[width=8.0cm,angle=270]{gmax_fit_e080_rs.ps} 
\caption{Maxima of the PDFs as a function of 
$\Gamma$ for three different packing fractions 
$\phi = 0.20$ (black squares), $\phi = 0.30$ 
(red circles), and $\phi = 0.40$ (green triangles).
A clear resonant peak occurs at $\Gamma$ = 1.22. 
In the inset the maximum of $g(r)$ is shown together 
with the exponential fit given in the text, omitting 
resonant values.}
\label{gMaxvsGamma1} 
\end{figure} 

\begin{figure}
\includegraphics[width=8.0cm,angle=270]{bumin_e080_rs.ps} 
\caption{Minima of pair interaction potentials as 
a function of $\Gamma$ for three different packing 
fractions $\phi = 0.20$ (black squares), 
$\phi = 0.30$ (red circles), and $\phi = 0.40$ 
(green triangles); a clear resonant 
minimum occurs at $\Gamma = 1.22$.}
\label{bumin_e080} 
\end{figure} 

In Fig.~(\ref{var_epsilon}), correlations 
at contact as a function of $\Gamma$ are 
shown for different restitution coefficients at
a fixed packing fraction of $\phi = 0.40$. 
The best defined resonant peak is the one 
shown in Fig.~(\ref{gMaxvsGamma1}) for 
$\epsilon = 0.80$, while for different 
restitution coefficients broader peaks 
appear centered at lower $\Gamma$ values. 
Here, as in the case of an inertial-elastic 
systems, when dissipation increases the 
resonant peak broadens. 

Therefore, in general, for $\Gamma$ values 
out of resonance the decreasing tendency 
of correlation at contact with increasing 
shaking strength dominates, together with 
an increase in correlation with decreasing 
elasticity, as can be seen in 
Fig.~(\ref{gMaxvsGamma1}). 

In the supplementary material section a 
set of tables containing all the explored 
experimental conditions are presented, 
showing for all cases agreement with the 
general tendency described in the particular 
cases analyzed previously in the text,
and in Figs.~(\ref{gMaxvsGamma1}) 
and~(\ref{var_epsilon}).
 
\begin{figure}
\includegraphics[width=8.0cm,angle=270]{var_epsilon_2.ps} 
\caption{Effective potentials corresponding 
to the PDFs obtained for 2D packing fraction 
$\phi = 0.40$, at four different values of 
the restitution coefficient, at five different 
values of forcing $\Gamma$.
(a): $\Gamma = 1.05$; 
(b): $\Gamma = 1.15$; 
(c): $\Gamma = 1.25$;
(d): $\Gamma = 1.35$;
(e): $\Gamma = 1.45$.} 
\label{var_epsilon} 
\end{figure} 

There are some other noteworthy features 
in the potentials found. In particular, 
they display a kink ---for small $\epsilon$, 
a small well--- for $r/\sigma = 2$, 
signaling the fact that the arrangements of 
three particles in line are specially efficient 
for energy dissipation; it is important to 
notice that no particular features are found 
for the next-nearest distance, 
$r/\sigma = \sqrt{3}$, implying that a 
closed-pack triangle of grains does not 
represent a strongly dissipative 
configuration. For small $\epsilon$, some 
of these potentials also show a wide 
well between $r/\sigma = 1$ and 
$r/\sigma = 2$, at a distance clearly 
larger than contact. These peculiarities
are also found in the original experiments, 
as can be seen for the line denoted 
``painted steel'' in Fig.~(3) of 
ref.~\cite{bordallo-favela_EPJE_2009}, which 
shows a wide potential well between 
$r/\sigma = 1$ and $r/\sigma = 2$;
and in the line denoted ``plastic'' 
in the same figure, where a small kink at 
$r/\sigma = 2$ can be seen.
It is however clear that, given that 
there are some differences between the 
layouts for the experiment and the simulation, 
that there is an experimental dispersion in sizes,
and also some other sources of experimental noise, 
one should not expect that the very detailed 
features found in the strongly controlled 
simulational environment should 
always be present in the experimental results.
Notice finally that the secondary attractive 
well developed (at $r/\sigma = 2$) 
is consistent with the original finding of 
the inelastic collapse singularity that 
apparently halted the simulations 
of McNamara and Young \cite{mcnamara_PRE_1996}.

\section{Validation of the Effective Potentials: Monte Carlo 
Simulations}

The basic assumption of the results 
presented up to now is that the 
stationary configurations of the 
forced-dissipative system studied 
can be understood in terms of an 
equilibrium 2D fluid with an effective
potential that emerges, trough the 
dynamics, from the forcing, density 
and dissipative characteristics of 
the system. To validate this idea, 
we need to see if in fact the
calculated potentials induce the 
same statics in the system. To 
do this we implement in this 
section Monte Carlo 
Simulations (MCS) of a 
corresponding 2D fluid, using a system 
containing the same number of disks 
as their corresponding EDMD experiments, 
affected by the effective potential 
obtained by EDMD; that is, the 
attractive interaction potential between 
pairs of particles as a function of interparticle 
distance is inserted numerically into a 
Monte Carlo simulation. 
It should be noted that a typical attractive effective interaction 
obtained by EDMD has a maximum range between 3 and 5 
particle diameters, and beyond that it can be considered null 
for any practical purpose; this allows us
to introduce a cut-off distance of $ 5\sigma$.
Initial conditions was set by randomly distributing 
a set of disks in the sites of an hexagonal matrix.
The side of the simulation 
cell is chosen to satisfy the desired packing fraction 
for a given number of particles (800 in all runs), 
using $L = \sqrt{N \pi r^2 / \phi}$. 
We used in most cases 20,000 Monte Carlos sweeps.
The PDFs for the granular fluid of disks 
were calculated, and, as can be seen in Fig.~(\ref{mc_e080_g145_p030}), 
they are quantitatively very similar, almost identical, to those
obtained from the original molecular dynamics. 
This accounts for the two equivalent descriptions: on the one hand a fluid
of disks in thermodynamic equilibrium interacting via a potential 
(representing a conservative force by definition); and on the other 
hand a granular gas in stationary state in which the energy 
lost due to inelastic collisions is injected back into the system 
trough sinusoidal vertical shaking.
In Fig.~(\ref{varias}) PDFs measured 
from the 3D-EDMS and from the effective-potential 
fluids in 2D Monte Carlo simulations are compared for different packing 
fractions, shaking accelerations, and restitution coefficients. 
They show excellent agreement in reproducing the pair distribution for both
diluted and concentrated fluids (Fig.~(\ref{varias}-a) and Fig.~(\ref{varias}-b)), 
from the very elastic to the very inelastic limits 
(Fig.~(\ref{varias}-c) and Fig.~(\ref{varias}-d)), 
and for cold or hot gases (high or low $\Gamma$ respectively 
in Fig.~(\ref{varias}-e) and Fig.~(\ref{varias}-f)).  
The very good agreement observed between the PDFs
measured in the fully 3D EDMD of a vertically shaken granular gas in 
a stationary out-of-equilibrium state, and those obtained 
using 2D Monte Carlo simulations of disks interacting trough the 
effective potential obtained in this way, in our opinion validates 
considering a granular gas excited trough shaking not as an ideal 
hard sphere gas but instead as a short-range interacting gas. 
Of course, this interacting granular gas can undergo condensation if 
the temperature (associated to the system energy injection rate) 
is lower than certain threshold value. 


\begin{figure}
\includegraphics[width=8.0cm,angle=270]{mc_e080_g145_p030.ps} 
\caption{Pair Distributions Functions obtained for a restitution 
coefficient $\epsilon = 0.80$, acceleration $\Gamma = 1.45$ 
and packing fraction $\phi = 0.30$, obtained from the original EDMD
(blue circles) and from a MC simulation with the derived effective 
potential (green crosses).}
\label{mc_e080_g145_p030} 
\end{figure} 

\begin{figure}
\includegraphics[width=9.0cm,angle=270]{mul2.ps} 
\caption{ Pair Distribution Functions obtained via Monte Carlo (green crosses), 
and EDMD (blue circles) methods. Panels (a) and (b) correspond to the minimum 
and maximum packing fractions used: $\phi$ = 0.20 and 0.40 respectively. 
Panels (c) and (d) correspond to the minimum and maximum restitution 
coefficients used: $\epsilon$ = 0.60 and 0.90. Panels (e) and (f) 
correspond to the minimum and maximum vertical accelerations used: 
$\Gamma$ = 1.05 and 1.45.}
\label{varias} 
\end{figure} 

\section{Conclusions and Final Comments}

We have simulated the quasi-2D granular gas produced 
by sinusoidal vertical shaking of a horizontal flat substrate
in which a set of spheres can bounce and roam. The
simulations are carried out using a fully 3D EDMD. 
Two-dimensional PDFs are obtained 
from the projected sphere center positions, for diverse sets of
parameters, and used as the input for an inversion 
procedure in order to get an effective pair potential 
that describes the apparent attractive interaction 
among pairs of particles. The inversion procedure 
involves the 2D Ornstein-Zernique equation using the 
Percus-Yevick closure; this of course assumes as hypothesis 
the thermodynamic equilibrium of the system. The
effective potentials thus obtained are later on used as
inputs for Monte Carlo simulations of a 2D gas of disks.
We have found very good agreement in the PDFs obtained from
both procedures. 
Then, a dissipative granular gas in which energy injection
maintains a steady state can be viewed as an equilibrium 
system of particles interacting through a potential well.
The Correlation at contact grows ---and consequently the 
effective pair interaction potential well becomes deeper---
as the restitution coefficient $\epsilon$ decreases,
as observed experimentally 
\cite{bordallo-favela_EPJE_2009} and in previous simulations 
\cite{perera-burgos_PRE_2010}. A fundamental exception to this trend 
does emerge, due to a resonant cooling down 
of the system that appears for the synchronization 
of the parabolic flights of the particles in the granulate,
given the appropriate combination of shaking amplitudes 
and restitution coefficients.

Among the issues that remain open with respect of the 
effective potentials discussed here, we would like to comment 
on two of them. The first one is
how far can we take this analogy for lower values of the forcing,
where condensation appears. The results found here
do show strong sensitivity of the effective potentials to density, 
and therefore it is not clear at all that, in a situation where
there is coexistence between a (crystallized) condensate and 
mobile and more diluted fluid, the same effective potential 
may be used for both. There is one other well known situation where 
effective forces should depend on local densities, and
still a single effective interaction potential manages to 
produce adequate phase diagrams; that is the case of 
depletion interactions in bidisperse mixtures of spherical 
colloids \cite{dijkstra_PRL_1998,dijkstra_PRL_1999}. 
For the present case we still 
do not have an answer. A second point of interest is whether or not
the high concordance obtained for the static behavior of the
two considered systems, as shown in the PDFs, will also 
extend to dynamical quantities, such as the diffusion coefficients
and Intermediate Scattering Functions. This of course will require 
the use of soft-disc simulations, since the definition of time scales
in Monte Carlo simulations is not completely obvious. 
We are now exploring these possible extension of the 
forced-dissipative to equilibrium analogy.


\begin{acknowledgments}
We acknowledge support from CONACyT though grant 221961.
S.\ V.-P.\ acknowledges the support from a CONACyT
Doctorate Fellowship 290611 and scholarship 290935.
\end{acknowledgments}

\newpage

\references

\bibitem{duran_2000}
J.\ Duran \emph{Sands, Powders, and Grains: An Introduction to the 
Physics of Granular Materials (Partially Ordered Systems)\/}, 
Springer (2000).

\bibitem{jaeger_RMP_1996}
H.\ M.\ Jaeger, S.\ R.\ Nagel and R.\ P.\ Behringer, 
Granular solids, liquids, and gases,
Rev.\ Mod.\ Phys. {\bf 68} 1259 (1996).

\bibitem{bizon_PRL_1998} 
Bizon, C.\ Shattuck, M. D.\ Swift, J. B.\ McCormick, W. D.\ Swinney, H. L.,
Patterns in 3D Vertically Oscillated Granular Layers: Simulation and Experiment, 
Phys.\ Rev.\ Lett. {\bf 80} 57 (1998).

\bibitem{roeller_PRL_2011} 
Roeller, K.\ Clewett, J. P. D.\ Bowley, R. M.\ Herminghaus, S.\ Swift, M. R., 
Liquid-Gas Phase Separation in Confined Vibrated Dry Granular Matter,
Phys.\ Rev.\ Lett. {\bf 107} 048002 (2011).

\bibitem{tatal_PRL_2000} 
Tata, B. V. R.\ Rajamani, P. V.\ Chakrabarti, J.\ Nikolov, A.\ Wasan, D. T.,
Gas-Liquid Transition in a Two-Dimensional System of Millimeter-Sized 
Like-Charged Metal Balls, 
Phys.\ Rev.\ Lett. {\bf 84} 3626 (2000).

\bibitem{knight_PRE_1998} 
Nowak, E. R.\ Knight, J. B.\ Ben-Naim, E.\ Jaeger, H. M.\ Nagel, S. R.,
Fluctuations in the Density of a Granular Material During Vibration,
Phys.\ Rev.\ E {\bf 57} 1971 (1998).

\bibitem{eshuis_Phys_Flu_2007}
Eshuis, P.\ van der Weele, K.\ van der Meer, D.\ Bos, R.\ Lohse, D.,
Phase diagram of vertically shaken granular matter,
Phys.\ Fluids {\bf 19} 123301 (2007).

\bibitem{reis_PRL_2007} 
P.\ M.\ Reis, R.\ A.\ Ingale and M.\ D.\ Shattuck,
Caging Dynamics in a Granular Fluid, 
Phys.\ Rev.\ Lett. {\bf 98} 188301 (2007).

\bibitem{olafsen_PRL_1998} 
J.\ S.\ Olafsen and J.\ S.\ Urbach,
Clustering, Order, and Collapse in a Driven Granular Monolayer, 
Phys.\ Rev.\ Lett. {\bf 81} 4369 (1998).

\bibitem{olafsen_PRE_1999} 
J.\ S.\ Olafsen and J.\ S.\ Urbach,
Velocity distributions and density fluctuations in a granular gas, 
Phys.\ Rev.\ E {\bf 60} R2468 (1999).

\bibitem{melby_JPCM_05}
P.\ Melby \emph{et al.\/},
The dynamics of thin vibrated granular layers,
J.\ Phys.: Condens.\ Matter {\bf 17} S2689 (2005).

\bibitem{vega_reyes_2008}
F.\ Vega-Reyes and J.\ S.\ Urbach, 
Effect of inelasticity on the phase transitions of a thin vibrated granular layer, 
Phys.\ Rev.\ E {\bf 78} 051301 (2008).


\bibitem{olafsen_PRL_2005} 
J.\ S.\ Olafsen and J.\ S.\ Urbach, 
Two-Dimensional Melting Far from Equilibrium in a Granular Monolayer,
Phys.\ Rev.\ Lett. {\bf 95} 098002 (2005).

\bibitem{reis_PRL_2006} 
P.\ M.\ Reis, R.\ A.\ Ingale and M.\ D.\ Shattuck, 
Crystallization of a Quasi-Two-Dimensional Granular Fluid,
Phys.\ Rev.\ Lett. {\bf 96} 258001 (2006).

\bibitem{alder_PR_1962} 
B.\ J.\ Alder and T.\ E.\ Wainwright,
Phase Transition in Elastic Disks,
Phys.\ Rev. {\bf 127} 359 (1962).

\bibitem{engel_PRE_2013} 
Engel, M.\ Anderson, J. A.\ Glotzer, S. C.\ Isobe, M.\ Bernard, 
E. P.\ Krauth, W., Hard-disk equation of state: First-order liquid-hexatic 
transition in two dimensions with three simulation methods,
Phys.\ Rev.\ E {\bf 87} 042134 (2013).

\bibitem{nie_EPL_2000}
X.\ Nie, E.\ Ben-Naim and S.\ Y. Chen,
Dynamics of vibrated granular monolayers,
Europhys.\ Lett. {\bf 51} 679 (2000).

\bibitem{perez-angel_PRE_2011} 
G.\ P\'erez-Angel and Y.\ Nahmad-Molinari,
Bouncing, rolling, energy flows, and cluster formation in a two-dimensional vibrated granular gas,
Phys.\ Rev.\ E {\bf 84} 041303 (2011).



\bibitem{perera-burgos_PRE_2010}
J.\ A.\ Perera-Burgos, G.\ P\'erez-\'Angel and Y.\ Nahmad-Molinari,
Diffusivity and weak clustering in a quasi-two-dimensional granular gas,
Phys.\ Rev.\ E {\bf 82} 051305 (2010)

\bibitem{bordallo-favela_EPJE_2009}
R.\ A.\ Bordallo-Favela \emph{et al.\/},
Effective potentials of dissipative hard spheres in granular matter,
Eur.\ Phys.\ J.\ E {\bf 28} 395 (2009).

\bibitem{derjaguin_APCURSS_1941}
B.\ Derjaguin and L.\ Landau,
Theory of the stability of strongly charged lyophobic sols and of the adhesion 
of strongly charged particles in solutions of electrolytes, 
Acta Phys.\ Chim.\ URSS {\bf 14} 633 (1941).

\bibitem{verwey_1948}
J.\ E.\ Verwey and J.\ T.\ G.\ Overbeek, 
\emph{Theory of the stability of lyophobic Colloids\/}, 
Elsevier, Amsterdam (1948).

\bibitem{hiemenz_1997} 
P.\ C.\ Hiemenz and M.\ Rjagopalan,
\emph{Principles of Colloid and Surface Chemistry\/},
3rd edition,
Marcel Dekker, New York (1997).

\bibitem{asakura_JCP_1954}
S.\ Asakura and F.\ Oosawa, 
On Interaction between Two Bodies Immersed in a Solution of Macromolecules,
J.\ Chem.\ Phys. {\bf 22} 1255 (1954).

\bibitem{asakura_JPS_1958}
S.\ Asakura and F.\ Oosawa (1958) 
Interaction between Particles Suspended in Solutions of Macromolecules, 
J.\ Pol.\ Sci. {\bf 33} 183 (1958).




\bibitem{ciamarra_PRL_2006}
M.\ P.\ Ciamarra, A.\ Coniglio and M.\ Nicodemi,
Dynamically Induced Effective Interaction in Periodically Driven Granular Mixtures,
Phys.\ Rev.\ Lett. {\bf 97} 038001 (2006).

\bibitem{poschel_2005}
T.\ P\"oschel and T.\ Schwager,
\emph{Computational Granular Dynamics: Models and Algorithms\/}, 
Springer (2005).

\bibitem{rapaport_2004}
D.\ C.\ Rapaport, 
\emph{The Art of Molecular Dynamics Simulation},
2nd edition, Cambridge University Press (2004).

\bibitem{perez_2008_Pramana}
G.\ Perez,
Numerical simulations in granular matter: The discharge of a 2D silo,
Pramana - J.\ Phys. {\bf 70} 989 (2008).

\bibitem{schwager_2008_EPJE}
T.\ Schwager, V.\ Becker and T. P\"oschel,
Coefficient of tangential restitution for viscoelastic speheres,
Eur.\ Phys.\ J.\ E {\bf 27} 107 (2008).

\bibitem{mcnamara_PRE_1996} 
S.\ McNamara and W.\ R.\ Young,
Dynamics of a freely evolving, two-dimensional granular medium, 
Phys.\ Rev.\ E {\bf 53} 5089 (1996).

\bibitem{mcnamara_PRE_1998}
S.\ McNamara and S.\ Luding,
Energy flows in vibrated granular media,
Phys.\ Rev.\ E {\bf 58} 813 (1998).

\bibitem{gonzalez_EPJST_2009}
S.\ Gonz\'alez, D.\ Risso and R.\ Soto, 
Extended event driven molecular dynamics for simulating dense granular matter,
Eur.\ Phys.\ J.\ Special Topics {\bf 179} (2009) 

\bibitem{fisher_JMP_1964} 
M.\ E.\ Fisher,
Correlation Functions and the Critical Region of Simple Fluids,
J.\ Math.\ Phys. {\bf 5} 944 (1964).

\bibitem{textbooks} 
The PY and other closures for the OZ equation are well covered in 
many textbooks, see for instance:
D.A.\ McQuarrie, \emph{Statistical Mechanics\/}, 
(University Science Books, Sau\-sa\-li\-to, 2000);
D.L.\ Goodstein,
\emph{States of Matter\/},
(Dover Publications, New York, 1985)

\bibitem{hill_AMM_00} 
J.\ M.\ Hill, M.\ J.\ Jennings, D.\ V.\ To and K.\ A.\ Williams,
Dynamics of an elastic ball bouncing on an oscillating plane and the oscillon, 
Appl.\ Math.\ Model. {\bf 24} 715 (2000).

\bibitem{dijkstra_PRL_1998}
M.\ Dijkstra, R.\ van Roij and R.\ Evans,
Phase Behavior and Structure of Binary Hard-Sphere Mixtures.
Phys.\ Rev.\ Lett. {\bf 81} 2268 (1998).

\bibitem{dijkstra_PRL_1999}
M.\ Dijkstra, R.\ van Roij and R.\ Evans,
Direct Simulation of the Phase Behavior of Binary Hard-Sphere Mixtures:
Test of the Depletion Potential Description,
Phys.\ Rev.\ Lett. {\bf 82} 117 (1999).





\end{document}
