Asymmetric and Long Range Interactions in Shaken Granular Media

TL;DR

This study uses computational modeling to reveal long-range attractive forces and interaction dynamics between intruders in a vibrated granular medium, highlighting the influence of shaking parameters and packing density.

Contribution

It introduces a detailed computational analysis of long-range interactions and effective forces between intruders in vibrated granular media, emphasizing the role of shaking parameters.

Findings

Emergence of long-range attractive forces between intruders.

Interaction forces depend on shaking angle, amplitude, and packing density.

No density gradients, but additive kinetic energy gradients observed.

Abstract

We use a computational model to investigate the emergence of interaction forces between pairs of intruders in a horizontally vibrated granular fluid. The time evolution of a pair of particles shows a maximum of the likelihood to find the pair at contact in the direction of shaking. This relative interaction is further studied by fixing the intruders in the simulation box where we identify effective mechanical forces, and torques between particles and quantify an emergent long range attractive force as a function of the shaking relative angle, amplitude, and the packing density of grains. We determine the local density and kinetic energy profiles of granular particles along the axis of the dimer to find no gradients in the density fields and additive gradients in the kinetic energies.

Tables

Table 1. Table 1: Effective temperature obtained by minimizing χ 2 superscript 𝜒 2 \chi^{2} combining measures of the mechanical formation energy and the for different densities, shaking amplitudes, and orientation angles. Values with missing α 𝛼 \alpha correspond to the angular averaged formation energies and probabilities.

$A_0$	$\varphi$	$\alpha$	$k_B{T}_{\text{eff}}$	$A_0$	$\varphi$	$\alpha$	$k_B{T}_{\text{eff}}$
$0.75$	0.6	$0$	0.35	$0.75$	0.75	$0$	0.24
0.75	0.6	$\pi /4$	0.64	0.75	0.75	$\pi /4$	0.42
0.75	0.6	$\pi /2$	7.4	0.75	0.75	$\pi /2$	1.21
0.75	0.6	$-$	0.61	0.75	0.75	$-$	0.46
1.5	0.6	$0$	0.35	1.5	0.75	$0$	0.16
1.5	0.6	$\pi /4$	0.46	1.5	0.75	$\pi /4$	0.26
1.5	0.6	$\pi /2$	$-0.47$	1.5	0.75	$\pi /2$	$0.95$
1.5	0.6	$-$	$0.39$	1.5	0.75	$-$	$0.21$

Equations

$$m_i\frac{d^2{\mathbf{x}}_i}{d{t}^2}=-{\gamma }_{s,i}\left(\frac{d{\mathbf{x}}_i}{dt}-{\mathbf{v}}_s(t)\right)+{\mathbf{F}}_i^c+{\mathbf{F}}_i^d+{\mathbf{F}}_i^r,$$

$$M\frac{d^2{\mathbf{X}}_i}{d{t}^2}=-{\gamma }_s\frac{d{\mathbf{X}}_i}{dt}+{\mathbf{F}}_i^c+{\mathbf{F}}_i^d+{\mathbf{F}}_i^r,$$

$${\mathbf{F}}_{ij}^c=k{\xi }_{ij}\theta \left({\xi }_{ij}\right){\hat{\mathbf{n}}}_{ij},$$

$${\mathbf{F}}_{ij}^d=\gamma \left({\mathbf{v}}_{ij}\cdot {\mathbf{n}}_{ij}\right)\theta \left({\xi }_{ij}\right){\mathbf{n}}_{ij},$$

$$F_r(d)=⟨\frac{{\mathbf{F}}_b({\mathbf{X}}_b)-{\mathbf{F}}_a({\mathbf{X}}_a)}{2}\cdot \hat{\mathbf{r}}⟩,$$

$$F_{t}(d)=\Big{\langle}\big{(}\bm{F}_{b}(\bm{X}_{b})-\bm{F}_{a}(\bm{X}_{a})\big{)}\cdot\hat{\bm{t}}\Big{\rangle},$$

$$F_r=F(\varphi )\left[\frac{\theta (A_0-A_0^c(\varphi ))}{{\sigma }_g}+B_0(\varphi )\right]{\left(\frac{{\sigma }_g}{r}\right)}^2{\cos }^2\alpha .$$

$${\mathcal{U}}_{\alpha }(d)={\int }_{\infty \to d}{\mathbf{F}}_{int}\cdot d\mathbf{l}={\int }_d^{\infty }{F}_r(r,\alpha )dr.$$

$$\delta e_k(h)=\frac{e_k(h)-⟨e_k⟩}{⟨e_k⟩},\delta t_k(h)=\frac{t_k(h)-⟨t_k⟩}{⟨t_k⟩}$$

$$\frac{d}{dt}{\mathbf{r}}_i(t)=\frac{1}{m_i}{\mathbf{p}}_i(t),$$

$$\frac{d}{dt}{\mathbf{p}}_i={\mathbf{F}}_i-{\gamma }_i\left(\frac{{\mathbf{p}}_i}{m_i}-{\mathbf{v}}_s\right)+{\mathbf{f}}_i,$$

$$\left(\frac{d}{dt}+\frac{{\gamma }_i}{m_i}\right){\mathbf{p}}_i={\mathbf{F}}_i+{\gamma }_i{\mathbf{v}}_s+{\mathbf{f}}_i={\mathbf{\Phi }}_i+{\mathbf{f}}_i.$$

$${\mathbf{p}}_i^{(h)}\left(t+t_0\right)={\mathbf{p}}_i(t_0)\exp \left[\frac{{\gamma }_i}{m_i}\left(t-t_0\right)\right],$$

$$\frac{d}{dt}{\mathbf{\pi }}_i(t)=\exp \left(\frac{{\gamma }_i}{m_i}t\right)\left({\mathbf{\Phi }}_i+{\mathbf{f}}_i\right),$$

