# Collective drifts in vibrated granular packings:
the interplay of friction and structure.

A. Plati  Department of Physics, University of Rome Sapienza, P.le Aldo Moro
2, 00185, Rome, Italy  Institute for Complex Systems - CNR, P.le Aldo Moro 2,
00185, Rome, Italy  A. Puglisi  Department of Physics, University of Rome
Sapienza, P.le Aldo Moro 2, 00185, Rome, Italy  Institute for Complex Systems
- CNR, P.le Aldo Moro 2, 00185, Rome, Italy  INFN, University of Rome Tor
Vergata, Via della Ricerca Scientifica 1, 00133, Rome, Italy

###### Abstract

We simulate vertically shaken dense granular packings with horizontal periodic
boundary conditions. A coordinated translating motion of the whole medium
emerges when the horizontal symmetry is broken by disorder or defects in the
packing and the shaking is weak enough to conserve the structure. We argue
that such a drift originates in the interplay between structural symmetry
breaking and frictional forces transmitted by the vibrating plate. A non-
linear ratchet model with stick-slips reproduces many faces of the phenomenon.
The collective motion discussed here underlies phenomena observed recently
with vibrofluidized granular materials, such as persistent rotations and
anomalous diffusion.

_Introduction._ — One of the major challenges of statistical mechanics,
nowadays, is understanding systems far from equilibrium, e.g. in the presence
of energy flows or dissipation [1]. A realm of physics where fluctuations and
lack of thermodynamic equilibrium are ubiquitous is soft matter [2] and, in
this context, an established test-ground for theories is provided by
vibrofluidized granular systems [3, 4, 5, 6, 7] where energy is continuously
injected from an external source and dissipated in friction. When isotropy is
broken, a fraction of this energy current can be exploited to realise a
ratchet effect [8, 9, 10, 11, 12, 13] similar to what is seen in active matter
[14]. A series of recent granular experiments suggest that spontaneous
persistent drifts emerge in the presence of isotropic disorder: this is
revealed by superdiffusion [15], which has been connected to transient drifts
with very long relaxation times [16, 17], or by the appearance of steady
rotations of disks in dense vibrated packings [18] and by the emergence of
collective motions under swirling excitations [19, 20]. In all these examples
random energy is converted into a steady flow, realising an interesting class
of “disordered engines”.

Here we aim at understanding the general ingredients underlying these effects,
focusing on a few simplified setups and - in the conclusions - on a model of
frictional ratchet. The mechanism investigated explains the aforementioned
phenomena as it represents their translational (instead of rotational) [15,
16, 17], collective (instead of individual) [18] and spontaneous (instead of
externally stimulated) [19, 20] counterpart.

_Numerical setup._ — We simulate, through an established Discrete Elements
Method (DEM), a system of $N$ spherical grains confined by hard walls and/or
periodic boundary conditions (PBC) with different geometries. Each grain has
radius $R_{i}$, mass $m_{i}$, position $r_{\alpha i}$, velocity $v_{\alpha i}$
and angular velocity $\omega_{\alpha i}$, where $\alpha=\\{x,y,z\\}$ and
$i=\\{1,\cdots,N\\}$. The grain-grain and grain-boundary contact forces follow
the Hertz Mindlin (HM) model, see details in the Supplemental Material (SM)
[21] . In all cases the setup is enclosed in a 3D box of height $10\times
2\max_{i}R_{i}$, energy is injected by vertical vibrations $z(t)=A\cos(2\pi
ft)$ of the upper and lower confining hard walls, and there is one or more
horizontal direction with infinite horizon (i.e. where movement is not
constrained by walls). In quasi-2D setups (cases in Fig. 1a-g, $N=60$) $x$ has
PBC while $y$ is confined by two parallel vertical walls of width $L=32$ mm
separated by a distance $d=2\max_{i}R_{i}$. In the full 3D case in Fig. 1h
($N=2600$) both $x$ and $y$ directions have PBC and the base of the box has
dimensions $98\times 98$ mm2. In the full 3D case in Fig. 1i ($N=2600$) the
box is a cylinder with a conical-shaped base, identical to what previously
used in experiments and simulations [15, 16, 17]. The shaking intensity is
measured by $\Gamma=A(2\pi f)^{2}/g$ ,where $g$ is the gravity acceleration.
Here we always vary $\Gamma$ by $A$ keeping $f$ constant at 100 Hz in the
quasi-2D simulation and at 200 Hz in the all other ones.

Figure 1: Simulation geometries. Quasi 2D vertical layers can be random
polydisperse (a) – with equal proportions of three species of radii
$R_{i}=\\{1.5,2.0,2.5\\}$ mm – or ordered monodisperse (b-g) – $R_{i}=R=2$ mm
– with eventually defects. The defect nomenclature is related to the position
of the grain that breaks the symmetry. The 3D setups are cubic (h) and cone-
base cylinder (i) (see [16] for details), both monodisperse.

Regarding the quasi-2D setups, we prepare the initial state of our packings in
two main ways: the first one consists of randomly pouring a polydisperse
assembly of grains in the container, initial velocities are zero but during
the pouring dynamics they acquire energy and rapidly reach a stationary
statistics; the second one is obtained by placing a monodisperse assembly on
an hexagonal lattice with the possibility to have vacancies in determined
sites, while the particles initial velocities are drawn from a Gaussian
distribution with zero mean and a variance small enough to keep the crystal
stable. For the 3D cases we simply pour monodisperse grains in the containers.
Since the latter do not satisfy the right proportions for crystallization the
resulting packings are fairly disordered. The monodisperse packings with
quasi-2D geometry (b-g in Fig. 1) can be without defects (g), or with four
symmetric defects (f), or with three defects placed to break in different ways
the symmetry of the crystal with respect to $z$ and $x$ (b-e). The
polydisperse quasi-2D packing (Fig. 1a) is inspired by a recent experimental
and numerical study [18] with shaken disks without PBC: for $\Gamma$ low
enough a persistent net angular velocity is seen for each disk. We verified
(see SM [21]) that such persistent rotational modes are present also in our
simulations with spheres, in all the geometries explored here, with both hard
walls and PBC. However, in this Letter, we focus on simulations with
horizontal PBC only, where a different phenomenon, namely the collective
horizontal persistent drift, superimposed on the previously observed
rotational modes, appears.

_Drifting disordered packings -_ We start from polydisperse packings in the
quasi-2D geometry. At low values of $\Gamma$ all the system moves coherently
with the center of mass (CM) along the $x$-axis i.e. the free direction
allowed by PBC, whose motion in time
$X^{CM}(t)=M_{\text{tot}}^{-1}\sum_{i}m_{i}x_{i}(t)$ is shown in Fig. 2a. We
note that, even with the same $\Gamma$, a different random packing can lead to
a stable drift with very different magnitudes or to an intermittent drift. In
panel b, considering the $x$ component of the CM’s velocity $V_{x}^{CM}(t)$,
we also verify that short time properties as the time variance
$\sigma^{2}\left(V_{x}^{CM}\right)$ are fully determined by the driving
parameters (small errorbars), while the slow cooperative dynamics
(characterized by its time average $\langle V_{x}^{CM}\rangle$), sensibly
depends on the packing configuration (large errorbars). Nevertheless, raising
$\Gamma$ up to values for which the system fluidizes (typically $\Gamma>10$),
the CM performs Brownian-like trajectories and $\langle V_{x}^{CM}\rangle$
vanishes (SM-Video1 [21]). We recall that a dense vibrofluidized granular
system has several time-scales, the smallest associated to fast vibrational
motion, the largest associated to slow rearrangements of the global contact
network. What we observe here is a rapid divergence of the largest time-scale
when $\Gamma$ is reduced below $\sim 3\div 10$. The value of $\langle
V_{x}^{CM}\rangle$, therefore, does not depend crucially on the trajectory’s
duration, provided it is longer than the small time-scales, e.g. $\gg 10^{-1}$
s.

Figure 2: a) $X^{CM}(t)$ (symbols) and $x_{i}(t)$ for all the grains $i$ in
the system (lines) for quasi-2D random packings: the grains move coherently
with the CM. Simulations with the same $\Gamma$ refer to different random
realizations. b) Comparison between $\sigma^{2}(V_{x}^{CM})$ and $\langle
V_{x}^{CM}\rangle$, obtained averaging over three independent random
realizations. c) Time averaged velocity of center of mass as a function of
$\Gamma$ for different monodispersed packings. Each point is mediated over
five independent realizations of the dynamics. Here, differently from random
packings, $\langle V_{x}^{CM}\rangle$ vanishes for $\Gamma\simeq 3$ because
moderate vibrations destroy the asymmetric configuration of defects. d) zoom
on low $\Gamma$s to highlight the sign changes with fixed structure.

_Effect of symmetries -_ The observation that average motion is erased by fast
particle rearrangements (fluidization) is a hint of its correlation with the
system’s spatial configuration. For this reason we study $\langle
V_{x}^{CM}\rangle$ for ordered packings with defects (Fig. 2c). Net drifts are
never observed in the FC and SD cases, whereas for all asymmetric
configurations they are. Moreover, when the defect configuration is mirrored
with respect to the $z$-axis, e.g. BL $\to$ BR or TR $\to$ TL etc., at a given
value of $\Gamma$, $\langle V_{x}^{CM}\rangle$ changes its sign remaining with
comparable magnitude. A reflection of the configuration with respect to $x$
always changes the magnitude of the drift but the sign reverses just for a few
values of $\Gamma$. Remarkably, a variation of $\Gamma$ with keeping the same
layout of defects brings to multiple inversions of $\text{sgn}(\langle
V_{x}^{CM}\rangle)$ in the region $0.8\leq\Gamma\leq 1.05$ zoomed in Fig. 2d.
In summary: i) defect asymmetry is needed to observe non-zero mean velocity of
the CM, ii) reflection w.r.t. $z$ inverts the direction of motion, iii) the
defect configuration alone does not define the verse of the motion, the
driving parameter determines it too (SM-Video2 [21]).

Figure 3: a-d) Comparison between the time evolution of the CM velocity and
the mean coordination number of the contact network
$Z(t)=N^{-1}\sum_{ij}\Theta(R_{i}+R_{j}+\delta-|\boldsymbol{r}_{ij}|)$ where
$0<\delta\ll R_{i}+R_{j}$ allows to detect nearest neighbour not in contact
due to vibrations. Trajectories are smoothed with a running average. a-b)
random packings (cases with upside-down triangles and circles in Fig. 2a). c)
3D cubic setup, here $V_{xy}^{CM}$ is two dimensional so we plot both the
modulus and the orientation $\theta_{xy}^{CM}$. d) 3D cylindrical geometry,
the collective motion is a rotation around the central axis so
$\Omega_{c}=N^{-1}\sum_{i}|\boldsymbol{r}_{i}|^{-2}(\boldsymbol{r}_{i}\times\boldsymbol{v}_{i})_{z}$.
e) Time averaged modulus of the normal force and compenetration between the
plate and the bottom particles for two $z$-reflected configurations of
defects. Site 4 is under the double defect (see Fig 1b-c). f-h) Scatter plots
with bins of the total $x$-tangential force modulus exerted by the plate VS
the total velocity of the bottom particles for ordered packings with defects
(f) and random packings (g-h) shaken at $\Gamma=0.7$. In the latter,
$\gamma^{\pm}=|\langle F^{T}_{x}|v_{x}\gtrless 0\rangle|/\langle
v_{x}|v_{x}\gtrless 0\rangle$ linearly interpolate the data clouds.