$${\mathbf{\pi }}_i\left(t+\Delta t\right)$$

$${\mathbf{p}}_i(t+\Delta t)={\Gamma }_i{\mathbf{p}}_i(t)+\tilde{\mu }{m}_i{\mathbf{\Phi }}_i+{\tilde{D}}_i{\mathbf{\zeta }}_i,$$

$${\mathbf{r}}_i(t+\Delta t)={\mathbf{r}}_i(t)+\frac{3{\mathbf{p}}_i(t)+{\mathbf{p}}_i(t+\Delta t)}{4m_i}\Delta t.$$

$$e_k=\frac{1}{N}\sum _i\frac{m_i}{2}{\left|{\mathbf{v}}_i\right|}^2,t_k=\frac{1}{N}\sum _i\frac{m_i}{2}{\left|{\mathbf{v}}_i-{\mathbf{V}}_{cm}\right|}^2,$$

Full text

††preprint: APS/123-QED

I Introduction

Granular matter has attracted recent attention because of its intrisically different nature as compared to equilibrium solid, liquid and gas states Jaeger, Nagel, and Behringer (1996); due to its ubiquity in the industry in the form of plastic beads, sand, as well as, edible grains, i.e., coffee, wheat; jamming and clogging processes, and stress distribution which lead to the formation of stable structures. Jamming in granular matter has been widely studied in the past Liu and Nagel (1998); Cates et al. (1998).

Shaken granular matter are actuated systems out-of-equilibrium due to a steady flux of energy from the container to the grains that is finally dissipated through contact interactions. Striking phenomena occur in mixtures of particles where shaking leads to species segregation ranging from clusters to stripes Ottino and Khakhar (2000); Kudrolli (2004); Aranson and Tsimring (2006). The best known manifestation of this granular separation is the so called Brazil nut effect Godoy et al. (2008); Ciamarra et al. (2006); Sanders et al. (2004). For horizontally driven matter, gravity is no longer a relevant parameter, and a mixture can phase separate into stripes orthogonal to the shaking direction Mullin (2000); Reis and Mullin (2002); Reis, Sykes, and Mullin (2006); Pica Ciamarra, Coniglio, and Nicodemi (2007), or even form clusters for swirling shakings Aumaître, Kruelle, and Rehberg (2001). Understanding the ability of granular mixtures to demix is of direct industrial relevance and can help to explain diverse phenomena, such as stratification in terrestrial environments or in asteroid or planetoid formation Opsomer et al. (2014).

As a paradigmatic example of the qualitative change that granular systems experience under shaking we have the case of granular flow though a bottleneck – like the flow of sand in an hourglass. Under the sole effect of a constant force, such as gravity, the granular flow can be spontaneously interrupted Masuda, Nishinari, and Schadschneider (2014); Lozano et al. (2012) if the exit is not wide enough. On the contrary, for granular matter under shaking the nature of the flux interruption dramatically changes and the system, once clogged, may spontaneously unclog Zuriguel et al. (2014).

In the non-equilibrium state induced by shaking, granular matter has been reported to fluidize for vertical and horizontal Salueña, Pöschel, and Esipov (1999); Pöschel, Schwager, and Salueña (2000); Ristow, Straßburger, and Rehberg (1997) forcing. Once granular matter fluidizes, large density fluctuations have been reported and such a situation leads to the opening of Casimir-like scenarios Kardar and Golestanian (1999) for the effective granular interactions Cattuto et al. (2006).

Band segregation for a binary mixture of horizontally agitated granular matter has been analyzed for attractive and anisotropic pairwise interactions between inclusions Ciamarra, Coniglio, and Nicodemi (2006). However, the details of the interactions were not directly addressed. Position tracking experiments of metallic spheres immersed in a dense poppy-seed granular monolayer have revealed not only an attraction between pair but an aligning interaction relative to the shaking direction Lozano et al. (2015). Following this experimental evidence we present this model that captures and quantifies the pair interaction between intrusive particles.

This paper is structured as follows. In Section II we present a model for both grains, and inclusions under horizontal vibration. In Section III, we study the free motion of an intruder pair pair in a shaken granular bed and identify the probability to locate particles as a function of their separation. In Section IV, we fix the inclusion pair in order to quantify the emergent effective radial and tangential forces; the mechanical formation energy of the dimers is then compared to the pair distance distribution previously obtained in Section IV; and finally the effect of the pair in the kinetic energy of the granular bed.

II Model

Recent experiments have studied the dynamics and structure of a pair of phosphor-bronze spheres immersed in a horizontally shaken monolayer of poppy seeds. Poppy seeds are kidney shaped granular particles with typical diameters ranging $0.5~{}\text{mm}$ to $1~{}\text{mm}$ with material density $\rho=0.2~{}\text{g cm}^{-3}$. Poppy seeds have a wide contact area with the surface of the tray, so their friction is large, and periodically displace with the vibrated tray. Bronze spheres, on the contrary, are smooth, and uniform in size with a diameter of $1.5~{}\text{mm}$, and material density Reis et al. (2004) $\rho=8.8~{}\text{g cm}^{-3}$.

Spheres easily rotate on the vibrated flat surface of the system Reis et al. (2004), whereas grains have a tendency to follow the moving tray. Inspired by this situation, and as sketched in Fig. 1, we consider a simplified two dimensional model for the dynamics of grains and inclusions which captures the essential features of the forced grain monolayer, and develop an integration scheme to computationally solve the equations of motion of the granular system.

II.1 Equations of motion

We treat grains as $N$ disks at positions $\bm{x}$ of diameter $\sigma$, drawn from a uniform distribution with average $\left\langle\sigma\right\rangle=\sigma_{g}$ and $10\%$ dispersion to account for experimental polydispersivity Lozano et al. (2015), in a 2d periodic box of side $L$. We define an averaged packing fraction of the monolayer as $\phi=N\pi\langle\sigma^{2}\rangle/(4L^{2})$. Each grain with a mass that depends on their diameter as $m=m_{0}\left(\sigma/\sigma_{g}\right)^{3}$, where $m_{0}$ is the mass of a grain of size $\sigma_{g}$. Grains are periodically driven with the tray velocity $\bm{v}_{s}(t)$. Bronze spheres, inclusions, are modeled as disks of diameter $\sigma=3\sigma_{g}/2$, and mass $M=50\ m_{0}$ at positions $\bm{X}$. The solid rotation of inclusions on the moving plate is compatible with the description of the position vector $\bm{X}$ relative to the laboratory frame, i.e., we assume that inclusions are not being displaced by the tray motion $\bm{v}_{s}(t)$. For each particle, either $\bm{x}$ or $\bm{X}$ we present its corresponding equation of motion,

[TABLE]

[TABLE]

where we introduce an oscillatory velocity of the tray, $\bm{v}_{s}(t)=A_{0}\omega\sin\omega t\hat{\bm{e}}_{x}$, characterized by its displacement amplitude, $A_{0}$, and frequency $\omega$. The motion of the grains relative to the tray introduces a friction force whose magnitude $\gamma_{s}\left(d\bm{x}_{i}/dt-{\bm{v}}_{s}\right)$ depends on the velocity of the grain relative to the tray, and a dissipation constant $\gamma_{s}$ proportional to the contact area of the disk $\sigma^{2}$.

The excluded volume interaction among particles, as well as the energy dissipated in particle collisions, is described through a linear spring-dashpot model for the interactions between particles. The total conservative and dissipative forces on particle $i$ are pairwise additive, ${\bm{F}}_{i}^{(c/d)}=\sum_{j}{\bm{F}}_{ij}^{(c/d)}$ that activate when grains overlap Hidalgo et al. (2009); Pöschel, Schwager, and Salueña (2000); Pöschel and Schwager (2005). The degree of compression between two disks $i$ and $j$ is quantified by $\xi_{ij}=\left|{\bm{x}}_{ij}\right|-\left(\sigma_{i}+\sigma_{j}\right)/2$. Accordingly, the elastic pairwise conservative force reads

[TABLE]

with $k$ an elastic constant (a material dependent parameter), and $\hat{\bm{n}}_{ij}$ the centre-to-centre unit vector. $\bm{F}_{ij}^{c}$ does not vanish only when the two disks overlap, i.e., $\xi_{ij}>0$ as accounted to by the Heavisde function, $\theta(\xi)$.

The dissipative pairwise interaction, responsible of the energy loss, reads

[TABLE]

with $\gamma$ a dissipation constant, $\bm{v}_{ij}=\bm{v}_{i}-\bm{v}_{i}$, and its direction opposes the direction of the relative velocity between pairs. This dissipation mechanism between pairs may lead to configurations of small overlaps with the elastic and dissipative forces canceling each other. In these situations pairs of particles tend to remain close to each other. To avoid such computational artifacts, as proposed in Pöschel and Schwager (2005), we use the total force of interaction between pairs to be always either repulsive or zero $|\bm{F}_{ij}|=\max(0,|\bm{F}^{c}_{ij}+\bm{F}^{d}_{ij}|)$.

Additionally, we introduce random forces, $\bm{F}^{r}$, to account for irregularities and vertical collisions. The source of noise is Gaussian for each component with zero mean and variance $\left\langle\hat{\bm{e}}_{\alpha}\cdot\bm{F}^{r}_{i}(t)\hat{\bm{e}}_{\beta}\cdot\bm{F}^{r}_{j}(t^{\prime})\right\rangle=2\Lambda_{\alpha}\delta_{ij}\delta_{\alpha\beta}\delta\left(t-t^{\prime}\right)$. The $\Lambda_{\alpha}$ parameter accounts for the asymmetry in the noise strength in the parallel and perpendicular directions relative to the tray shaking.

In the following sections we will consider the tray oscillation period $\tau^{-1}=2\pi\omega$, the average grain size $\sigma_{g}$, and its mass $m_{0}$ as time, length and mass units respectively. In terms of these magnitudes we define the elastic constant $k$, the dissipation constants $\gamma$, and $\gamma_{s}$, and $\Lambda_{x}$ with standard values present in the literatureSalueña, Pöschel, and Esipov (1999). The integration of the equations is performed by a stochastic integrator inspired on Allen and Tildesley (1989) and detailed in Appendix A.

In Appendix B we characterize a bed of granular particles and analyze their kinetic energies. We observe that the packing fraction of the grain monolayer has an impact in the inclusion energy loss above $\phi=0.6$ whereas beyond $\phi\approx 0.75$ the dynamics of the grains slows down and the the displacement statistics of grains show the caging effects introduced by the formation of a quasi-crystalline structure of the grains. In the rest of this paper we consider granular beds within the density range $\phi\in[0.6,~{}0.75]$.

III Free moving dimer

In order to analyze the emergent features induced by a forced granular system on embedded intruders, we will analyze the dynamics of a pair of intruders. To this end, we consider the free evolution of a dimer of inclusion particles in a granular bed. In subsection III.1 we define the relative coordinates that describe the pair motion. In subsection III.2 we extract the relative arrangement of the inclusion pair from its evolution as a function of the granular shaking amplitudes and packing fractions. Finally, in III.3 we analyze the the averaged probability distribution of the pair relative positions to measure the pair separated a distance $d$ and the relative orientation for both touching and distant pairs.