_Sensitivity to small configuration changes -_ Given the observed complicate
entanglement of external parameters and drifts, we now focus on fluctuations
of the drift during the same trajectory. Such fluctuations are more evident
and frequent in random packings: Fig. 3a-b compares the coordination number of
the contact network $Z(t)$ with the instantaneous drift $V_{x}^{CM}(t)$: they
often change together. Such changes do not correspond to a significant
rearrangement of particles (changes of $Z(t)$ can be smaller than $1\%$),
confirming that the collective motion is sensible to small deformations of the
contact network. Collective drifts are observed also in the 3D realistic
setups (see [16, 17] and SM [21]) and similar correlations are present here
too (Fig. 3c-d). Obviously, the coordination number is not directly related to
the packing asymmetry. A more appropriate quantity eluded our analysis,
perhaps because of the observed challenging scenario: on one hand, the drift
may significantly change at fixed $\Gamma$ for weak deformations of the
contact network (as seen for random packings) and, on the other hand, it can
change direction by varying the driving parameters with a fixed structure (as
occurs in the ordered packings with defects). Analogous difficulties have been
discussed recently for glassy systems [35, 36, 37, 38, 39, 40]. In our study,
the athermal nature of the system (that needs a mechanical driving to reach a
stationary state) is a peculiar feature that has no counterpart in thermal
glasses [41].

_Bulk asymmetries originate a ratchet effect -_ We now look for observables
that mediate between bulk structure and dynamics. The profiles, along $x$ of
the mean normal pressure and the mean compenetration between the plate and the
bottom particles are shown in Fig. 3e, for an ordered packing with defects.
Each profile is asymmetric and properly inverted under a $z$-reflection:
defects in the bulk actually influence the way in which the boundaries
interact with the source of energy. This may be expected, since both normal
pressure and compenetration affect the tangential component of the force
exerted by the vibrating plate (see SM [21]): what is remarkable in our
opinion is its global dynamical consequence. Recent analytical results [42]
for an idealized model of vibrofluidized dense granular material suggest long-
range spatial correlations that may justify the global influence of defects. A
deeper insight is provided by the external tangential force provided by the
plate $F^{T}_{x}$ ($x$-component), which directly affects the CM:
$\dot{V}^{\text{CM}}_{x}=M_{\text{tot}}^{-1}F^{T}_{x}(r_{xj},r_{zj},v_{xj},\omega_{yj})$
where the index $j$ refers to the bottom particles (for ordered packings
$j\in[1,9]$). We characterise the dependence of $F^{T}_{x}$ upon
$v_{x}=\sum_{j}v_{xj}$ by the scatter plot in Fig. 3f for two $z$-reflected
ordered configurations with defects. Such dependence is typical of a
frictional force but, quite surprisingly, its intensity depends on the sign of
the velocity. Moreover, also this frictional asymmetry is inverted under
$z$-reflection of packing. Comparing with Fig. 2d clarifies that the drift
occurs towards lower friction. Panels g-h show that the same kind of bias is
more pronounced for random packings, indeed $|\langle V_{x}^{CM}\rangle|$ is
far larger. All this makes clear that structural disorder and defects
introduce an asymmetric interaction with the plate. The role of the
$\omega_{yj}$s is discussed in the SM [21].

Figure 4: a) Time evolution of $F_{c}=\mu|F^{N}|-|F^{T}_{x}|$. Diamonds mark
the instants for which $F_{c}\sim 0$ corresponding to the sliding condition
for the tangential force between the plate and the bottom grains. b)
Trajectory of the $X^{CM}$ compared with the horizontal displacement of the
bottom layer accumulated in the sliding instants ($\Delta X^{b}_{\text{sl}}$)
and in all the other ones ($\Delta X^{b}_{\text{ot}}$). c) Stick-slip in the
difference of the mean $x$-coordinates of the two lowest horizontal layers in
a granular simulation. Diamonds refer to sliding instants. These three panels
are obtained with BL and $\Gamma=1.25$. d) Mean velocity of the CM as a
function of the horizontal (triangles) and vertical (diamonds) size of the
system for different values of the defect concentration. The vertical size
effect is shown for $c=1$ but we verified that its relevant behavior does not
depend on $c>0$. e) Scaling of $\langle V_{x}^{CM}\rangle$ and the diffusivity
of the fluctuation around the drift $D$ (see SM [21]) as a function of the
horizontal system size for random packings. Each point is obtained averaging
over ten independent realizations of the packing, error bars are standard
deviations.

Simulations with a higher acquisition rate allow to study the grain
trajectories on timescales shorter than a driving period $1/f$. In Fig. 9a, we
show $F_{c}(t)=\mu|F^{N}|-|F^{T}_{x}|$ where $\mu$ is the dynamic friction
coefficient. As explained in the SM [21], this difference is zero if the grain
and the plate surfaces slide on each other. We observe roughly periodical
cycles (with the periodicity of the driving plate), with a part of each period
where $F_{c}\sim 0$. In Fig. 9b, we compare the trajectory of the CM with the
total horizontal displacement accumulated by the bottom layer considering just
the sliding instants and all the remaining ones. From this comparison we
conclude that the main contribution to the irreversible motion occurs when
$F_{c}\sim 0$. From Fig. 9c we also understand how the mean horizontal motion
originated at the bottom of the system is transfered to the top: there we plot
the difference $\Delta X$ between the mean $x-$coordinate of the lowest layer
and of the one just above it, noting a clear stick-slip dynamics. We have
verified that, as suggested by the marked points in the graph, this distance
slowly increases during the sliding of the bottom layer and then suddenly
decreases when the sliding condition ceases.

The provided analysis suggest the interpretation of the phenomenon under study
as a ratchet effect originated from structural asymmetries [43, 44]. To
further support this, a variation of the well known periodically rocked
ratchet [45] is proposed in the SM [21]. Inspired by the stick-slip dynamics
between the layers of the granular packing we have replaced the asymmetric
periodic potential with a confining but slipping one. Our model reproduces
many features of the phenomenon under study (see SM [21]).

_System’s size effects-_ The results for 2D packings are obtained in a
relatively small system in order to establish a direct connection with former
experiments [18]. The presence of analogous behavior also in larger and more
realistic 3D systems signals the generality of the phenomenon. Nevertheless, a
study of the system’s size effects in the simplified geometry is non-trivial
since the origin of the drift and the resistance against it can depend both on
bulk volume (through the concentration of defects) and external surfaces
(through friction and energy input). We performed simulation of ordered
packings with defects where the original ”module” of $N\sim 60$ grains is
replicated $\mathcal{N}$ times along both the $x$ and the $z$ axis. Each
replica can be with or without defects so that a concentration of defects $c$
can be defined as the number of asymmetric modules over $\mathcal{N}$. For the
horizontal size-scaling we considered three different values of defect
concentration $c=\\{1/\mathcal{N},0.5,1\\}$. From Fig. 9d we see that
increasing the system size along $x$ with a vanishing concentration of defects
($c=1/\mathcal{N}$) leads to the weakening of the drift. For finite defect
concentration instead, we observe a non zero-asymptotic $\langle
V_{x}^{CM}\rangle$ whose magnitude depends on $c$. Regarding the vertical
size-scaling, we note that the drift is suppressed for large $\mathcal{N}$.
This fact can be reasonably explained considering that, in absence of lateral
wall in the $x$-direction, the pressure at the base of the system increases
with its height reducing the global mobility. We have accordingly verified
that for packings higher than $\mathcal{N}=5$, the sliding condition
$\mu|F^{N}|\sim|F^{T}_{x}|$ is never satisfied. The scenario for large random
packings is different: the global asymmetry of many disordered granular
patches each pushing the whole system in a different direction is expected to
decrease with their number. Indeed we observe a reduction of the drift as
$1/\sqrt{\mathcal{N}}$ when the system size is scaled horizontally (Fig. 9e).
Interestingly, also the diffusivity of the fluctuations around the average
motion decreases with $\mathcal{N}$. We define such a quantity as the
diffusion coefficient $D$ of the fast component of the CM dynamics that can be
measured from the power spectral density of $V^{CM}_{x}(t)$ as shown in the SM
[21]. For $\Gamma>1$ it follows a scaling $D\sim 1/\mathcal{N}$ while for
$\Gamma<1$ it seems to reach a minimum at large $\mathcal{N}$. In the former
case the quantity $D/\langle V_{x}^{CM}\rangle^{2}$ that is is the typical
time after which the drift (even very small) becomes visible is independent
from $\mathcal{N}$. This observation confirms that the drift is more relevant
than a finite-size effect and poses an interesting challenge for the recently
discovered Thermodynamic Uncertainty Relations [46, 47].

_Conclusion._ — We have presented a numerical study of a realistic granular
contact model with several ideal experiments to put in evidence the existence
of a random-to-direct energy conversion based upon the concurrent breaking of
time and space symmetries. As in many ratchet-like phenomena, the same
geometry may lead to opposite drifts, depending on energy injection. Inspired
by the interplay of friction and structure emerged in the numerical analysis,
we introduce a novel ratchet model with asymmetric interactions and stick-
slips. The drifts discovered and explained underlie several phenomena observed
recently in shaken granular media [15, 16, 17, 18] and are expected to be a
general feature of soft matter systems, such as active matter, crowd dynamics
etc.

###### Acknowledgements.

The authors are indebted to Marco Baldovin and Rafael Díaz Hernández Rojas for
fruitful scientific discussions. The authors acknowledge the financial support
from the MIUR PRIN2017 project 201798CZLJ