III.1 An inclusion dimer

In a periodically driven granular system the energy is introduced by the external forcing through both $A_{0}$, and $\omega$. On the one hand, the shaking amplitude $A_{0}$ introduces a length scale that may lead to structural system deformations of typical length $\sim 2A_{0}$. On the other hand, the shaking frequency $\omega$ modifies the dissipation rate of the granular system. In this paper we control the tray vibration by changing the shaking amplitude, $A_{0}$, and fixing its frequency.

A doublet of intruders of diameter $\sigma$, $\{\bm{X}_{a},~{}\bm{X}_{b}\}$, defines a dimer whose center to center vector is ${\bm{r}}=\bm{X}_{a}-\bm{X}_{b}$. The position vector ${\bm{r}}=(r_{x},~{}r_{y})$ is equally described by its modulus $r$ and angle $\alpha$ relative to the shaking direction, $\bm{r}=r\left(\cos\alpha,~{}\sin\alpha\right)$. We define a “parallel” (“perpendicular”) configuration of the dimer for a configuration $\alpha\approx 0$ ($\pi/2$). For simplicity we also use the distance between the surfaces of inclusions $d=r-\sigma$ to define the dimer distance.

III.2 Probability Landscapes

In order to understand the dynamics of a dimer of intruders in a shaken granular bed we have prepared a large collection of initial conditions with $d\in\left[4,7\right]$ and $\alpha\in\left[0,2\pi\right]$ to sample between $10^{6}$ ($\phi=0.6$) and $6\cdot 10^{6}$ ($\phi=0.75$) cycles for systems of box size $L=32\sigma_{g}$ at shaking amplitudes $A_{0}=0.75\sigma_{g}$ and $A_{0}=1.5\sigma_{g}$.

We define $\mathcal{P}(r_{x},~{}r_{y})$ as the probability to measure the dimer in a configuration with $r_{x}\in\left[r_{x},r_{x}+\Delta x\right]$ and $r_{y}\in\left[r_{y},~{}r_{y}+\Delta y\right]$, where $\Delta x=\Delta y=0.05\sigma_{g}$, that results from simulations sampling at each shaking cycle. From this probability, we can extract $-\ln\mathcal{P}\left(r_{x},r_{y}\right)$, which identifies the most probable configurations of the dimer freely evolving in the granular bed. If the system were in equilibrium, the logarithm of the probability $-\ln\mathcal{P}$ would be proportional to the energy landscape $\mathcal{U}(r_{x},r_{y})$ of the system. We leave the comparison between $\mathcal{P}$ and the dimer mechanical formation energy, $U$, for section IV.3.

We associate minima of the magnitude $-\ln P$ to more probable configurations (effective attractions), and maxima to more unlikely configurations (effective repulsions). Probability landscapes $-\ln\mathcal{P}(r_{x},r_{y})$, see Fig. 2, identify $r\approx\sigma$ as the most probable pair distance for all studied granular beds. Additionally, all four baths show a high degree of anisotropy at contact, $r\approx\sigma$, where $-\ln\mathcal{P}(\sigma,\alpha)$ is always smaller for $\alpha\approx 0$ than for $\alpha\approx\pi/2$. This clear asymmetry shows a strong tendency of dimer to align in the shaking direction.

The structure of the granular bed is clearly manifested in the probability landscape as a sequence of concentric oscillating rings. At long distances from contact, contour lines manifest the anisotropy of the interaction. Moreover, the probability landscape decays faster along the shaking direction than perpendicular to it.

Fig. 3 further details the probability landscape for $A_{0}=1.5\sigma_{g}$, and packing fraction $\phi=0.75$. The granular structure is rapidly lost in the parallel direction whereas it shows stronger persistence for configurations where the dimer is oriented perpendicular to the tray shaking direction. At long distances, we observe that perpendicular configurations are favored and an island of repulsion emerges centered at $r=8\sigma_{g}$, and $\alpha=0$.

III.3 Marginal probability densities

We can provide a more compact description of the distance and angular dependence that characterize the dimer interaction through the corresponding marginal probability densities. These quantities are also more attractive experimentally because they require less statistics. Specifically, we quantify the distance behavior of $\mathcal{P}$ through $P(r)$, the averaged angular distribution, obtained integrating $\mathcal{P}(r,\alpha)$ over the angular configurations of the dimer.

Fig. 4 displays the angular averaged probability density, $P(d)$, as a function of $d=r-\sigma$, so that $d=0$ corresponds to the contact position $r=\sigma$. We identify three different features in $P(d)$. First, at contact ($d=0$), $-\ln P(d)$ reaches its lowest value; hence it is the preferred and most probable configuration. Second, the oscillations of $-\ln P(d)$ at short distances clearly show the structure of the bed since the periodicity corresponds to the granular size. The surroundings of minima are locally stable and correspond to an integer number of grains in the internal region of the inclusion dimer. Third, the oscillations are superimposed on a monotonous decay towards zero. The decay is stronger the more compact the system is – the signal is stronger at $\phi=0.75$. Within the available range of distances, we have tested that the decay in $-\ln P(d)$ is compatible with an algebraic tail -\ln\big{(}P(d)\big{/}P(\infty)\big{)}\sim d^{-1}, shown in dashed lines in Fig. 4.

To quantify anisotropy we compute the conditional probability densities $P_{<}(r^{\star};\alpha)$, and $P_{>}(r^{\star};\alpha)$. Obtained as integrals of $\mathcal{P}(r,\alpha)$ over $r$ in the ranges $r\in\left(\sigma,\sigma+r^{\star}\right)$, and $r\in\left(\sigma+r^{\star},\infty\right)$ for $P_{<}$ and $P_{>}$ respectively. $P_{<}(r^{\star};\alpha)$ quantifies the anisotropy of the configurations of the dimer at short distances whereas $P_{>}(r^{\star};\alpha)$ quantifies long distance anisotropy.

The angular distributions displayed in Fig. 5 have been obtained for $r^{\star}=2\sigma_{g}$. At short distances, $P_{<}(2\sigma;\alpha)$, we observe a preference for parallel configurations. The anisotropy in the configurations has a major dependence on the shaking amplitude, while it shows a mild dependence on packing density. At larger distances, $d>2\sigma_{g}$, on the contrary, the dimer tends to orient perpendicular to the shaking direction. In this case the degree of anisotropy is relatively weak except for large densities and shaking amplitudes, $\phi=0.75$ and $A_{0}=1.5\sigma$, where the anisotropy increases revealing the instability island clearly appreciable in the 2d heat map in Fig. 3.

If the system were to be in equilibrium, the interpretation of the probabilities would be straightforward. The minimum, and slow decay of the radial probability $-\ln P(d)$ would translate into an energy minimum at $d=0$, and an emergent attractive long range interaction $U\sim d^{-1}$. The results in the anisotropy shown in Fig. 5(a) would imply the emergence of an aligning torque at short distances; a torque that would align the dimer towards the shaking direction. However, the system is not in equilibrium, and thus we cannot relate the probability of the dimer configurations to the interaction energies in the system. In the following section we compute mechanically the effective interaction forces and formation energies of a dimer following similar procedures to the ones used by Harder et al. (2014); Ni, Cohen Stuart, and Bolhuis (2015); Zaeifi Yamchi and Naji (2017).

IV Fixed inclusions

In this section we quantify the mechanical interactions between inclusions that arise from the granular shaking. Fist in IV.1 we fix the inclusions at relative coordinates $(r,~{}\alpha)$ and measure their relative mean radial and tangential interactions. Second in IV.2 we propose a model that captures the long distance behavior of the radial force as a simple function of the shaking amplitude $A_{0}$, the relative angle $\alpha$, and the inclusions separation distance $d$. Third in IV.3 we define the mechanical formation energy of the dimer and compare it to the configurations’ probabilities for the freely evolving dimers measured in Section III.3 Finally in IV.4 we compute local profile measures of the granular bed along the axis of the dimer to provide further insight in the mechanism underlying the effective interactions between inclusions.

IV.1 Radial and tangential forces

We fix the pair of inclusions at positions $\bm{X}_{a}$ and $\bm{X}_{b}$ but we let the velocity of the inclusions evolve in time since particle collisions depend on the relative velocity between interacting particles. The forces given by the granular bed on the inclusions are $\bm{F}_{a}$ and $\bm{F}_{b}$ and define the effective interaction between inclusions. Computing the relative force $\bm{F}_{b}(\bm{X}_{b})-\bm{F}_{a}(\bm{X}_{a})$ we extract the relevant information of the relative force acting on the dimer. Finally we project along the radial and tangential directions,

[TABLE]

[TABLE]

where unit vectors $\hat{\bm{r}}$ and $\hat{\bm{t}}$ are an orthonormal basis of $\mathbb{R}^{2}$ with $\hat{\bm{r}}=\left(\bm{X}_{b}-\bm{X}_{a}\right)/\left|\bm{X}_{b}-\bm{X}_{a}\right|$ and, after rotating $\pi/2$ we obtain $\hat{\bm{t}}=\mathcal{R}_{\pi/2}\hat{\bm{r}}$. The brackets $\left\langle\cdot\right\rangle$ denote the average of the magnitude $(\cdot)$ over cycles and independent realizations of the system.

Inclusions have been placed at configurations with surface to surface distance the range $d\in\left[0,8\right]$ and angles relative to shaking $\alpha\in\left[0,\pi/2\right]$. We average on various sets with fixed parameters $\alpha,d,A_{0},\phi$ to extract the mean values $F_{r}$, and $F_{t}$.

Within this definition of $\hat{\bm{r}}$, and $\hat{\bm{t}}$, the projection $F_{r}$ defines the behavior of the interaction in the radial direction with attractive effective forces captured by negative values of the relative force, $F_{r}<0$, and repulsive otherwise. The transverse effective force, $F_{t}$, indicates the emergence of a neat torque acting on the dimer with clock-wise torques identified by $F_{t}<0$, and counterclock-wise otherwise. For our dimer configurations in the first quadrant of the plane, clock-wise tangential forces align the dimer in the shaking direction.

The relative radial force $F_{r}$ in Fig. 6 reflects that the effective interaction for a pair of inclusions at contact is attractive. We observe that the magnitude of the force at contact also depends on the orientation of the dimer relative to the shaking direction. The interaction is stronger the more aligned –to the shaking– the dimer is. As distance increases, the force reaches a local maximum before $d=\sigma_{g}$ and then $F_{r}$ oscillates with an amplitude that reflects the structure of the grain suspension.

The magnitude and sign of the radial force is sensitive to the orientation of the dimer with respect to the tray shaking direction, $\alpha$. For perpendicular configurations the decaying oscillations in the force reflect the structure of the granular bed for long distances as it is clear from the persistence of the force oscillations. As the angle is reduced, $\alpha<\pi/2$, the oscillation structure in the radial force is gradually destroyed at long distances while a neat attractive force remains for several inclusion diameters. Overall, shaking destroys the structure of the granular bed and induces an effective attraction that increases as the dimer aligns relative to the shaking.