## References

  * Marconi  _et al._ [2008] U. M. B. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani, Phys. Rep. 461, 111 (2008).
  * Nagel [2017a] S. R. Nagel, Rev. Mod. Phys. 89, 025002 (2017a).
  * Kadanoff [1999] L. P. Kadanoff, Rev. Mod. Phys. 71, 435 (1999).
  * Jaeger  _et al._ [1996] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Rev. Mod. Phys. 68, 1259 (1996).
  * Rodhes [1997] M. Rodhes,  _Principles of powder technology_ (John Wiley and Sons, 1997).
  * Nagel [2017b] S. R. Nagel, Rev. Mod. Phys. 89, 025002 (2017b).
  * de Arcangelis  _et al._ [2019] L. de Arcangelis, E. Lippiello, M. Pica Ciamarra, and A. Sarracino, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 377, 20170389 (2019).
  * Reimann [2002] P. Reimann, Physics Reports 361, 57 (2002).
  * Hänggi and Marchesoni [2009] P. Hänggi and F. Marchesoni, Rev. Mod. Phys. 81, 387 (2009).
  * Costantini  _et al._ [2007] G. Costantini, U. Marini Bettolo Marconi, and A. Puglisi, Phys. Rev. E 75, 061124 (2007).
  * Eshuis  _et al._ [2010] P. Eshuis, K. van der Weele, D. Lohse, and D. van der Meer, Phys. Rev. Lett. 104, 248001 (2010).
  * Gnoli  _et al._ [2013] A. Gnoli, A. Petri, F. Dalton, G. Pontuale, G. Gradenigo, A. Sarracino, and A. Puglisi, Phys. Rev. Lett. 110, 120601 (2013).
  * Balzan  _et al._ [2011] R. Balzan, F. Dalton, V. Loreto, A. Petri, and G. Pontuale, Phys. Rev. E 83, 031310 (2011).
  * Di Leonardo  _et al._ [2010] R. Di Leonardo, L. Angelani, D. Dell’Arciprete, G. Ruocco, V. Iebba, S. Schippa, M. P. Conte, F. Mecarini, F. De Angelis, and E. Di Fabrizio, Proceedings of the National Academy of Sciences 107, 9541 (2010), https://www.pnas.org/content/107/21/9541.full.pdf .
  * Scalliet  _et al._ [2015] C. Scalliet, A. Gnoli, A. Puglisi, and A. Vulpiani, Phys. Rev. Lett. 114, 198001 (2015).
  * Plati  _et al._ [2019] A. Plati, A. Baldassarri, A. Gnoli, G. Gradenigo, and A. Puglisi, Phys. Rev. Lett. 123, 038002 (2019).
  * Plati and Puglisi [2020] A. Plati and A. Puglisi, Phys. Rev. E 102, 012908 (2020).
  * Moukarzel  _et al._ [2020] C. F. Moukarzel, G. Peraza-Mues, and O. Carvente, Phys. Rev. Lett. 125, 028001 (2020).
  * Zhao and Pöschel [2021] S.-C. Zhao and T. Pöschel, Physics of Fluids 33, 081701 (2021), https://doi.org/10.1063/5.0056143 .
  * Zhao and Pöschel [2022] S.-C. Zhao and T. Pöschel, Phys. Rev. E 105, L022902 (2022).
  * [21] See Supplemental Material http://link.aps.org/supplemental/XXX for more informations and additional numerical analysys, which includes Refs. [22-34].
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics (1995).
  * Lam [a] http://lammps.sandia.gov/ (a).
  * Zhang and Makse [2005] H. P. Zhang and H. A. Makse, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 72, 1 (2005).
  * Silbert  _et al._ [2001] L. E. Silbert, D. Ertaş, G. S. Grest, T. C. Halsey, D. Levine, and S. J. Plimpton, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 64, 1 (2001).
  * Brilliantov  _et al._ [1996a] N. V. Brilliantov, F. Spahn, J. M. Hertzsch, and T. Pöschel, Physical Review E 53, 5382 (1996a).
  * Popov [2010] V. L. Popov,  _Contact Mechanics and Friction_ (Springer-Verlag, Berlin, 2010).
  * Lam [b] https://docs.lammps.org/pair_granular.html (b).
  * Di Renzo and Di Maio [2004] A. Di Renzo and F. P. Di Maio, Chemical Engineering Science 59, 525 (2004).
  * T. Pöschel [2005] T. S. T. Pöschel,  _Computational Granular Dynamics_ (Springer, Berlin, 2005).
  * Rackl and Hanley [2017] M. Rackl and K. J. Hanley, Powder Technology 307, 73 (2017).
  * Plati  _et al._ [2021] A. Plati, L. de Arcangelis, A. Gnoli, E. Lippiello, A. Puglisi, and A. Sarracino, Phys. Rev. Research 3, 013011 (2021).
  * Wilkie [2004] J. Wilkie, Phys. Rev. E 70, 017701 (2004).
  * Stukowski [2010] A. Stukowski, Modeling and Simulation in Materials Science and Engineering 18, 10.1088/0965-0393/18/1/015012 (2010).
  * Díaz Hernández Rojas  _et al._ [2021] R. Díaz Hernández Rojas, G. Parisi, and F. Ricci-Tersenghi, Soft Matter 17, 1056 (2021).
  * Coulais  _et al._ [2014] C. Coulais, R. P. Behringer, and O. Dauchot, Soft Matter 10, 1519 (2014).
  * Tong and Tanaka [2018] H. Tong and H. Tanaka, Phys. Rev. X 8, 011041 (2018).
  * Cubuk  _et al._ [2015] E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rottler, D. J. Durian, E. Kaxiras, and A. J. Liu, Phys. Rev. Lett. 114, 108001 (2015).
  * Hentschel  _et al._ [2019] H. G. E. Hentschel, I. Procaccia, and S. Roy, Phys. Rev. E 100, 042902 (2019).
  * Schwartzman-Nowik  _et al._ [2019] Z. Schwartzman-Nowik, E. Lerner, and E. Bouchbinder, Phys. Rev. E 99, 060601(R) (2019).
  * Gradenigo  _et al._ [2010] G. Gradenigo, A. Sarracino, D. Villamaina, T. S. Grigera, and A. Puglisi, J. Stat. Mech. 2010, L12002 (2010).
  * Plati and Puglisi [2021] A. Plati and A. Puglisi, Scientific Reports 11, 14206 (2021).
  * van den Broek  _et al._ [2009] M. van den Broek, R. Eichhorn, and C. V. den Broeck, EPL (Europhysics Letters) 86, 30002 (2009).
  * Kumar  _et al._ [2008] K. V. Kumar, S. Ramaswamy, and M. Rao, Phys. Rev. E 77, 020102(R) (2008).
  * Bartussek  _et al._ [1994] R. Bartussek, P. Hänggi, and J. G. Kissner, Europhysics Letters (EPL) 28, 459 (1994).
  * Barato and Seifert [2015] A. C. Barato and U. Seifert, Phys. Rev. Lett. 114, 158101 (2015).
  * Horowitz and Gingrich [2020] J. M. Horowitz and T. R. Gingrich, Nat. Phys. 16, 15 (2020).
  * Brilliantov  _et al._ [1996b] N. V. Brilliantov, F. Spahn, J. M. Hertzsch, and T. Pöschel, Physical Review E 53, 5382 (1996b).

SM: SUPPLEMENTAL MATERIAL

##  I Model for contact dynamics

We develop our DEM simulations through the LAMMPS package [22, 23] using the
Hertz-Mindlin model [24, 25, 48] to solve the dynamics during the particle-
particle and particle-container contact. This visco-elastic model takes into
account both the elastic and the dissipative response to the mutual
compression between the grains. In addition, it includes in the dynamics not
only the relative translational motion but also the rotational one. The
contact forces are described by Eqs. (1) below.

$\displaystyle{\boldsymbol{F}}^{N}_{ij}=\sqrt{R_{ij}^{\text{eff}}}\sqrt{\xi_{ij}(t)}\left[(k_{n}\xi_{ij}(t)-m^{\text{eff}}_{ij}\gamma_{n}\dot{\xi}_{ij}(t))\cdot{\bf
n}_{ij}(t)\right]$ (1a)
$\displaystyle{\boldsymbol{F}}^{T}_{ij}=\begin{cases}{\boldsymbol{F}}^{\text{hist}}_{ij}\quad\text{if}\quad|{\boldsymbol{F}}^{\text{hist}}_{ij}|\leq|\mu{\boldsymbol{F}}^{N}_{ij}|\\\
-\dfrac{|\mu{\boldsymbol{F}}^{N}_{ij}|}{|{\bf g}^{T}_{ij}(t)|}\cdot{\bf
g}^{T}_{ij}(t)\quad\text{otherwise}\\\ \end{cases}$ (1b)
$\displaystyle{\boldsymbol{F}}^{\text{hist}}_{ij}=-\sqrt{R_{ij}^{\text{eff}}}\left[k_{t}\displaystyle\int\sqrt{\xi_{ij}(t^{\prime})}\textbf{ds}_{ij}(t^{\prime})+m^{\text{eff}}_{ij}\gamma_{t}\sqrt{\xi_{ij}(t)}{\bf
g}^{T}_{ij}(t)\right]$ (1c)

These equations are written for two spherical particles with radius $R_{i}$,
$R_{j}$, mass $m_{i}$,$m_{j}$ position ${\bf r}_{i}$, ${\bf r}_{j}$,
translational velocity ${\bf v}_{i}$, ${\bf v}_{j}$ and rotational velocity
${\boldsymbol{\omega}}_{i}$, ${\boldsymbol{\omega}}_{j}$. The relative
velocity is defined as ${\bf g}_{ij}=(\dot{{\bf
r}}_{i}-{\boldsymbol{\omega}}_{i}\times R_{i}{\bf n})-(\dot{{\bf
r}}_{j}+{\boldsymbol{\omega}}_{j}\times R_{j}{\bf n})$ where ${\bf
n}=\left({\bf r}_{i}-{\bf r}_{j}\right)/\left|{\bf r}_{i}-{\bf r}_{j}\right|$
defines the normal unitary vector of the contact surface; we call ${\bf
g}_{ij}^{N}$ and ${\bf g}_{ij}^{T}$ the two components of ${\bf g}_{ij}$
respectively normal and parallel to ${\bf n}$. The instantaneous overlap
between the two grains is represented by $\xi_{ij}(t)=R_{i}+R_{j}-|{\bf
r}_{i}-{\bf r}_{j}|$ and its derivative is $\dot{\xi}_{ij}(t)=|{\bf
g}^{N}_{ij}|$. Regarding the effective parameters we have
$R_{ij}^{\text{eff}}=R_{i}R_{j}/(R_{i}+R_{j})$ and
$m_{ij}^{\text{eff}}=m_{i},m_{j}/(m_{i}+m_{j})$; here is important to mention
that the wall surfaces are treated as spheres with infinite mass and radius.
During the contact, the particles are subjected to a normal force
${\boldsymbol{F}}^{N}_{ij}$ and a tangential one ${\boldsymbol{F}}^{T}_{ij}$;
both these components have an elastic and a dissipative contribution
multiplied by the couples of parameters $k_{n}$, $k_{t}$ and $\gamma_{n}$,
$\gamma_{t}$, respectively. In the normal force ${\boldsymbol{F}}^{N}_{ij}$ we
can see an elastic term that descends from the Hertzian theory of contact
mechanics [27] characterized by a non-linear dependence on the displacement.
The elastic part of the tangential force consists on a memory term that takes
into account the past history of the tangential displacement whose
infinitesimal contribute is given by $\textbf{ds}_{ij}(t)$. Regarding the
cases in Eq. 1b, they implement the Coulomb friction condition with a
coefficient $\mu$. When it is satisfied
($|{\boldsymbol{F}}^{\text{hist}}_{ij}|\leq|\mu{\boldsymbol{F}}^{N}_{ij}|$),
the typical solid-on-solid sliding force is realized on the contact surface.
For more details and possible extensions of the Hertz-Mindlin contact forces
implemented in LAMMPS see [28].

For the purpose of the main text it is important to note that the tangential
force is proportional to the normal overlap $\xi_{ij}$ below the Coulomb
threshold and to the normal force modulus $|F_{ij}^{N}|$ above it. This
motivate the analysis of Fig. 4e of the main text.

##  II Details and calibration of the simulation

Here we discuss how the parameters of Eqs. 1 and the simulation time step has
been set in relation to the properties of the materials in play. We will
provide also some additional details on the numerical setup to complete the
description done in the main text. For all the simulated geometries we
consider steel grains in a container made of plexiglass (PMMA). The elastic
coefficients for two different materials in contact can be directly derived as

$k_{n}^{ij}=\frac{4}{3}\left(\frac{1-\nu_{i}^{2}}{Y_{i}}+\frac{1-\nu_{j}^{2}}{Y_{j}}\right)^{-1}\quad
k_{t}^{ij}=8\left(\frac{2-\nu_{i}}{G_{i}}+\frac{2-\nu_{j}}{G_{j}}\right)^{-1}$
(2)

where $Y$ is the Young modulus, $G=\frac{Y}{2(1+\nu)}$ the shear modulus and
$\nu$ is the Poisson ratio [24, 27, 29]. When dealing with just one species we
define $k_{n/t}=k^{ii}_{n/t}$. Direct formulas for $\gamma_{n}$ and
$\gamma_{t}$ are lacking but it is a common strategy to choose them verifying
_a posteriori_ the good agreement with experimental data [30] . A widely
accepted criterion to fix the simulation time step is to choose it as a
fraction of the Rayleigh time namely the time that a superficial acoustic wave
takes to cross a single grain. It is related to the characteristics of the
material in this way:

$t_{\text{ray}}=\dfrac{\pi
R_{\text{min}}\sqrt{2\left(1+\nu\right)\rho/Y}}{0.163\nu+0.8766},$ (3)

where $\rho$ is the grain density and $R_{\text{min}}$ the radius of the
littlest grain in the system [31]. With this choice the contact dynamics is
properly resolved because $t_{\text{ray}}$ is usually much smaller than the
typical collision times. In our simulation the simulation time step is
$dt=0.2t_{\text{ray}}$

In Tab. 1 we report the numerical values of the parameters where apex $s$ is
referred to steel and $p$ to plexiglass. In order to reduce obstructions of
the collective drift in the $x$ direction, we removed the tangential
interaction between the vertical walls and the grains of the 2D setup by
setting the related $k_{t}^{sp}$, $\gamma_{t}^{sp}$ and $\mu^{sp}$ to zero
(the values reported in the table refer to the interaction with the bottom
plate). It is worth noting that we reduce the Young modulus of steel and PMMA
of three and two order of magnitude respectively (their real values are
$Y^{s}\simeq 200$ Gpa and $Y^{p}\simeq 3$ Gpa). The reason is the need to have
a bigger $dt$ in order to decrease the simulation time (note that
$t_{\text{ray}}\propto Y^{-1/2}$). This is a general strategy in DEM
simulation where it is possible to have a comparison with experimental data
and verify that the softening of the grains doesn’t affect the phenomena under
interest. We successfully applied this procedure in the numerical study of the
3D cylindrical system [plati2019, 17, 32] where an experimental reference was
present so we find reasonable to use this method also in the more ideal setups
analyzed in the main text (i.e. the 2D and the cubic 3D ones). Also the
numerical values used for $\gamma_{n}$ and $\gamma_{t}$, that cannot be fixed
from the material properties, are justified _a posteriori_ with this same
criterion.

$Y^{s}$ | 210 Mpa
---|---
$\nu^{s}$ | 0.293
$\mu^{ss}$ | 0.5
$k_{n}^{ss}$ | 153$\times 10^{6}$ Pa
$k_{t}^{ss}$ | 190$\times 10^{6}$ Pa
$\gamma_{n}^{ss}$ | 3.1$\times 10^{4}$ $(\text{sm})^{-1}$
$\gamma_{t}^{ss}$ | 9.3$\times 10^{3}$ $(\text{sm})^{-1}$
$Y^{p}$ | 50 Mpa

$\nu^{p}$ | 0.37
---|---
$\mu^{sp}$ | 0.5
$k_{n}^{sp}$ | 617$\times 10^{5}$ Pa
$k_{t}^{sp}$ | 725$\times 10^{5}$ Pa
$\gamma_{n}^{sp}$ | 1.2$\times 10^{7}$ $(\text{sm})^{-1}$
$\gamma_{t}^{sp}$ | 1.0$\times 10^{5}$ $(\text{sm})^{-1}$
$t_{\text{ray}}$ | 6.75$\times 10^{-5}$ s
$dt$ | 1.35$\times 10^{-5}$ s

Table 1: Numerical values for the material properties and the coefficients of
the visco-elastic interaction. We also report the time step $dt$ used for the
simulations. This numerical values specifically refers to the 2D and the 3D
cubic setups. The vertical walls of the 2D setup have
$k_{t}^{sp},\gamma_{t}^{sp},\mu^{sp}=0$ so the don’t exert any tangential
force on the grains. This is the reason why in the Newton equation of the CM
reported in the main text only the bottom particles are considered. For the
cylindrical system we refer to the Supplemental Materials of [plati2019].

Figure 5: a) Trajectories in the $xy$ plane of six grains during a time
interval $\Delta t$ of a simulation in the cylindrical setup shaken at
$\Gamma=39.8$. The overall dynamics can be described as the superimposition of
a slow rigid-body rotation on individual fast vibratons. The black dot
represents the dimension of a grain in scale with the graph: we see that the
amplitude of the fast fluctuations is a small with respect to the grain
diameter b) Trajectories in the $xy$ plane of a substet of 37 grains (lines)
and the center of mass (symbols) in the cubic geometry shaken at
$\Gamma=10.5$.

##  III Drfit trajectories in 3D realistic setups

Collective drifts in realistic 3D setups are the main focus of two former
publications. There, a connection between this phenomenon in the cylindrical
setup and the occurrence of anomalous diffusion for a tracer has been
established [plati2019]. We also proposed phenomenological models to reproduce
the relevant properties of the drifts [17]. In the main text of this paper we
then concentrate on 2D setups, while 3D goemetries are used as test-ground for
the analysis of correlations between structure and dynamics. Nevertheless, in
order to show explicitly how the phenomenon under study is manifested in these
more realistic conditions, we report here two plots of drift trajectories in
the cylindrical setup and in the cubic one (Fig. 5).

##  IV Persistent Motion of Single Particles’ Angular Velocity

In this section we show that the phenomenon discovered in [18] (i.e.
persistent rotation of vertically shaken disks around their own axis) is
observed also in all the setups considered in our study.

Figure 6: Trajectories of the single grains absolute angles around the
$y$-axis for the following cases: polydispersed packing with hard (a) and
periodic (b) boundary conditions on $x$, monodispersed ordered packing with no
defects (c), ordered packings with the TR (d) and RL (e) specular
configurations of defects, ordered packing with the symmetric SD configuration
of defects (f). All the data refer to simulations with $f=100$ and
$\Gamma=0.7$ except the first panel with $\Gamma=0.85$. To see the
correspondent behavior of the collective motion see Fig. 2 of the main text.

We concentrate just on the main observable related to this behavior namely the
trajectories of the absolute angle $\theta_{\alpha
i}(t)=\int_{0}^{t}\omega_{\alpha i}(t^{\prime})dt^{\prime}$ (Fig. 6). Unlike
disks, spherical grains can rotate around the three axis $\alpha=\\{x,y,z\\}$.
For a direct comparison between the single particle dynamics and the
collective one, we suggest to compare Fig. 6 with Fig. 2 of the main text. We
start by the most similar condition with respect to [18] considering a 2D
polydispersed vertical layer confined with hard walls both in the $x$ and the
$y$ direction. In Fig. 6a we plot $\theta_{yi}(t)$ for all the 60 grains in
the system. We can see that the single particle rotation around the $y$ axis
(i.e. the equivalent that we would have had with disks) shows a persistent
motion for almost all the grains in the system. Here, the presence of hard
walls in the $x$ direction prevents the occurrence of the collective drift.
Nevertheless, in Fig. 6b, we report the same analysis performed in the
equivalent setup with PBC on $x$ that allow the formation of the global
translation: here is clear that the single particle persistent motion and the
collective one can be observed at the same time. Panel c of the same figure
shows the results obtained for the ordered monodispersed system with no
defects. We see that there is not any evident long memory effect in the single
particle rotational dynamics consistently with the lack of a collective
translational drift. Figs. 6d-e refer to two ordered packings with specular
defects and panel f to the system with a configuration of defects symmetric
w.r.t. $z$. Also in all these cases the persistent motion of the grains
rotation is observed but we note two additional remarkable facts: i) the
distribution of the trajectory verses for TR and TL is approximately specular
ii) the SD case shows long memory effects just in the single particle dynamics
(we recall that no net collective drift is observed for SD). This last
behavior is consistent with the claim of [18] that the rotational transport is
triggered by local symmetry breakings of the granular packing. The SD
configuration is indeed globally symmetric but creates an asymmetric cage of
nearest neighbours for many bulk particles. Regarding the mirroring of the
trajectories going from TR $\to$ TL, it suggests an intriguing relation
between the distribution of the $\omega_{yi}$s, the global asymmetry of the
packing and consequently the persistent translational motion. In the section
below we sketch some ideas on that.

The trajectories of $\theta_{\alpha i}(t)$ for a subset of grains in the 3D
setups are reported in Fig. 7. We note that persistent single grains rotations
are present in all the cases shown. In the cylindrical setup (panel a and b)
we concentrate on $\theta_{zi}(t)$ instead of $\theta_{yi}(t)$ because, in
this geometry, despite the vertical vibrations, the equivalent hard wall that
is parallel to the free direction is the lateral one. Regarding the 3D cubic
setup where bot the $x$ and the $y$ directions are free we observe persistent
dynamics for both $\theta_{yi}(t)$ and $\theta_{xi}(t)$.

Figure 7: Trajectories of the single grains absolute angles for the
cylindrical setup with $\Gamma=19.5$ (a) and the cubic one with $\Gamma=10.5$
(b-c). In both cases $f=200$.

Figure 8: Trajectories of the single grains absolute angle in the bottom layer
for TR/TL configuration with PBC at $\Gamma=0.7$ (a), $\Gamma=1.25$ (c) and
$\Gamma=0.7$ without PBC (d). In this last panel, the inset contains all the
trajectories in order to show that the absence of persistent rotations regards
mainly the bottom layer. In panel b, in order to check the pure rolling
condition, we show the ratio of the averaged velocity at the contact point and
the averaged translational one of the grains center for both TL and TR. All
the data refer to simulations with $f=100$

###  IV.1 Connection with global translation

The mirroring of the $\omega_{yi}$ distribution going from TR to TL hinted by
Fig. 6d-e suggests to study the single particle angular velocity of the bottom
particles. Indeed, the total tangential force exerted by the bottom plate to
the system depends on these degrees of freedom too; so they can be very
important in order to understand the origin of the global dynamical asymmetry
discussed in the paper (see Fig. 3e-h of the main text). In Fig. 8a we show
the $\theta_{yi}$s of the bottom particles at $\Gamma=0.7$ for the two
mirrored configurations. We see that, referring to the orientation of the
y-axis of Fig. 1 in the main text, all the bottom particles of TR and TL
rotate respectively in the negative (anti-clockwise) and positive (clockwise)
verse. Remarkably, these verses correspond to the ones that, in case of pure
rolling, would involve a translation with the same sign of the collective
motion observed for the two configurations (i.e. negative for TR and positive
for TL). From this analysis one could be tempted to conclude that the
collective translation is simply the consequence of the coordinated rolling of
the bottom particles acting as the ”wheels” of the packing. Nevertheless, the
real situation is quite more complex. First of all we can see from Fig. 8b,
where the ratio between $R\omega_{yi}$ and $v_{xi}$ of the bottom particles is
reported, that the condition of pure rolling ($\langle
R\omega_{yi}\rangle/\langle v_{xi}\rangle\simeq 1$) is only partially
satisfied signalling that also sliding is occurring. Moreover, going toward
higher $\Gamma$ we don’t observe anymore such a correspondence between the
rolling of the bottom particles and the collective translation. Indeed, in
Fig. 8c we can see that the bottom $\theta_{yi}$s for $\Gamma=1.25$ rotate in
both the verses within the same configuration. We verified that in this case
the ratio $\langle R\omega_{yi}\rangle/\langle v_{xi}\rangle$ is even more
different from the unity (not shown). It is worth noting that also in [18] two
different regimes of rotations at low and high $\Gamma$ were found. As a final
analysis we verified that, in simulations of TR and TL _without_ periodic
boundary conditions along the x-axis, the coordinated rolling of the bottom
particles disappears and their motion is not persistent (Fig. 8d). However,
the long memory effect of the rotations is very strong on the upper layers
(see the inset).

After all these remarks, we can say that the slow collective translations and
the persistent single particle rotations influence each other. But neither of
these phenomena can be explained as the simple consequence of the other. In
all the cases in which the collective motion is observed (i.e. when a free
direction of motion is present) the persistent rotations are observed too.
When the global translation is hindered by hard boundaries we still have
steady rotations but with different properties with respect to the ones
observed in presence PBC. As an example, we have seen that the rotational
dynamics of the bottom particles for TR and TL qualitatively changes according
to the presence or absence of PBC (compare panel a and d of Fig. 8).

To conclude, we point out that our numerical study of long memory effects in
the single particles rotational degrees of freedom for spherical grains in
both 2D and 3D geometries represents a considerable extension of the
experimental of work reported in [18]. We have provided many evidences about
the robustness of this phenomenon and suggested some interesting connections
with the collective motion on which the main text is focused on. Deepening the
mechanisms underlying this interplay of memory effects, its origin rooted in
the structure and the resulting asymmetric interaction with the bottom plate
represents a promising perspective that we reserve for future works.

##  V A stick-slip ratchet model

Here we present the stick-slip ratchet model mentioned in the main text. Its
aim is to further support the interpretation of the phenomenon under study as
a ratchet effect. Our model describes the horizontal dynamics of the granular
medium by reducing it to two coarse-grained variables $x$ and $x_{u}$,
identifying the average position of the lower and the upper region
respectively. The dynamics of these coupled variables is defined in such a way
to mimic the behaviour of two adjacent layers of the packing (see the main
text). Then, we find reasonable to introduce in this model four crucial
ingredients i.e. frictional forces, sinusoidal driving, spatial asymmetry and
stick-slips. Considering for simplicity an overdamped dynamics, we write down:

$\displaystyle\dot{x}=-\frac{\mathcal{K}(\Delta x)}{\gamma}\Delta
x+A_{x}\cos(2\pi ft)+\sqrt{2D}\eta(t)$ (4a)
$\displaystyle\dot{x}_{u}=s\frac{\mathcal{K}(\Delta x)}{\gamma_{u}}\Delta x$
(4b)

where $\Delta x=x-x_{u}$, $\eta(t)$ is a Gaussian white noise with zero mean
and unitary variance: $\langle\eta(t)\eta(t^{\prime})\rangle=\delta(t-t’)$,
and $s=0,1$ denotes the state of the system (stick $s=0$, slip $s=1$). The
reminiscence of viscous friction in the underdamped dynamics is given by the
mobility of the two variables $\gamma^{-1}$ and $\gamma_{u}^{-1}$. Regarding
$\mathcal{K}(\Delta x)=\left[1+\text{sign}(\Delta x)\epsilon\right]k$, it
represents a non-isotropic stiffness whose degree of asymmetry is tuned by the
adimensional parameter $-1<\epsilon<1$. From Eq. (4)a we see that the external
driving (that directly acts just on the bottom layer) is modelled as a noise
superimposed to a sinusoidal. Indeed, we find reasonable to assume that the
effective force exerted in the horizontal direction keeps the same periodicity
$1/f$ of the vertically vibrated plate but with a different amplitude $A_{x}$.
The presence of noise accounts from the randomness introduced by the vertical
$\rightarrow$ horizontal conversion of the excitation (it is mediated by the
non linearities of the plate-grain and grain-grain contact forces). We also
point out that in dense vibrofluidized granular systems, the power spectral
density of orthogonal variables with respect the external forcing often
exhibits a peak at the driving characteristic frequencies (Fig. 10c) (see also
experimental spectra in [15]). The non-isotropic stiffness $\mathcal{K}(\Delta
x)$ is the simplest way to introduce a spatial symmetry breaking in our model
without directly defining a biased external force. It is further justified by
a preliminary numerical analysis that suggests an effective asymmetric force
between layers also in the DEM simulations (Fig. 9a).

In order to implement stick-slips we use the following prescription. We define
the high-strain region (HS) to be $\Delta x>\frac{\gamma
F^{*}}{k(1+\epsilon)}\cup\Delta x<\frac{-\gamma F^{*}}{k(1-\epsilon)}$ and the
low-strain (LS) region to be $|\Delta x|<\frac{\gamma F^{*}}{Mk}$ with
$M>\max(1+\epsilon,1-\epsilon)$, while the rest of possible values of $\Delta
x$ are in the medium-strain region (MS). The state variable $s$ goes into
$s=1$ when $\Delta x$ goes into HS, while it goes into $s=0$ when $\Delta x$
enters LS. In MS, $s$ does not change, it depends upon its past (all
simulations start in MS with $s=0$). If $s=0$ the upper variable does not
move, while by tuning $\gamma_{u}^{-1}\gg\gamma^{-1}$ we have that as $s$ is
switched to 1 a very fast slip of $x_{u}$ towards $x$ occurs. Looking at the
behavior of two adjacent layers in a granular packing (Fig. 9b), we realize
that the stick-slips of our model are a crude but reasonable approximation of
the one observed in the DEM simulations.

Eqs. (4) have been numerically integrated [33] and the stick-slip behavior of
$|\Delta x|$ is compared with the one of two adjacent DEM-simulated layers in
Fig. 9c, while Fig. 9d shows a trajectory of $x,x_{u}$ in the deterministic
case ($D=0$). We can clearly observe a drift+oscillation dynamics of both
variables.

Figure 9: a) Evidence of an effective asymmetric force between two adjacent
layers of a granular packing in DEM simulations with $\Gamma=1.25$. Points are
obtained from a scatter plot with bins, lines are linear interpolation of the
data cloud. We note that the asymmetry is correctly inverted if the defects
are reflected with respect to the $z$-axis. b) Zoom of the average position of
the bottom layer and the one just above it as a function of time during a
stick-slip. We note that the velocity of the upper layer undergoes a sudden
increase once a sufficient displacement between the two layers is reached.
Black lines gives an idea of how the stick-slips implemented by our model
would behave in this condition. c) Stick-slip dynamics for the difference of
the mean $x$-coordinates of two adjacent horizontal layers in a granular
simulation and in the relative motion $\Delta x$ from the model. For these
last two panels we considered a BL configuration with $\Gamma=1.25$. d)
Typical trajectories of the model with drift and oscillations. e) Signal to
noise ratio of $\dot{x}_{u}$ as a function of $A_{x}$. For $D=10^{-3}$ the
typical shape of the curve is close to the deterministic case. Values of the
fixed parameters: $k=1$ s-2, $\gamma=1$ s, $f=2$ s-1, $\gamma_{u}^{-1}=500$ s,
$M=10$, $F^{*}=1$.

In Fig. 9e we show the signal to noise ratio of $\dot{x}_{u}$ as a function of
$A_{x}$. We note that after an amplitude threshold (related to $F^{*}$), the
system exhibits a net motion: by setting $F^{*}\gg A$ the system is always
trapped with a zero mean velocity. Also in the opposite limit $F^{*}=0$ (no
stick-slip) there is no drift. The behavior of $\dot{x}_{u}$ VS $A_{x}$
reveals a rich phenomenology, resembling granular simulation (Fig. 2c-d) of
the main text: we have drift’s sign changes as $A_{x}$ changes with the same
$\epsilon$ and reflection of the curves as $\epsilon\to-\epsilon$. For
$\epsilon=0$ there is no mean drift, as for symmetric packings. We stress that
the essential phenomenology of the model in the regime where
$\gamma^{-1}\ll\gamma_{u}^{-1}$, $f^{-1}\ll\gamma_{u}^{-1}$,
$A_{x}/f\gg\sqrt{D/f}$ is contained in the deterministic dynamics, in contrast
with the ancestor model [45]. A detailed study of the model, however, goes
beyond the scope of this work.

To summarize, our model demonstrates that frictional forces, sinusoidal
driving, spatial asymmetry and stick-slips are minimal ingredients to realize
the specific kind of ratchet effect that reflects the behaviour of dense
vibrated granular packings. It consists of a rectification of unbiased
fluctuations whose verse is determined by the interplay of internal
asymmetries and external forcing.

##  VI Size effects

Here we report some supplemental informations about the study of the system’s
size effects. In order to better clarify the different ways in which we
enlarged the size of the granular packings, we provide in Fig. 10a some
snapshots the simulations. We show both ordered packings with defects at
different concentration $c$ and random packings.

As mentioned in the main text, the characterization of the drift properties
for large random packings requires a bit more efforts with respect what is
done for ordered packings with defects. This is due to the fact that we need
to estimate the asymptotic behaviour of $D$ in addition to $\langle
V_{x}^{CM}\rangle$. In order to explain the procedure used to estimate
$D$/$\langle V_{x}^{CM}\rangle^{2}$, we need to point out that the global
dynamics of the packings can be described as a sum of two Ornstein-Uhlenbeck
processes with two well separated characteristic times [17]. The slow
component of the motion represents the global drift with very long
correlations times (i.e. it decorrelates after many rearrangements that are
rare events) while the fast component represents the fluctuations around the
slow drift originated from the local vibrations of the grains. As we expect
from this modelization and shown in Fig. 10b, the power spectral density of
the CM velocity takes the form of two superimposed Lorentzian curves:

$S(\bar{f})=\lim_{t\to\infty}\frac{1}{2\pi
t}\left|\int^{t}_{0}V_{x}^{CM}(t^{\prime})e^{2i(\pi\bar{f})t^{\prime}}dt^{\prime}\right|^{2}=L_{1}(\bar{f})+L_{2}(\bar{f})$
(5)

where
$L_{i}(\bar{f})=\frac{D_{i}}{2\pi}\left[1+(2\pi\bar{f}\tau_{i})^{2}\right]^{-1}$.
The low frequency decay observed in Fig. 10b for $\bar{f}\lesssim 0.1$ Hz
corresponds to the $1/\bar{f}^{2}$ decay of the slow component of the motion
while for $\bar{f}\gtrsim 0.1$ Hz the spectrum is dominated by the Lorentzian
representing the fast component. The diffusion coefficient of the fast
fluctuations around the drift is then proportional to the height of the mid
frequency plateau. Indeed, for a single Lorentzian the diffusion coefficient
is $D_{i}=2\pi L_{i}(0)$. Operatively, we considered only the portion of the
spectrum at $\bar{f}\gtrsim 0.1$ Hz and we fitted it with a single Lorentzian
$L_{1}(\bar{f})$ (see Fig. 10c). Then, $D$ is estimated as $D=D_{1}=2\pi
L_{1}(0)$.

Figure 10: a) Snapshots of the packings used for the study on the system’s
size effect. We consider cases where $\mathcal{N}=6$. In the vertical packing
the highest module is without defects because we find that density
fluctuations in the vertical direction tends to destroy the configurations of
defects in that region of the system. b) Typical shape of the velocity power
spectral density on the CM. We plot the case of a random packing with
$\mathcal{N}=24$ and $\Gamma=1.5$. c) Fitted procedure for the estimate of $D$
from the plateau of the fast component of the spectra fitted with a single
Lorentzian. For each $\mathcal{N}$ we considered ten independent realizations
of random packings. The values of $D$ reported in the main text are calculated
from the average values of the plateau whose heights are highlighted with
dashed black curves in the figure. The peaks at $\bar{f}\sim 100$ Hz are
related the period of the external driving.

##  VII Direct visualization of the simulations

In addition to this supplemental materials, we also provide two videos of the
simulations. Video1.mp4, shows a subset of a trajectory for a random packing
with PBC at high and low $\Gamma$. In the first case the system is completely
fluidized, in the second one it performs a global translation on the
$x$-direction. Particles rapidly vibrate in their cage and follow the
collective motion on a slow time scale. Video2.mp4 contains videos for ordered
packing with defects. Two configurations for two values of $\Gamma$ are
considered in order to see that a sign change of the drift verse occurs
reflecting the defects configuration on the $z$-axis but also varying $\Gamma$
with the same defects configuration. Both the videos and the snapshots in Fig.
1 of the main text has been realized with OVITO [34].


© 2024 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

# Multi-task Learning-based Joint CSI Prediction and Predictive Transmitter
Selection for Security