In addition to the radial interaction, a torque captured as a relative tangential force $F_{t}$ appears for configurations with neither $\alpha\neq 0$, nor $\alpha\neq\pi/2$. We have observed a maximum value in the tangential force strength for $\alpha=\pi/4$, which we select to display in Fig. 6. This non-central force is present in the system for short separation distances, $d\lesssim 2\sigma_{g}$ and orients the pair along the shaking direction. At $\phi=0.75$, and $A_{0}=1.5$ a finite dealigning torque clearly emerges at long distances $d\gtrsim 2\sigma_{g}$. The nature of the tangential force shown in Fig. 6 is compatible with the anisotropy found in the probability $\mathcal{P}(r_{x},r_{y})$ for moving inclusions reported in Figs. 2 and 5 in the previous section III.

IV.2 Long distance behavior of the radial force

Fig. 4 shows that the effective force between inclusions generically displays an attractive, long-range interaction in length scales larger than $10\sigma_{g}$. Even though we are not able to analyze the functional dependence of this asymptotic decay over several decades, the observed decay is consistent with an algebraic dependence. In granular media, long range forces arising from their non-equilibrium fluctuations have already been described in vibrated systems Cattuto et al. (2006).For the weak decaying signal of the radial force measure for configurations aligned close to the perpendicular direction an algebraic fitting curve of the form $F_{2}(\sigma_{g}/d)^{2}$ has proven to be numerically more robust than an exponential decay which is largely interfered by the strong oscillations, see Fig. 7(a). For an algebraic fitting, the strength of the long range interaction is given by the force $F_{2}$ and presented in Fig. 7.

To capture the full dependence on $\alpha$, and $A_{0}$ on the long range interaction amplitude, $F_{2}$, we have systematically swept the space of parameters $\alpha$, and $A_{0}$. The structure of the bed of grains is appreciable in both Fig. 7 (a) and (b); at perpendicular configurations the emergent long range interaction vanishes, $F_{2}\approx 0$. At low values of the shaking amplitude, the grain bed structure persists at lower angles –denoted in dashed lines. As the angle of the dimer decreases the effect of the external forcing increases and the local granular structure is lost. In this regime a clear interaction between particles develops, and increases as the dimer aligns towards the forcing direction.

To quantify the dependence on the orientation $\alpha$ we present $F_{2}$ as a function of $\cos^{2}\alpha$, $F_{2}=B_{2}\cos^{2}\alpha+B_{0}$, the lowest power of $\cos\alpha$ allowed by the dimer and shaking symmetry. Results in Fig. 7(c) confirm a quadratic dependence of the cosine on $F_{2}$. The analysis of $B_{2}(A_{0})$ in Fig. 7(d) a $B_{2}$ that remains essentially constant for shaking amplitudes below a value $A_{0}^{c}$ followed by a linear increase for $A_{0}>A_{0}^{c}$,

Overall, we model the radial dimer interaction by means of the simple function,

[TABLE]

For a system at packing fraction $\phi=0.75$ we obtain $A_{0}^{c}(\phi=0.75)=\left(0.6\pm 0.15\right)\sigma_{g}$ as the shaking amplitude above which the interaction increases linearly with $A_{0}$, $B_{0}(\phi=0.75)=0.3$, and $F=-67\left\langle t_{k}\right\rangle/\sigma_{g}$, where we have prescribed an a priori dependence on the packing density on each of the parameters.

IV.3 Mechanical formation energies

We benefit from the previous section, and in particular from the radial force displayed in Fig. 6 to define the mechanical formation energy of a dimer aligned an angle $\alpha$ with respect to the imposed external driving, $\mathcal{U}_{\alpha}(d)$. Specifically, we define the dimer formation energy as the total mechanical work to bring a pair of inclusions from infinity to a relative distance $d$ along a line at constant angle,

[TABLE]

In equilibrium the energy of formation is derived from the probability to measure the pair at a given configuration $(d,~{}\alpha)$ by the simple expression $\mathcal{U}=-k_{B}T\ln{\mathcal{P}}$. This relation implies, in particular, that the energy of formation scales with the system temperature, $k_{B}T$. Since the granular bed is subject to an energy flux, induced by the tray shaking, we cannot take the previous equality for granted. We subsequently test the validity of this previous relation through an effective bath temperatureSong, Wang, and Makse (2005), $k_{B}T_{\mbox{eff}}$ that could connect the formation energy to the probability to measure a dimer configuration.

To this end, we fit the formation energy of the dimer to the probability to measure the pair of moving inclusions in that dimer configuration. To extract an effective $k_{B}T_{\mbox{eff}}$ for each set of formation energies $\{\mathcal{U}_{\alpha}(d_{i})\}$ and probabilities $\{\mathcal{P}(d_{i},\alpha)\}$ defined at distances $\{d_{i}\}$ we minimize the chi squared function $\chi^{2}=\sum_{i}\left[\mathcal{U}_{\alpha}(d_{i})-\left(k_{B}T_{\mbox{eff}}\ln\mathcal{P}(d_{i},\alpha)+U_{0}\right)\right]^{2}$, with $k_{B}T_{\mbox{eff}}$, and $U_{0}$ the fitting parameters.

Fig. 8 compares the formation energy derived from the relative force and from the angular constrained probabilities with the best fit for $k_{B}T_{\mbox{eff}}$. The corresponding temperatures depend in particular on $\alpha$, which indicates it is impossible to describe the dimer behavior as if it was immersed in an effective equilibrium bath, see Table 1. Specifically, the obtained formation energies display qualitatively different dependence on the dimer distance when compared to $-k_{B}T_{\mbox{eff}}\ln P(d)$. For perpendicular configurations the formation energy of the dimer presents a maximum with $\mathcal{U}>0$ that decays to zero. The minus logarithm of the probability, instead, presents an absolute minimum and then increases asymptomatically to zero. These qualitative differences do not allow to describe this system in terms of a $T_{\mbox{eff}}$.

Alternatively, we can extract an effective temperature from the angular averaged probability, $P(d)$, see Fig. 4. To this end, we introduce the angular averaged mechanical formation energy, $\left\langle\mathcal{U}(d)\right\rangle_{\alpha}$, obtained integrating $\left\langle F_{r}\right\rangle_{\alpha}\equiv 2/\pi\int_{0}^{\pi/2}d\alpha~{}F_{r}(d,\alpha)$. Fig. 9 shows that the formation energy obtained by the two routes shows semiquantitative agreement, and also, by construction, we can assign a $T_{\mbox{eff}}$ for each system. Hence, even if the comparison may be misleading, it is possible to interpret the angular average quantities in terms of an effective equilibrium.

If the functional proposal holds, see Eq. 7), the angular averaged formation energy coincides with the one obtained from the general situation at $\alpha=\pi/4$, and that $-\ln P(d)\sim d^{-1}$ asymptotically. This functional dependence allows a fit between $-\ln P(d)$, and $\left\langle\mathcal{U}\right\rangle_{\alpha}$ at long distances.

IV.4 Measure profiles in the granular bed

We analyze the impact that the inclusions have in the granular bed. Fixing the inclusions allows us not only to measure the relative forces and torques between them, but also quantify the spatial dependence of the physical properties of the grains. We focus on the relative variations of the granular properties along the axis that joins the two dimers. To measure this profile, we consider a rectangular box of size $L\times\sigma_{g}$ that encloses the inclusions in the axis of the dimer. We subsequently slice this box in cells of width $\Delta x=0.2\sigma$ and height $\Delta y=0.2\sigma$ and use them to build the corresponding histograms.

We compute the local density, $\phi(x^{\prime},y^{\prime})$, and kinetic energies $e_{k}(x^{\prime},y^{\prime})$, and $t_{k}(x^{\prime},y^{\prime})$ at each cell with center at $(x^{\prime},y^{\prime})$ and average over oscillation cycles for different initial conditions of the sea of grains. The result is then integrated in the vertical dimension to obtain the dependence on the distance from the surface of an inclusion $h$. The pair of inclusions define an internal region which is identified as $-d/2<h<0$, and an external region $h>0$. To gain statistics, we combine the results from the profiles centered on each inclusion. For visualization purposes we plot the negative region from surface to surface ($-d<h<0$).

[TABLE]

In Fig. 10 we display the relative excess of density, and kinetic energies for $\phi=0.75$, $A_{0}=1.5\sigma_{g}$ relative to the bulk values (9). We compute the profiles for parallel and perpendicular configurations of the dimer and two different separation distances, $d=2\sigma$ and $d=6\sigma_{g}$. For positive values of $h$, the behavior of each profile does not depend on the dimer separation $d$. To confirm this dependence we have computed the profile for a sole inclusion in the system and obtained the same external profiles.

We do not observe a significant deviation of the averaged density along the axis. Moreover, in the region between inclusions, the local density equals to the averaged density outside –except for the density at contact that decays due to the sampling size and the curvature of the disks. However, both kinetic energies –absolute, and relative– are not kept constant along $h$. On the one hand, the excess relative kinetic energy $\delta t_{k}$ has a positive value at contact and relaxes to zero. On the other hand, the excess of absolute kinetic energy departs from negative value and relaxes to zero. For orientations of the dimer perpendicular to the external shaking, $\delta e_{k}$ jumps to a positive value and follows the decay of $\delta t_{k}>0$ and relaxes to 0, whereas for parallel configurations $\delta e_{k}$ does not follow the relative kinetic energy but slowly increases to zero keeping its negative sign.

Finally, we have analyzed the additivity of one-inclusion profiles for the inner region of the dimer, $h<0$. From the disturbance profile generated by a single particle $\delta p(h)$, as presented in (9), we compute the total disturbance profile in the inner region between independent inclusions located at $h=0$, and $h=-d$ as $\delta p(-h)+\delta p(d-h)$. We show the additivity of the profiles by comparing the hollow squares, and dots to the solid lines in Fig. 10, resulting from the addition of one-inclusion profiles and the computation of the profiles for two-inclusion systems, respectively.

V Conclusions

In summary, we have introduced a mechanical model to capture the effective behavior of inclusion in a horizontally shaken granular bed. The analysis of the behavior of an inclusion pair has shown they interact with a granular mediated non-central, and long range interactions. We have considered both a freely moving pair and have computed the forces for a pair at a fixed distance. These two approaches have provided complementary insight.

We have extracted the two dimensional probability density of the relative vector between inclusions and measured a preference of dimers at contact to align parallel to the external forcing. Consistently, the force measurements with fixed inclusions have revealed the emergence of an attractive force with a maximum for particles at contact plus the emergence of an aligning torque at short distances. The detailed computation of the forces between the two inclusions as a function of the angle they form with the driving force indicates that inclusions tend to remove the ordered structure of the granular bed for $\alpha<\pi/2$ and moderate amplitudes of the shaking $A_{0}>\sigma_{g}/2$. When the granular structure around the inclusions becomes weak in the internal region of the pair, the long range attractive force between the pair becomes apparent. We have characterized the interaction in terms of the parameters of the dimer configuration and shaking of the granular bed with a simple model of interaction at a distance (7). We have matched the anisotropy in the probability $\mathcal{P}(r_{x},r_{y})$ to the appearance of a neat torque in for fixed inclusions. The torque aligns close inclusions $d\lesssim 2\sigma_{g}$ parallel to the shaking whereas at $d\gtrsim 4\sigma_{g}$ and large shaking amplitudes and granular densities an emergent torque aligns inclusions perpendicularly to the shaking. The appearance of this dealigning torque could be the mechanism leading towards the segregation of mixtures of grains and inclusions into bands perpendicular to the shaking directionMullin (2000); Reis and Mullin (2002); Reis et al. (2004); Ciamarra, Coniglio, and Nicodemi (2006). Finally, we have captured the kinetic energy profiles in the granular bed in the immediacies of the inclusion dimer and observed a differentiated behavior of the absolute kinetic energy on the orientation of the dimer. Additionally we observe that the perturbation of the granular kinetic energies in the internal region of the dimer results from the addition of two independent facing inclusion surfaces