Shashi Bhushan Kotwal1,2, Chinmoy Kundu2 Sudhakar Modem3, Holger Claussen4,
and Lester Ho5 {1shashi.kotwal<EMAIL_ADDRESS>{2chinmoy.kundu, 4holger.claussen<EMAIL_ADDRESS>13 Department of
EE, Indian Institute of Technology Jammu, Jammu & Kashmir, India  1School of
E&CE, Shri Mata Vaishno Devi University, Katra, Jammu & Kashmir, India
24Wireless Communications Laboratory, Tyndall National Institute, Dublin,
Ireland  4School of Computer Science and Information Technology, University
College Cork, Cork, Ireland  4Trinity College Dublin, Dublin, Ireland

###### Abstract

In mobile communication scenarios, the acquired channel state information
(CSI) rapidly becomes outdated due to fast-changing channels. Opportunistic
transmitter selection based on current CSI for secrecy improvement may be
outdated during actual transmission, negating the diversity benefit of
transmitter selection. Motivated by this problem, we propose a joint CSI
prediction and predictive selection of the optimal transmitter strategy based
on historical CSI by exploiting the temporal correlation among CSIs. The
proposed solution utilizes the multi-task learning (MTL) framework by
employing a single Long Short-Term Memory (LSTM) network architecture that
simultaneously learns two tasks of predicting the CSI and selecting the
optimal transmitter in parallel instead of learning these tasks sequentially.
The proposed LSTM architecture outperforms convolutional neural network (CNN)
based architecture due to its superior ability to capture temporal features in
the data. Compared to the sequential task learning models, the MTL
architecture provides superior predicted secrecy performance for a large
variation in the number of transmitters and the speed of mobile nodes. It also
offers significant computational and memory efficiency leading to a
substantial saving in computational time by around 40 percent.

###### Index Terms:

LSTM, multi-task learning, physical layer security, time-series prediction,
transmitter selection.

##  I Introduction

Transmitter selection (TS) is an effective way to achieve diversity gain in a
multi-transmitter system. Based on available channel state information (CSI),
various transmitter selection schemes are implemented to utilize this
diversity gain to improve the secrecy of communication networks [1, 2, 3, 4,
5]. However, these TS techniques are only applicable in slow-fading channels.
In a practical mobile communication scenario, such as vehicle-to-everything,
high-speed train, unmanned aerial vehicles, etc., the acquired CSI rapidly
becomes obsolete due to the Doppler effect [6].

In a rapidly varying channel, the actual CSI during transmission may differ
from its historical value applied to select a transmitter. Transmission by a
wrongly selected transmitter due to outdated CSI can degrade the secrecy
performance of a system employing TS [6, 7]. To alleviate this problem,
transmitter selection is proposed based on predicted CSI in [6]. The CSI is
predicted using a traditional linear autoregressive (AR) stochastic model. The
results show that transmitter selection based on predicted channels
significantly improves the bit error compared to selection using outdated CSI.
The accuracy of the traditional AR approaches is bounded by the amount of
historical CSI data and the computational complexity is proportional to the
square of AR model order [8].

In contrast to the traditional AR models, machine learning (ML) models can
provide more accurate CSI prediction by exploiting the inherent
characteristics and hidden relationships within historical CSI data via non-
linear approaches [8, 7, 9, 10, 11]. Using a Long Short-Term Memory
(LSTM)-based deep learning ML model, highly accurate CSI is predicted which is
subsequently used for relay selection in case of outdated CSI in [7].

CSI prediction and predictive transmitter selection (transmitter or relay) are
different tasks. The task of transmitter selection in [6] or the predictive
relay selection in [7] rely on a two-step sequential approach (CSI prediction
task followed by transmitter selection task). In contrast, in [12] a single-
step approach is proposed using a convolutional neural network (CNN) to
predict selected transmitter directly from historical CSI values. Directly
predicting selected transmitter from the historical CSI can avoid two-step
computation, however, it does not solve the problem of CSI prediction. CSI
prediction is also required for supporting various link-adaptive schemes, such
as adaptive modulation and coding, power adaptation, etc. and performance
prediction. Additionally, CSI prediction can reduce the need for channel
estimation, pilot transmissions, and feedback overhead, thus improving the
channel usage efficiency.

The articles [6, 7, 12] have focused on a single task at a time. However,
simultaneously performing both CSI prediction and predictive transmitter
selection provide a suitable solution for meeting mobile networks’
computational and memory constraints. It can reduce the training data
requirement compared to that required when CSI prediction and predictive
transmitter selection tasks are implemented in sequence. The multi-task
learning (MTL) framework can offer the solution for simultaneously learning
multiple tasks in ML [13, 14, 15]. MTL framework leverages the inherent
relation among training data across different tasks, allowing the model to
exploit these relations and potentially improve the performance of all tasks.
By sharing the common parameters across different tasks, the training model
complexity can also be reduced in MTL.

Motivated by the above discussion, a system with multiple transmitters and
eavesdroppers in a mobile communication scenario is considered for joint CSI
prediction and predictive transmitter selection for secrecy enhancement. It is
to be noted that articles [6, 7, 12] though considered predictive TS under
outdated CSI condition, these articles did not consider secrecy. Moreover,
none of these articles implemented multiple tasks jointly. Although [9]
considered secrecy, the work was restricted to the CSI prediction task only.
The main contributions of our work are outlined below.

  * •

The joint CSI prediction and predictive transmitter selection problem is
solved using MTL-based ML architecture that exploits the temporal relations in
the CSI using LSTM model.

  * •

We show that while both the LSTM and CNN architectures can be employed for
predicting the CSI and the optimal transmitter jointly, LSTM outperforms CNN
due to its superior ability to capture temporal features in the data.

  * •

We show that the MTL-based model has better predicted secrecy over a large
relative node speed and number of transmitters as compared to the sequential
task learning model.

  * •

In contrast to the sequential learning model, the MTL-based joint CSI
prediction and transmitter selection architecture is computationally and
memory efficient. It can save computation time by around 40 percent.

The rest of the paper is organized as follows. The system model is described
in Section II. The LSTM model for joint CSI prediction and predictive
transmitter selection is described in Section III, results are presented in
Section IV and conclusions are provided in Section V.

Notation: $|\cdot|$ and $||\cdot||_{2}$ denote magnitude and $L_{2}$ norm
operators, respectively. $\mathbb{E}\left[\cdot\right]$ and $R_{hh}(\cdot)$
denote the expectation and autocorrelation operator, respectively.
$\mathbf{T}$ denotes the transpose operator and
$\left[x\right]^{+}=\max\\{0,x\\}$ where $\max\\{0,x\\}$ denotes maximum of
$0$ and $x$.

##  II System Model

Figure 1: A system with $K$ transmitters and $N$ eavesdroppers for joint CSI
prediction and predictive transmitter selcetion.

We consider a mobile wiretap communication system as in Fig. 1 consisting of
$K$ number of transmitters $\mathrm{S}^{(k)}$, where $k\in\\{1,\ldots,K\\}$, a
legitimate destination D, and $N$ passive eavesdroppers $\mathrm{E}^{(n)}$,
where $n\in\\{1,\ldots,N\\}$. It is assumed that all nodes are equipped with a
single antenna. The $k$-th channel between $\mathrm{S}^{(k)}$ and $\mathrm{D}$
for each $k\in\\{1,\ldots,K\\}$ is assumed to undergo an independent and
identically distributed (i.i.d.) fading. Due to the Doppler effect in the
mobile scenario, the destination and eavesdropping channels are time-varying.
In this paper, the time-varying channels are modelled based on Clarke’s model
[16, 17]. Thus, the impulse response of the time-varying channel
$\mathrm{S}^{(k)}$-$\mathrm{D}$ can be stated with $M$ number of propagation
paths as

$\displaystyle{h}_{\mathrm{D}}^{(k)}(t)$
$\displaystyle=\frac{1}{\sqrt{M}}\left(Y_{CD}^{(k)}(t)+jY_{SD}^{(k)}(t)\right),$
(1)

where $Y_{CD}^{(k)}(t)$ and $Y_{SD}^{(k)}(t)$ are the quadrature components
and are expressed as

$\displaystyle Y_{CD}^{(k)}(t)$ $\displaystyle=\sum_{m=1}^{M}\cos\left(2\pi
f_{d}t\cos\beta_{m}^{(k)}+\phi_{m}^{(k)}\right),$ (2) $\displaystyle
Y_{SD}^{(k)}(t)$ $\displaystyle=\sum_{m=1}^{M}\sin\left(2\pi
f_{d}t\cos\beta_{m}^{(k)}+\phi_{m}^{(k)}\right),$ (3)

$f_{d}={f_{c}v}/{c}$ is the maximum Doppler frequency, where $f_{c}$ is the
carrier frequency, $v$ is the relative velocity between the transmitters and
the receivers including legitimate destination and eavesdroppers, and $c$ is
the speed of the light; and $\beta_{m}^{(k)}=({2\pi m+\theta_{m}^{(k)}})/{M}$
and $\phi_{m}^{(k)}$ are the angle of arrival and initial phase of the $m$-th
propagation path, respectively, at the destination. Both $\theta_{m}^{(k)}$
and $\phi_{m}^{(k)}$ are statistically independent and uniformly distributed
over $\left[-\pi,\pi\right]$ for each $m$. The autocorrelation function of the
channel complex envelope obtained from (1) is $R_{hh}(\tau)=J_{0}(2\pi
f_{d}\tau)$, where $J_{0}(\cdot)$ is the zero-order Bessel function of the
first kind [17]. The correlation in time between CSI samples is introduced due
to the nature of the time functions defined in (2) and (3). The coherence time
$T_{c}$ of all the channels is $T_{c}={0.423}/{f_{d}}$ [18].

The instantaneous SNR of the channel $\mathrm{S}^{(k)}$-$\mathrm{D}$ is given
by
$\gamma_{\mathrm{D}}^{(k)}(t)=\bar{\gamma}_{\mathrm{D}}|h_{\mathrm{D}}^{(k)}(t)|^{2}$
where $\mathbb{E}[|h_{\mathrm{D}}^{(k)}(t)|^{2}]=1$,
$\bar{\gamma}_{\mathrm{D}}=\alpha_{\mathrm{D}}P_{T}/{(\sigma_{\mathrm{D}})}^{2}$
is the average received SNR at D, $P_{T}$ is the transmitted power,
$\alpha_{\mathrm{D}}$ is the path loss factor inversely proportional to the
distance between $\mathrm{S}^{(k)}$ and D, and $({\sigma_{\mathrm{D}}})^{2}$
is the average noise power at D. Similar to ${h}_{\mathrm{D}}^{(k)}(t)$, the
time-varying channel between $\mathrm{S}^{(k)}$ and $\mathrm{E}^{(n)}$,
$h_{\mathrm{E}}^{(k,n)}(t)$, for each $k\in\\{1,\ldots,K\\}$ and
$n\in\\{1,\ldots,N\\}$ is assumed to undergo an i.i.d. fading. The impulse
response of the channel $\mathrm{S}^{(k)}$-$\mathrm{E}^{(n)}$ with $L$ number
of propagation paths can also be modelled similar to
${h}_{\mathrm{D}}^{(k)}(t)$ as shown in (1). With
$\mathbb{E}[|{h}_{\mathrm{E}}^{(k,n)}(t)|^{2}]=1$, the average SNR at
$\mathrm{E}^{(n)}$ is expressed as
$\bar{\gamma}_{\mathrm{E}}=\alpha_{\mathrm{E}}P_{T}/({\sigma_{\mathrm{E}}})^{2}$
where $\alpha_{\mathrm{E}}$ is the path loss factor inversely proportional to
the distance between $\mathrm{S}^{(k)}$ and $\mathrm{E}^{(n)}$, and
$({\sigma_{\mathrm{E}}})^{2}$ is the average noise power at
$\mathrm{E}^{(n)}$.

To consider the worst possible eavesdropping scenario, we assume that the
eavesdroppers are colluding and perform maximal ratio combining (MRC) of the
received signals to extract confidential information. The combined
instantaneous SNR $\gamma_{\mathrm{E}}^{(k)}(t)$ at the eavesdroppers
corresponding to the $k$-th transmitter is therefore given by

$\displaystyle\gamma_{\mathrm{E}}^{(k)}(t)$
$\displaystyle=\bar{\gamma}_{E}\sum_{n=1}^{N}|h_{\mathrm{E}}^{(k,n)}(t)|^{2},$
(4)

The optimal transmitter that provides the maximum instantaneous secrecy rate
among all the transmitters is selected. The instantaneous secrecy rate of the
system with optimal transmitter selection is defined in bits per channel use
(bpcu) as [4]

$\displaystyle C_{S}^{(k^{*})}(t)$
$\displaystyle=\max_{k\in{1,\ldots,K}}\left[\log_{2}\left(\frac{1+\gamma_{\mathrm{D}}^{(k)}(t)}{1+\gamma_{\mathrm{E}}^{(k)}(t)}\right)\right]^{+},$
(5)

where $k^{*}$ is the optimal selected transmitter. The ergodic secrecy rate
(ESR) of the system for the optimal transmitter selection is then given as [5]

$\displaystyle C_{\mathrm{erg}}$
$\displaystyle=\mathbb{E}\left\\{C_{S}^{(k^{*})}(t)\right\\}.$ (6)

##  III Joint CSI Prediction and Predictive Transmitter Selection

In this section, we formulate the problem of joint CSI prediction and
predictive transmitter selection. We then propose an LSTM-based MTL algorithm
for the solution.

###  III-A Problem Formulation

The joint CSI prediction and predictive transmitter selection is formulated as
an MTL problem where CSI for future time instances and optimal transmitter
that maximizes the instantaneous secrecy rate for those time instances are
jointly predicted through a single ML network based on historical CSI. To
formally define the problem statement, we first define the CSI vector that is
used to propose the ML algorithm as

$\displaystyle\mathbf{h}(t)$
$\displaystyle=\Big{[}\mathbf{h}_{\mathrm{D}}(t),\mathbf{h}_{\mathrm{E}}(t)\Big{]}^{\mathbf{T}},$
(7)

where

$\displaystyle\mathbf{h}_{\mathrm{D}}(t)$
$\displaystyle=\Big{[}h_{\mathrm{D}}^{(1)}(t),\ldots,h_{\mathrm{D}}^{(K)}(t)\Big{]},$
(8) $\displaystyle\mathbf{h}_{\mathrm{E}}(t)$
$\displaystyle=\Big{[}h_{\mathrm{E}}^{(1,1)}(t),\ldots,h_{\mathrm{E}}^{(1,N)}(t),h_{\mathrm{E}}^{(2,N)}(t),\ldots,h_{\mathrm{E}}^{(K,N)}(t)\Big{]},$
(9)

and
$|\mathbf{h}(t)|=\Big{[}|\mathbf{h}_{\mathrm{D}}(t)|,|\mathbf{h}_{\mathrm{E}}(t)|\Big{]}^{\mathbf{T}}$
is the vector of length $K+(K\times N)$ with magnitudes of current CSI samples
at time $t$. The proposed problem can then be formally stated as

$\displaystyle\\{|\mathbf{\hat{h}}(t+1)|,\hat{k}_{1}^{*}\\},\ldots,\\{|\mathbf{\hat{h}}(t+J)|,\hat{k}_{J}^{*}\\}$
$\displaystyle=f\left(|\mathbf{h}(t-T+1)|,\ldots,|\mathbf{h}(t)|\right),$ (10)

where $T$ is the total length of the past CSI samples, $J$ is the length of
the predicted CSI samples, $f(\cdot)$ is an arbitrary function that maps the
$T$ historical CSI samples to the future estimates of the channels
$|\mathbf{h}(t+1)|,\ldots,|\mathbf{h}(t+J)|$ as
$|\mathbf{\hat{h}}(t+1)|,\ldots,|\mathbf{\hat{h}}(t+J)|$ and also jointly
predicts $J$ future selected transmitter indexes $k_{1}^{*},\ldots,k_{J}^{*}$
as $\hat{k}_{1}^{*},\ldots,\hat{k}_{J}^{*}$. The problem is then to design an
ML model that implements the function $f(\cdot)$ in (III-A). The solution is
provided using the LSTM-based MTL algorithm.

###  III-B Model Training Methodology and Multi-task loss function

The prediction of the CSI based on historical CSI is implemented as a time
series prediction problem [19]. We utilize the temporal dependency capturing
feature of LSTM to predict the CSI of the destination and eavesdroppers’
channels. Additionally, we use an LSTM-based time series classifier for
predicting the selected transmitter. For both time series prediction and time
series classification, our proposed LSTM-based model exploits the temporal
correlation in the wireless channel.

The proposed model, named as LSTM-J model, as shown in Fig. 2, utilizes the
MTL framework where a single LSTM layer is used followed by two Fully
Connected Layers each serving different tasks [13]. The Fully Connected
Layer-I assists in CSI prediction whereas, the Fully Connected Layer-II
assists in transmitter selection by providing its output to a softmax layer
that converts the output into the probability distribution over the classes of
the transmitters for transmitter classification.

We provide past CSI data as the input to the network as shown in (III-A) to
train the model depicted in Fig. 2. The input data for the training is
arranged in $N_{data}$ training sets where each set consists of three matrices
$\mathbf{H_{train}}$, $\mathbf{H_{target}}$, and $\mathbf{K_{target}}$. The
matrix $\mathbf{H_{train}}^{K(1+N)\times
T}=\left[|\mathbf{h}(t-T+1)|,\ldots,|\mathbf{h}(t)|\right]$ corresponds to the
historical input data of destination and eavesdropping channel, the matrix
$\mathbf{H_{target}}^{K(1+N)\times
J}=\left[|\mathbf{h}(t+1)|,\ldots,|\mathbf{h}(t+J)|\right]$ corresponds to the
target future CSI values that follow the historical values
$\mathbf{H_{train}}$, and the vector $\mathbf{K_{target}}^{1\times J}$
corresponds to the transmitters selected for $J$ future time instances
according to the target future CSI values $\mathbf{H_{target}}^{K(1+N)\times
J}$. The training is completed in mini batches where each batch contains
$N_{b}$ training sets.

To design the multi-task loss function $L_{mt}$ for training, we choose
individual loss functions, mean square error (MSE) loss function $L_{p}$
corresponding to the CSI prediction task and cross-entropy loss function
$L_{c}$ corresponding to the transmitter classification task and produce a
weighted sum of these two loss functions as MTL loss function. This is
mathematically defined as [14]

$\displaystyle L_{mt}=wL_{p}+(1-w)L_{c},$ (11)

where $0<w<1$ is the weight assigned to $L_{p}$ and

$\displaystyle
L_{p}=\frac{1}{JK(N+1)}\sum_{n_{b}=1}^{N_{b}}\sum_{j=1}^{J}||\mathbf{e}_{n_{b}}^{(j)}||_{2}^{2},$
(12) $\displaystyle
L_{c}=-\sum_{n_{b}=1}^{N_{b}}\sum_{k=1}^{K}p_{n_{b},k}\ln\left(q_{n_{b},k}\right).$
(13)

In (12),
$\mathbf{e}_{n_{b}}^{(j)}=(|\mathbf{\hat{h}}(t+j)|_{n_{b}}-|\mathbf{h}(t+j)|_{n_{b}})$
is the element-wise difference in the predicted CSI
$|\mathbf{\hat{h}}(t+j)|_{n_{b}}$ and target CSI $|\mathbf{h}(t+j)|_{n_{b}}$
for the $j$-th time instance in the $n_{b}$-th training set where
$j\in\\{1,\ldots,J\\}$ and $n_{b}\in\\{1,\ldots,N_{b}\\}$. In (13),
$p_{n_{b},k}$ as well as $q_{n_{b},k}$ represent the predicted probability and
actual probability, respectively, assigned to the $k$-th transmitter in the
$n_{b}$-th training set for the transmitter classification purposes. The
objective of the training is to minimize $L_{mt}$.

###  III-C Comparison with other MTL and sequential ML models

We compare the performance of the proposed LSTM architecture with that of the
CNN-based MTL architecture. In the CNN architecture, the LSTM layer in Fig. 2
is replaced by the combination of a one-dimensional CNN layer, a Rectified
Linear Unit (ReLU) layer, a Layer Normalization layer, and a Global Average
Pooling layer [20]. In the CNN architecture, one-dimensional convolution
filters stride across time to extract temporal features from the available
CSI, enabling the prediction of future CSI and transmitter selection. The
incorporation of non-linearity due to the ReLU layer in CNN architecture
enables the network to learn complex representations of the data. We also
compare the performance of LSTM-J model in Fig. 2 with that of the model in
Fig. 3, denoted as LSTM-S model, where two individual networks, one for CSI
prediction task and the another for transmitter selection task based on
predicted CSI, are employed sequentially.

Figure 2: Layered architecture of the LSTM-J model trained for predicting both
CSI and selected transmitter jointly. Figure 3: Layered architecture of the
LSTM-S model for predicting CSI and selected transmitter sequentially.

###  III-D Performance Metrics

We evaluate the performance of proposed ML models for their CSI prediction
accuracy and the accuracy in transmitter selection. The CSI prediction
accuracy is measured using normalised mean square error (NMSE) between future
estimates of the channels and the actual channels, expressed as [21]

NMSE
$\displaystyle=\frac{\mathbb{E}\\{\sum_{j=1}^{J}\left(|||\mathbf{\hat{h}}(t+j)|-|\mathbf{h}(t+j)|||_{2}^{2}\right)\\}}{\mathbb{E}\\{\sum_{j=1}^{J}||\mathbf{h}(t+j)||_{2}^{2}\\}}.$
(14)

The transmitter selection accuracy is evaluated as the ratio of the number of
instances when $k^{*}=\hat{k}^{*}$ to the total number of instances, expressed
as a percentage.

Using
$\\{|\mathbf{\hat{h}}(t+1)|,\hat{k}_{1}^{*}\\},\ldots,\\{|\mathbf{\hat{h}}(t+J)|,\hat{k}_{J}^{*}\\}$
in (5), we predict the instantaneous secrecy rate
$\hat{C}_{S}^{(k^{*})}(t+1),\ldots,\hat{C}_{S}^{(k^{*})}(t+J)$. The prediction
accuracy of the CSI and selected transmitter determines how close the
estimated ESR $\hat{C}_{\mathrm{erg}}$ is from the actual ESR
$C_{\mathrm{erg}}$.

##  IV Numerical Results

TABLE I: Parameters for the dataset generation and training. Parameter | Value
---|---
Carrier centre frequency $(f_{c})$ | 2GHz
Sampling frequency $(f_{s})$ | 1KHz
Channel fading | Rayleigh
Number of paths | 100
Doppler spectrum |  Clarke’s model [17]
Historical CSI values $(T)$ | 10
Predicted CSI values $(J)$ | 1
Training data observations | 49990
Test data observations | 149990
Mini batch size $(N_{b})$ | 500
Training epochs | 5
Training iterations $(n\\_iter)$ | 500
Learning rate | $\frac{0.005}{(1+0.005\times n_{iter})}$
Optimizer | Adam
Number of hidden units (LSTM) | 200
Number of filters (CNN) | 50
Filter size (CNN) | 6

In this section, we evaluate the performance of the LSTM-based MTL
architecture proposed in Fig. 2. We also compare its performance with that of
a CNN-based MTL architecture. Both these MTL architectures are compared with
the sequential task learning architecture shown in Fig. 3. We use a typical
personal computer with an Intel® Core ™ i7-8700 CPU and 16 GB RAM for training
purposes. We generate a time series CSI dataset using (1) with specifications
given in Table I. In all the figures, LSTM-J and CNN-J denote outputs from the
MTL network architecture as proposed in Fig. 2. On the other hand, LSTM-S and
CNN-S represent the outputs from the sequential network architectures shown in
Fig. 3. The CSI prediction error is shown in terms of NMSE defined in (14) and
transmitter selection accuracy is shown as selection accuracy in the figures.

In Figs. 4(a) and 4(b), we show NMSE and the transmitter selection accuracy,
respectively, for both the LSTM and CNN-based architectures by varying the
weight $w$ in the MTL loss function. The figure tells how to choose $w$ for
the proposed MTL problem. We observe that a higher weight assigned to $L_{p}$
decreases the CSI prediction error, however, if $w$ is more than 0.9, the
transmitter selection accuracy drops. This suggests that the optimal weight
for the problem is around $w=0.9$. However, the CSI prediction error and the
transmitter selection accuracy remain relatively similar in the range of
$0.1\leq w\leq 0.9$. Hence, all the subsequent figures are plotted with
$w=0.9$.

(a) CSI prediction error.

(b) Transmitter selection accuracy.

Figure 4: Performance with $w$ when $\bar{\gamma}_{\mathrm{D}}=30$ dB,
$\bar{\gamma}_{\mathrm{E}}=10$ dB, $K=3$, $N=2$, $v=10$ meters/second.

In Figs. 5(a) and 5(b), we compare the CSI prediction error and the
transmitter selection accuracy, respectively, for both the LSTM and CNN-based
architecture by varying relative speed between transmitter and receiver $v$
from 5 to 50 meters/second. As seen in Fig. 5(a), CSI prediction error
increases with increasing $v$ for all the networks. With increasing $v$, the
Doppler frequency increases, lowering the coherence time of the channel. This
weakens the correlation between CSI samples across multiple coherence
intervals, making it difficult for both LSTM and CNN to learn the channel
parameters. We find that LSTM is better suited for CSI prediction in mobile
scenarios as compared to CNN as it has significantly less CSI prediction
error. This is attributed to the LSTM’s better ability to learn long-term
dependencies in sequential data than CNN [22].

Due to reduced correlation between CSI samples with $v$ as described above, in
Fig. 5(b), we see that the transmitter selection accuracy for both LSTM and
CNN-based architectures also degrades with the increasing $v$. A better
accuracy performance of LSTM is also observed as compared to CNN for a large
mobility range. Finally, we observe that both MTL and sequential architectures
provide similar CSI prediction errors. The transmitter selection accuracy of
the MTL and sequential architectures are also comparable at higher $v$,
however, the same for CNN-J is better than CNN-S at lower $v$.

(a) CSI prediction error.

(b) Transmitter selection accuracy.

Figure 5: Performance with speed when $\bar{\gamma}_{\mathrm{D}}=30$ dB,
$\bar{\gamma}_{\mathrm{E}}=10$ dB, $K=4$, $N=2$.

In Figs. 6(a) and 6(b), the CSI prediction error and the transmitter selection
accuracy, respectively, are plotted for the LSTM and CNN-based architecture by
varying the number of transmitters $K$. In Fig. 6(a), we find that the
increasing $K$ has a negligible effect on CSI prediction accuracy in all
schemes. This observation is reasonable as transmitters are independent of
each other. We also note that both MTL and sequential architectures provide
similar CSI prediction error with $K$.

In Fig. 6(b), the selection accuracy degrades with increasing $K$. This is due
to the growth in the learnable parameters with $K$ making it difficult to
predict the selected transmitter accurately. We observe that while CNN-J is
always better than CNN-S, LSTM-J outperforms LSTM-S after a certain number of
transmitters. The advantage of MTL-based architectures (LSTM-J and CNN-J) in
transmitter selection accuracy is that they can predict selected transmitters
directly based on the historical CSI. The potential errors in sequential
networks (LSTM-S and CNN-S) come from the fact that the individual tasks are
performed sequentially. The accuracy of transmitter classification is thus
dependent on the accuracy of CSI prediction.

(a) CSI prediction error.

(b) Transmitter selection accuracy.

Figure 6: Performance with the number of transmitters when
$\bar{\gamma}_{\mathrm{D}}=30$ dB, $\bar{\gamma}_{\mathrm{E}}=10$ dB, $N=2$,
and $v=10$ meters/second.

In Figs. 7(a) and 7(b), we depict the estimated ESR of proposed MTL-based
networks (LSTM-J and CNN-J) versus $v$ and $K$, respectively. The estimated
ESR is obtained by utilizing the estimated CSI and the predictive optimal
transmitter produced by the MTL-based network in (5). The ESR of the MTL-based
network is compared with that of the hypothetical perfect system with the
perfectly estimated CSI and accurately selected optimal transmitter. It is
clearly visible from both figures that the LSTM-based network is better than
the CNN-based network in following the perfect system performance.

(a) When $K=4$.

(b) When $v=10$ meters/second.

Figure 7: ESR with $v$ and $K$ when $\bar{\gamma}_{\mathrm{D}}=30$ dB,
$\bar{\gamma}_{\mathrm{E}}=10$ dB, $N=2$.

In Fig. 7(a), the ESR curves of both LSTM and CNN networks deviate
significantly from the perfect system in the high-mobility region (high $v$).
The ESR corresponding to MTL and sequential networks are close to each other
for both LSTM and CNN. These observations are consistent with that in Fig.
5(a) and 5(b). When performance is considered with respect to $K$ in Fig.
7(b), it is found that MTL-based network (CNN-J and LSTM-J) outperforms the
sequential network beyond a certain $K$ in respective LSTM and CNN systems.
This observation is consistent with that in Fig. 6(a) and 6(b).

A comparison of training and testing time required for the MTL and sequential
model is shown in Table II. It is found that the training and testing time is
reduced by 40.8 and 39.9 percent, respectively, in LSTM-J as compared to
LSTM-S. A similar observation is also found in CNN-J versus CNN-S. Further,
MTL-based architecture is trained only once, hence, it is also memory
efficient as compared to the sequential task learning model which requires
multiple trainings [23].

TABLE II: Comparison in training and testing time (seconds). Network | Training time | Testing time
---|---|---
LSTM-J | 146.1 | 17.7
LSTM-S | 246.7 | 29.0
CNN-J | 102.5 | 10.1
CNN-S | 137.1 | 4.5

##  V Conclusion

In this paper, we propose a novel approach of joint CSI prediction and
predictive transmitter selection tasks using historical CSI in mobile
scenarios. This approach utilizes LSTM-based MTL architecture that exploits
temporal correlation in CSI and the inherent relationship between tasks.
Compared to the CNN architecture, LSTM provides better performance. Compared
to the sequential task learning model MTL architecture provides superior
predicted secrecy performance for large variations in the number of
transmitters and mobile node speeds. It also delivers better computational and
memory efficiency. The LSTM-based MTL architecture is found to reduce the
computational time by around 40 percent as compared to its sequential task
learning architecture.

## References

  * [1] C. Kundu, T. M. N. Ngatched, and O. A. Dobre, “Relay selection to improve secrecy in cooperative threshold decode-and-forward relaying,” in _Proc. IEEE Global Communications Conference_ , Washington, DC, USA, Dec. 2016, pp. 1–6.
  * [2] C. Kundu and M. F. Flanagan, “Ergodic secrecy rate of optimal source selection in a multi-source system with unreliable backhaul,” _IEEE Wireless Communications Letters_ , vol. 10, no. 5, pp. 1118–1122, Feb. 2021.
  * [3] S. Tripathi, C. Kundu, O. A. Dobre, A.Bansal, and M. F. Flanagan, “Recurrent neural network assisted transmitter selection for secrecy in cognitive radio network,” in _Proc. IEEE Global Communications Conference_ , Taipei, Taiwan, Dec. 2020, pp. 1–6.
  * [4] S. B. Kotwal, C. Kundu, S. Modem, A. Dubey, and M. F. Flanagan, “Transmitter selection for secrecy in a frequency selective fading channel with unreliable backhaul,” in _Proc. IEEE Vehicular Technology Conference_ , Helsinki, Finland, Jun. 2021, pp. 1–7.
  * [5] S. B. Kotwal, C. Kundu, S. Modem, and M. F. Flanagan, “Transmitter selection for secrecy in frequency-selective fading with multiple eavesdroppers and wireless backhaul links,” _IEEE Transactions on Vehicular Technology_ , vol. 73, no. 1, pp. 860–875, Jan. 2024.
  * [6] K. E. Baddour, C. C. Squires, and T. J. Willink, “Mobile channel prediction with application to transmitter antenna selection for alamouti systems,” in _Proc. IEEE Vehicular Technology Conference_ , Montreal, QC, Canada, Feb. 2006, pp. 1–6.
  * [7] W. Jiang and H. D. Schotten, “Predictive relay selection: A cooperative diversity scheme using deep learning,” in _Proc. IEEE Wireless Communications and Networking Conference_ , Nanjing, China, May 2021, pp. 1–6.
  * [8] J. Yuan, H. Q. Ngo, and M. Matthaiou, “Machine learning-based channel prediction in massive MIMO with channel aging,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 5, pp. 2960–2973, May 2020.
  * [9] T. Peng, R. Zhang, X. Cheng, and L. Yang, “LSTM-based channel prediction for secure massive MIMO communications under imperfect CSI,” in _Proc. IEEE International Conference on Communications_ , Dublin, Ireland, Jun. 2020, pp. 1–6.
  * [10] E. Tekgul, J. Chen, J. Tan, F. Vook, S. Özen, and A. Jajoo, “Deep learning-based channel state information prediction with incomplete history,” in _Proc. IEEE Wireless Communications and Networking Conference_ , Austin, TX, USA, May 2022, pp. 447–452.
  * [11] S. R. Mattu, L. N. Theagarajan, and A. Chockalingam, “Deep channel prediction: A DNN framework for receiver design in time-varying fading channels,” _IEEE Transactions on Vehicular Technology_ , vol. 71, no. 6, pp. 6439–6453, Jun. 2022.
  * [12] J. Kim, J. Oh, E.-R. Jeong, and J. Joung, “Convolutional neural network-based transmit antenna selection for UAV-ground station communications with time-varying channels,” _ICT Express_ , vol. 10, no. 1, pp. 90–96, Feb. 2024.
  * [13] R. Caruana, “Multitask learning,” Ph.D. dissertation, Carnegie Mellon University, 1997.
  * [14] R. Cipolla, Y. Gal, and A. Kendall, “Multi-task learning using uncertainty to weigh losses for scene geometry and semantics,” in _Proc. IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , Los Alamitos, CA, USA, Jun. 2018, pp. 7482–7491.
  * [15] A. Rago, G. Piro, G. Boggia, and P. Dini, “Multi-task learning at the mobile edge: An effective way to combine traffic classification and prediction,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 9, pp. 10 362–10 374, Jun. 2020.
  * [16] R. H. Clarke, “A statistical theory of mobile-radio reception,” _The Bell System Technical Journal_ , vol. 47, no. 6, pp. 957–1000, Jul. 1968.
  * [17] C. Xiao, Y. R. Zheng, and N. C. Beaulieu, “Novel sum-of-sinusoids simulation models for rayleigh and rician fading channels,” _IEEE Transactions on Wireless Communications_ , vol. 5, no. 12, pp. 3667–3679, Dec. 2005.
  * [18] B. Sklar, “Rayleigh fading channels in mobile digital communication systems. I. characterization,” _IEEE Communications Magazine_ , vol. 35, no. 9, pp. 136–146, Sep. 1997.
  * [19] G. Liu, Z. Hu, L. Wang, J. Xue, H. Yin, and D. Gesbert, “Spatio-temporal neural network for channel prediction in massive MIMO-OFDM systems,” _IEEE Transactions on Communications_ , vol. 70, no. 12, pp. 8003–8016, Dec. 2022.
  * [20] J. Wang, K. Yu, J. Bu, Y. Lin, and S. Han, “Multi-classification of UWB signal propagation channels based on one-dimensional wavelet packet analysis and CNN,” _IEEE Transactions on Vehicular Technology_ , vol. 71, no. 8, pp. 8534–8547, Aug. 2022.
  * [21] Y. Huangfu, J. Wang, R. Li, C. Xu, X. Wang, H. Zhang, and J. Wang, “Predicting the mumble of wireless channel with sequence-to-sequence models,” in _Proc. IEEE Annual International Symposium on Personal, Indoor and Mobile Radio Communications_ , Istanbul, Turkey, Nov. 2019, pp. 1–7.
  * [22] N. Mohammadi Foumani, L. Miller, C. W. Tan, G. I. Webb, G. Forestier, and M. Salehi, “Deep learning for time series classification and extrinsic regression: A current survey,” vol. 56, no. 9, pp. 1–45, Apr. 2024.
  * [23] C. Zhang, P. Patras, and H. Haddadi, “Deep learning in mobile and wireless networking: A survey,” _IEEE Communications Surveys & Tutorials_, vol. 21, no. 3, pp. 2224–2287, Mar. 2019.