Acknowledgements.

The authors acknowledge I. Zuriguel for a detailed description of the experimental setup and insightful discussions. This work has been supported by MINECO and DURSI projects PGC2018-098373-B-100 and 2017 SGR 884, respectively. JC was supported by Spanish MECD fellowship FPU13/01911. IP also acknowledges SNF Project No. 200021-175719 for funding.

Appendix A Stochastic Dynamics of Grains

We describe the dynamics of granular particles at positions $\bm{r}$ and momenta $\bm{p}$

[TABLE]

where $m$ is the mass of each granular particle. The dynamics are given by the following Langevin set of equations

[TABLE]

where $\gamma$ are the dissipation constants, $\bm{v}_{s}$ the velocity of the substrate, and $\bm{F}$, and $\bm{f}$ the deterministic and random forces, respectively.

To solve equation (11) we separate the terms that depend on the momenta $\bm{p}$ on LHD separate the deterministic $\bm{\Phi}=\bm{F}+\gamma\bm{v}_{s}$, and random forces $\bm{f}_{i}$ on the RHS

[TABLE]

The solution to equation (12) is computed as $\bm{p}=\bm{p}^{(h)}+\bm{p}^{(p)}$ where homogeneous solution is,

[TABLE]

and a particular solution, $\bm{p}^{(p)}(t+t_{0})$, of the form, $\bm{\pi}(t)\exp\left(-\gamma/mt\right)$, where $\bm{\pi}(t)$ satisfies the following random process,

[TABLE]

and the solution

[TABLE]

To compress the result we define a new set of constants. The dissipation of moment is introduced by $\Gamma_{i}=\exp\left(-\gamma_{i}\Delta t/m_{i}\right)$; the mobility coefficient for the deterministic force $\tilde{\mu}_{i}=\left(1-\Gamma\right)/\gamma_{i}$; the effective diffusion coefficient $\tilde{D}_{i}=\left(\Gamma_{i}m_{i}/\gamma_{i}\right)^{1/2}\left(1-\Gamma_{i}^{2}\right)^{1/2}$. Now the evolution of $\bm{p}$ is simplified to,

[TABLE]

with a Gaussian random noise with zero mean and variance $\left\langle\zeta_{i,\alpha}(t)\zeta_{j,\beta}(t^{\prime}t)\right\rangle=2\Lambda_{i,\alpha}\delta_{ij}\delta_{\alpha\beta}\delta(t-t^{\prime})$. The evolution of the grains is simplified to

[TABLE]

where we first integrate equation (10) with $\bm{p}(t)$ to $t+\Delta t/2$, and second we integrate the positions $\bm{r}(t+\Delta t/2)$ to $t+\Delta t$ with momenta given by $(\bm{p}(t)+\bm{p}(t+\Delta t)/2$ from eq. (16).

Appendix B Granular measures

In order to characterize the granular bed we introduce the measure of kinetic energies. For a system under shaking we can construct either an absolute, and a relative kinetic energy,

[TABLE]

with $\bm{V}_{cm}$ the centre of mass velocity. The usual magnitude Song, Wang, and Makse (2005) defined in granular systems is $t_{k}$ but the introduction of external bodies gives a certain relevance to $e_{k}$. Here we discuss the advantages of using each of them to scale the interaction forces obtained in Section IV.

The relative kinetic energy $t_{k}$ does not depend on the shaking amplitude $A_{0}$ and it is best suited to compare forces in systems at different shaking amplitudes and granular properties such as MSD.

The absolute kinetic energy $e_{k}$ includes the flux of energy that the shaking introduces to the grains, and it is relevant for the interactions with external bodies. Even though both energies are suitable to define an energy scale we choose the averaged $t_{k}$ for the sake of commensurability.

In Fig. 11 we present the dependence of the averaged energies measured at different densities, forcing amplitudes, and $\Lambda_{x}$, for a system with fixed granular parameters: $\Lambda_{x}/\Lambda_{y}=2$, $k=2.5\ 10^{4}$,  $\gamma_{s}=20$, and $\gamma=33$ in system units.

We observe a decay of the kinetic energy as the density of grains increases. Dissipation in particle-particle collision events in granular systems introduces an energy loss and it is responsible for the energy dissipation of grains as density increases. The difference in kinetic energies is solely due to the external shaking and it is confirmed in Fig. 11(b) where the the kinetic energy difference remains approximately constant in the considered packing fraction range. We determine a monolayer to be dense when the energy loss is relevant, in this case for $\phi>0.6$. Since we are interested in a dense system we determine our lower bound for the packing range to be $\phi=0.6$

The upper bound for the packing fraction must correspond to a system within the same mobility regime. For this purpose we measure the jump distribution of grains, $P(\Delta)$ as the probability density distribution of the displacement $\Delta$ of grains after a time unit.

Fig. 12 shows a Gaussian jump distributions for $\phi\leq 0.75$ at $\Lambda_{x}=50$. For this reason we set the packing upped limit for exploration in this paper at $\phi=0.75$.

At larger densities, for example, $\phi=0.8$ the jump distribution presents a pronounced peak at $|\Delta|\approx 0$ that indicates the caging effect of the neighboring grains and it is no longer Gaussian. Mean Squared Displacements, not shown here, also confirm the caging of grains, and thus the system is not in a dense fluidized regime.