%~Mouliné par MaN_auto v.0.37.1 (6604bb9b) 2025-11-12 13:24:22
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_Sano_20250510.bib}

\usepackage[group-digits = integer]{siunitx}
\usepackage{subcaption}

\newcommand{\solid}{\mathrm{solid}}
\newcommand{\granular}{\mathrm{granular}}
\newcommand{\tot}{\mathrm{tot}}
\newcommand{\Glass}{\mathrm{Glass}}
\newcommand{\Aluminum}{\mathrm{Aluminum}}
\newcommand{\Lead}{\mathrm{Lead}}

\newcommand{\dif}{\mathop{}\!{\operatorfont{d}}}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\bracks}{[}{]}

%%%---------------------------------------------------------------------------

%%% FIGURES

\makeatletter
\renewcommand\p@subfigure{}
\makeatother

%%%---------------------------------------------------------------------------

%%% SYMBOLS

%%% \lesssim
\newcommand{\lesssim}{\mathrel{\raisebox{-2.5pt}{$\stackrel{\raisebox{-1pt}{\scalebox{1}[.8]{<}}}{\sim}$}}}

%%% \gtrsim
\newcommand{\gtrsim}{\mathrel{\raisebox{-2.5pt}{$\stackrel{\raisebox{-1pt}{\scalebox{1}[.8]{>}}}{\sim}$}}}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}

\title{Density wave propagation of momentum and the formation of arches in the undulation of a vertically vibrated thin granular layer}
\alttitle{Propagation des ondes de densité de l'impulsion et formation d'arcs dans l'ondulation d'une couche granulaire mince soumise à des vibrations verticales}

\author{\firstname{Osamu} \lastname{Sano}}
\email{sano@cc.tuat.ac.jp}
\address{Department of Applied Physics, Tokyo University of Agriculture \& Technology, Tokyo 184-8588, Japan}

\begin{abstract}
Experimental and theoretical investigations are made on the aspect ratio of the height $\delta$ to wavelength $\lambda$ of the undulations generated in a dry granular layer confined between two parallel vertical planes under vertical vibration of frequency $f$ and amplitude $a$. We elucidated the propagation of the density wave along the layer, which suggests one or two order smaller macroscopic elastic constant ascribed to a configuration change of the constituent particles. Our experiment also suggests the relation $\lambda \propto 1/(p-p_c) $ irrespective of the size of the container, where $p$ ($\propto fa$) is the momentum given by the external forcing, which is reminiscent of the de Broglie wave known in quantum mechanics. We further performed a numerical simulation on the onset and development of undulations using our simplified model, in which horizontal dilatation generated by mutual intrusion of particles near the bottom of the container develops to an arch formation of the granular layer. Influence of the friction of the granular material is examined.
\end{abstract}

\begin{altabstract}
Des études expérimentales et théoriques sont menées sur le rapport d'aspect entre la hauteur $\delta$ et la longueur d'onde $\lambda$ des ondulations générées dans une couche granulaire sèche confinée entre deux plans verticaux parallèles sous l'effet d'une vibration verticale de fréquence $f$ et d'amplitude $a$. Nous avons élucidé la propagation de l'onde de densité le long de la couche, ce qui suggère une constante élastique macroscopique inférieure d'un ou deux ordres de grandeur, attribuée à un changement de configuration des particules constitutives. Notre expérience suggère également la relation $\lambda \propto 1/(p-p_c) $, quelle que soit la taille du conteneur, où $p$ ($\propto fa$) est l'impulsion donnée par la force externe, qui rappelle l'onde de de Broglie connue en mécanique quantique. Nous avons ensuite réalisé une simulation numérique sur l'apparition et le développement des ondulations à l'aide de notre modèle simplifié, dans lequel la dilatation horizontale générée par l'intrusion mutuelle des particules près du fond du conteneur se développe pour former une arche dans la couche granulaire. L'influence du frottement du matériau granulaire est examinée.
\end{altabstract}

\keywords{\kwd{Vibrated granular layer} \kwd{undulation} \kwd{density wave propagation} \kwd{momentum transfer} \kwd{de Broglie wave}}
\altkeywords{\kwd{Couche granulaire vibrée} \kwd{ondulation} \kwd{propagation des ondes de densité} \kwd{transfert d'impulsion} \kwd{onde de de Broglie}}

\dateposted{2025-11-24}
\begin{document}

\maketitle

\vspace{4\baselineskip}

\section{Introduction}\label{introduction}

Since the pioneering work by Douady et al.~\cite{Douady}, lots of investigations on the undulation have been made. They performed experiments on vertically vibrating glass spheres, and examined the subharmonic instabilities, such as period-doubling instability, fragmentation of the layer, and kink-type topological defects, etc. They describe a regular pattern of defects above a critical excitation as ``the bouncing motion of the particle layer... occurs at twice the vibration period $T$''. They reported the dilated fluid layer part and solidified part that alternate spatio-temporally, and the minimum width of the arch created at the solidified part of the layer $l$ is characterized by $H/l \approx 0.16(\Gamma-4.2)$, where $H$ is the layer thickness, and $\Gamma=a(2\pi f)^2/g$ ($a$ and $f$ are the amplitude and frequency of external forcing, respectively, and $g$ is the acceleration of gravity). Similar spatio-temporal behavior of patterns which arise in nonlinear dissipative systems far from equilibrium by external forcing has been reported, such as transverse bending~\cite{Goldshtein}, arches~\cite{Lan}, and kinks~\cite{Moon}, etc. These pattern changes have common features of wavy deformations characterized by temporal arch formation of an almost constant thickness layer with integer or half-integer number of waves along the layer, and alternating ridge-foot positions with a period twice of the forcing period, which we call ``undulation'' in our investigation~\cite{Ugawa2003, Kanai, Sano2005, Sano2010, Sano2016}.

Quite a number of experimental or numerical studies have been made on the planar pattern formation of thin horizontal granular layer vibrated vertically, which elucidated regular wavy patterns, such as stripes, triangles, squares, hexagons~\cite{Melo1994, Melo1995, Metcalf, Bizon, Umbanhowar2000} as well as spatially localized excitations (oscillons)~\cite{Umbanhowar1996}, quasi-crystal patterns~\cite{Sano1999}. These wavy patterns are reminiscent of parametrically excited surface waves of fluid~\cite{Faraday1831}, and are successfully classified in terms of frequency $f$ and $\Gamma$. The pattern formation termed ``ripples'' in vertically vibrated granular material confined between two narrowly spaced parallel vertical plates~\cite{Clement, Aoki, Luding, Ugawa2002} also attracts attention, which focuses the vertical cross-sectional view of the above-mentioned planar patterns. A lot of dispersion relations are proposed depending on the layer height, friction, air viscosity, etc. reflecting the dissipative nature of the system. Meanwhile, a granular meltdown~\cite{Gotzendorfer} and the solid-fluid transition corresponding respectively to the compressed phase and dilated phase of the granular layer in a vertically vibrated granular material seem to add complexity on the characteristics of granular material~\cite{Sano2010, Sano2011, Sano2012, Sano2014}. In the latter, the presence of density wave has been reported in thicker granular layer, in which a ray theory of the refraction of density wave due to the vertically nonuniform distribution of grains is applied to determine the wavelength.

In this paper, we shall provide experimental evidences of the density wave in a thin layer of granular material confined in a container between two parallel vertical plates. The vertical momentum from the vibrating container (proportional to $fa$) is transported along the layer, which induces undulations (Section~\ref{secDensityWave}). In the latter, linear dependence of the aspect ratio $\delta/\lambda$ on $fa$ ($\delta$ and $\lambda$ are respectively the height and wavelength of the undulation arch), which has been shown in our previous paper (\cite[Figure~14]{Sano2016} and Section~\ref{linear}), is confirmed by the analogy to the kinetic theory of gasses (Section~\ref{kinematic}). Numerical simulation based on the dilatation due to the intrusion of the particles adjacent to the bottom wall into particles in the upper layers, thereby effect of the friction between particles is taken into account (Section~\ref{numerical}), agrees qualitatively with the experimental findings. Our study reveals a difference between material constants like Young's modulus and friction coefficient in a solid state and those associated with configuration change of granular material, which may provide a clue to deepen the physics of granular material.

\section{Experiment}\label{experiment}

Experimental apparatus is basically the same as our previous studies~\cite{Sano2010, Sano2011, Sano2012, Sano2014, Sano2016}. A rectangular container made of transparent acrylic resin with horizontal dimensions $L \times W$ and height $H$ was mounted vertically on an electromagnetic shaker. The container was oscillated sinusoidally with a frequency $f$ and an amplitude $a$, so that the position of the container bottom $z$ is specified as $z =z_0 + a \sin(2\pi ft+\phi)$. Here we have chosen the Cartesian coordinate axes $x$ and $z$ in the horizontal and vertical directions, respectively. The pattern formation of a vertically oscillated layer of dry granular material of grain size $d \approx \qty{1}{[mm]}$ filled in a container of horizontal dimension $L=90 \sim \qty{200}{[mm]}$ and width $W=3 \sim \qty{10}{[mm]}$ was observed from the side by a high-speed video camera with sampling time 500 to~2000 frames per second. All experiments were performed under atmospheric pressure. In our experimental setting, the height of the container $H$ is sufficiently large to prevent the escape of granular material, and $d \lesssim W \ll L$ assures the two-dimensional behavior of granular material between two narrowly spaced vertical walls. In a certain range of frequency, amplitude and phase of collision to the bottom wall, etc., stationary bending-wave-like structure (or undulation) of the granular layer, and hence characterized by eigen-modes oscillation, was recognized~\cite{Ugawa2003,Kanai, Sano2005, Sano2016}.

\begin{figure}[h!]
\includegraphics[height=50mm]{UndSesGlass}
\hfill
\includegraphics[height=50mm]{UndLead}

\vskip .25\baselineskip

\includegraphics[width=\textwidth]{TimeSequence}

\vskip .2\baselineskip

\includegraphics[width=.48\textwidth]{Undulationf33a17}
\hfill
\includegraphics[width=.48\textwidth]{f2Ripplef14a49}
\caption{Characteristics of undulations. (a)~Sesame seeds (vibrated at $f=\qty{30}{[Hz]}$, $a=\qty{2.02}{[mm]}$ in a container of $L=\qty{91}{[mm]}$), and (b)~glass spheres of $d=\qty{1}{[mm]}$ (vibrated at $f=\qty{30}{[Hz]}$, $a=\qty{3.40}{[mm]}$ in a container of $L=\qty{140}{[mm]}$). Two snapshots in respective Figures~(a) and~(b) are taken at time interval $T$ apart, where $T$ is the period of external forcing. (c)~Undulations observed for lead spheres of $d=\qty{1}{[mm]}$ in different size of the container $L$ (shown at respective pictures) due to the external forcing of almost the same $fa$. (d)~Time sequence of undulations of lead spheres ($d=\qty{1}{[mm]}$) observed in a section (width of strip \qty{4.4}{[mm]}) fixed in space ($f=\qty{49}{[Hz]}$, $a=\qty{1.04}{[mm]}$). (e)~Time sequence of the lower boundary of the layer viewed at its nodal position of undulations ($f=\qty{33}{[Hz]}$, $a=\qty{1.7}{[mm]}$) and (f)~that of $f/2$-ripples ($f=\qty{14}{[Hz]}$, $a=\qty{4.9}{[mm]}$), where dotted curves show the position of the container bottom. In the former, collision of the layer with the bottom wall is ``soft'', while it is ``hard'' in the latter.}\label{universal}
\end{figure}

\goodbreak

\begingroup
From our previous experiments, we have revealed the following universal features of undulations of dry granular layer under external forcing of the container.
\begin{enumerate}[leftmargin=*]
\item Thin layers of dry granular material show wavy deformation with period $2/f$ in a certain range of $f$ and $a$, irrespective of the kinds of granular material and the size of container. Figure~\ref{universal} shows examples of undulations observed in (a)~sesame seeds, (b)~glass spheres, and (c)~lead spheres in the container of different horizontal sizes $L$. Two snapshots in respective figures~(a) and~(b) are taken at time interval $T$ ($=1/f$) apart, where $T$ is the period of external forcing. Figure~\ref{universal}(d) shows the time sequence observed in a section fixed in space, which reveals period $2T$ behavior, (e)~time sequence of the bottom wall and the lower boundary of the layer at a fixed horizontal position of undulations, and (f)~that of ripples.
\item Lower boundaries of the granular layer show clear arches irrespective of the layer thickness and materials, while the deformation of upper boundaries is blurred as the thickness increases (see Figure~\ref{universal}(d)). In sufficiently thicker layers, upper boundaries are nearly flat.
\item Stable undulations are maintained when the ends of the layer meet perpendicularly to the vertical side walls of the container. We denote the mode $S_n$ if it satisfies $n \lambda=L$, and $A_n$ if $(n-1/2)\lambda =L$, where $n$ is a natural number and $\lambda$ is the wavelength. Figure~\ref{universal}(c) shows undulations observed for lead spheres of $d=\qty{1}{[mm]}$ in the different lengths of containers $L$ with almost the same $fa$, which exhibits almost the same wavelength.
\item Stationary undulations are maintained when the lower boundary of the layer {viewed at its nodal region} meets the downward motion of the container bottom, followed by pushing up with the upward motion of the latter (sometimes called ``soft collision'', see Figure~\ref{universal}(e)). On the contrary, if the layer in the downward motion meets the bottom wall at (nearly) the strongest upward motion of the container, it develops to ``ripples'' (sometimes called ``hard collision'', see Figure~\ref{universal}(f)).
\item When the wavelength does not exactly satisfy the matching condition mentioned above, the wavelength is modified by forming a solidified part along the layer, which effectively plays a role of boundary wall, or sometimes creating sub-arches to adjust the main structure of undulations~\cite{Sano2016}.
\end{enumerate}
These characteristics are obtained mainly by visualization. In the following, we show the dynamical aspects of undulations.
\endgroup

\begin{figure}[h]
\includegraphics[width=80mm]{UndulationExp2}
\caption{Observation of undulations in $f$-$a$ plane. Data point of the experiment by Douady et al.~(estimated from reference~\cite{Douady}) is also added, although the details of experimental conditions may not be the same.}\label{f-a-diagram}
\end{figure}

\subsection{Observed region of undulations}

Figure~\ref{f-a-diagram} is an $f$-$a$ diagram of the observed region of undulations. Results for lead, aluminum and glass spheres are superposed. Here, differences of the layer thicknesses, and the sizes of the container, etc.\ are disregarded except for the condition on the geometry of the container $3\lesssim W/d \lesssim 10$ that assures the two-dimensionality in the vertical cross section. Generally speaking, observed regions of undulations in the former two materials overlap indistinguishably, while that of the glass spheres seems to shift to larger $fa$ area. For comparing quantitative characteristics such as wavelength and wave height, etc., however, we employ the data obtained in the container of similar sizes, and lower mode undulations such as $A_1, S_1, A_2$ are eliminated (otherwise stated) to avoid the influence of the side walls.

\subsection{Container size dependence of the wavelength of the eigen-mode}

As has been mentioned, undulations are interpreted as eigen-modes of buckling of the granular layer. This naturally raises a question on the dependence of the container size on the wavelength. In order to clarify the latter, we performed experiments using a container of varied horizontal sizes, and compared the wavelengths of undulation. Figure~\ref{universal}(c) is an example of snapshots showing the dependence of the container size $L$ on the undulation wavelength $\lambda$. These waves are observed under the conditions that are close to eigen-modes in a given container. We can recognize the agreement between $\lambda$ of $S_2$ mode in $L=\qty{100}{[mm]}$ (Figure~\ref{universal}(C1)) and $\lambda$ of $S_4$ mode in $L=\qty{200}{[mm]}$ (Figure~\ref{universal}(C4)) within the experimental accuracy. In view of the exact eigen-value problem of the buckling of a continuum elastic plate, the container size $L$ required for $A_3$ is \qty{125}{[mm]}, which looks nearly satisfied in the $L=\qty{120}{[mm]}$ case (Figure~\ref{universal}(C2)). Similarly, $S_3$ mode requires $L=\qty{150}{[mm]}$, which is marginally satisfied in the $L=\qty{140}{[mm]}$ case (Figure~\ref{universal}(C3)). The latter two cases suggest that the condition on the eigen-value problem of the granular layer is relaxed by creating a solidified part within the granular material which effectively plays the role of temporarily built wall of the container, so that the relation $n=L/\lambda$ in $S_n$ case (or similarly $n=L/\lambda+1/2$ in $A_n$ case) is approximately satisfied. The adjustment of the wavelength and the transition from one mode to the other are fully described in our previous paper~\cite{Sano2016}. These experiments suggest that the wavelength is determined by the product $fa$ irrespective of the container size provided that near eigen-value conditions are satisfied.

\subsection{Aspect ratio of arch height to wavelength}\label{linear}

In this subsection, we shall pay attention to the aspect ratio of the height $\delta$ to wavelength $\lambda$ of undulation arch (see the definition sketch Figure~\ref{universal}(d)). Figure~\ref{AspectRatio}(a) shows the relation between $\delta / \lambda$ and non-dimensional acceleration $\Gamma$ ($\equiv 4\pi^2 f^2 a/g$), while Figure~\ref{AspectRatio}(b) shows the one between $\delta / \lambda$ and $fa$, measured in the containers of different sizes ($L, W$) and external forcings ($f, a$). The number of the layer $N$, which is specified in the initial state at rest, gives an approximate height of the layer $h$ ($\approx \bracks[\big]{1+(N-1)\sqrt{3}/2}d$). In accordance with our previous work~\cite{Kanai}, effect of the horizontal depth $W$ may not be negligibly small for $W \lesssim 3d$, while the layer of height $h \gtrsim 10d$ may suffer splitting of layers between the group of particles near bottom region and that in upper region depending on the magnitude of vertical acceleration~\cite{Sano2014}, which restricts our attention to $W\gtrsim 4d$ and $4d \lesssim h \lesssim12d$. Figure~\ref{AspectRatio}(a) shows a relation between $\delta/\lambda$ and the non-dimensional acceleration $\Gamma$, which exhibits a dispersion in a broad range. On the other hand, the same data points are well-fitted by a single line
\begin{equation}\label{eq_2_1}
\frac{\delta}{\lambda}=\alpha \bracks[\big]{fa -(fa)_c},
\end{equation}
where $\alpha = \qty{0.00316}{[s/mm]}$ and $(fa)_c=\qty{16.4}{[mm/s]}$ for the case of lead spheres of $d=\qty{1}{[mm]}$ as shown in Figure~\ref{AspectRatio}(b). Note that we reported $\alpha=\qty{0.00365}{[s/mm]}$ and $(fa)_c=\qty{20.5}{[mm/s]}$ for the lead sphere case in our previous paper~\cite{Sano2016}, the difference of which is attributed to the choice of data points. In Figure~\ref{f-a-diagram}, observed undulations are all included, while we have chosen undulations that are close to eigen-modes in the present case. Furthermore, to avoid the uncertainty ascribed to the influence of the side-walls of the container, only undulations that are close to eigen-mode with several waves are chosen. In addition, the waves adjacent to the side-walls are not employed for the measurement of $\delta$ and $\lambda$. The data obtained for the case of glass and aluminum beads of the same sizes, as well as the one for the sesame seeds of similar sizes (volume-equivalent size is \qty{1.7}{[mm]}), are also shown, which seems to support the linear dependence on $fa$ like eq.~\eqref{eq_2_1}. A closer look of the data shows that $\alpha$ and $(fa)_c$ for lead spheres and aluminum spheres are almost indistinguishable to each other. In fact, our sand-pile measurement using the same materials gives the angle of repose $\theta_r = \ang{29}$ and $\ang{27}$ for lead spheres and aluminum spheres, respectively. Considering the ambiguity of determining the boundary surface of granular material, variation of above-mentioned angle $\Delta \theta = \ang{2}$, and hence the difference of $(fa)_c$ between 16.4 and 19.9 shown in the parentheses may be within an experimental error. On the other hand, the data for the glass beads seems to be considerably separated from the former two materials, in spite of having similar angle of repose. This point will be discussed in the later sections. A remarkably good fitting with a single line implies that the structure of undulation is not remarkably influenced by the sizes of the container $L,\, W$ and the total number of particles provided that the quasi-two dimensionality in the vertical plane is maintained.

\begin{figure}[t]
\includegraphics[width=.48\textwidth]{Gam-dependenceLASG}
\hfill
\includegraphics[width=.48\textwidth]{fa-dependenceLASG}\
\caption{Container-size dependence of undulation arch of the thin layer of lead, aluminum, glass spheres of $d=\qty{1}{[mm]}$, and sesame seeds of volume-equivalent size $d=\qty{1.7}{[mm]}$. (a)~Relation between $\lambda / \delta$ and $\Gamma$. (b)~Relation between $\lambda / \delta$ and $fa$, where dotted lines and the critical values in the parentheses are least square fitting of the experimental data of respective materials obtained in the container of $L=\qty{140}{[mm]}$. To avoid the uncertainty ascribed to the influence of the side-walls of the container, only the data obtained in near eigen-mode undulations with several waves are shown.}\label{AspectRatio}
\end{figure}

The differences of pattern characterization with $\Gamma$ or $fa$ may come from the different physics between the 2D planar patterns observed in a thin vibrated granular layer and the vibration induced patterns of granular layer in thin vertical 2D space. The former is characterized by fluid-like behavior similar to Faraday waves on the surface of water, while the latter is characterized by bending of solid layer, which we shall discuss in later sections.

Equation~\eqref{eq_2_1} is rewritten as
\begin{equation}\label{eq_2_2}
\lambda = \frac{\delta}{\alpha\parens[\big]{fa -(fa)_c}},
\end{equation}
which is reminiscent of the de Broglie wave $\lambda=h/p$ in quantum mechanics ($h$ is the Planck's constant, and $p$ is the momentum)~\cite{deBroglie} in that the wavelength is inversely proportional to the momentum. Of course, our granular material is a dissipative system which reflects the introduction of threshold value $(fa)_c$, in addition to the extremely large difference in spatial and temporal scales between the two. We shall consider the above-mentioned analogy between microscopic and macroscopic worlds in the appendix.

\subsection{Density wave propagation and undulations in a vibrating thin granular layer}\label{secDensityWave}

Ten-layer glass beads of a diameter $d = \qty{1.0}{[mm]}$ vibrated vertically with frequency $f = \qty{40}{[Hz]}$ (hence the period of oscillation of the container $T$ is \qty{25}{[ms]}) and amplitude $a = \qty{1.2}{[mm]}$ in a container $W/d \approx 1$ was observed by a high-speed video camera with sampling time 1/1000\,\unit{[s]}. The pattern period of undulations in the present case is $2T$ ($=\qty{50}{[ms]}$). We show the density wave propagation along the vertically vibrated undulating layer in Figure~\ref{DensityWaveGlass}. In each figure, the baseline is chosen at the bottom of the container, and two successive images of particles with time separation \qty{1}{[ms]} are superposed. Contours of particles at one shot are shown by red open circles, on which those \qty{1}{[ms]} later are superposed by blue open circles. Thus, the particles with no relative motion exhibit clear blue circles, while those with relative motion exhibit blurred circles in red and blue. A closer look at these figures shows the presence of regions with no relative motion, which is regarded as ``solid'' domain, while the regions with relative motion are regarded as ``fluid'' domain. The two domains change sizes and positions along the undulating layer.

In Figure~\ref{DensityWaveGlass}(a) to~(i), one half period of undulation pattern is shown, where the numerical number put on respective figures is the time in millisecond. The container is in the upward sinusoidal motion from~(a) to~(e), while it is in downward sinusoidal motion from~(e) to~(i). Arrows at the lower left of respective figures schematically show the motion of the container, and the double circles in (a)(e) and (i) show that the container is in the equilibrium position. In the first frame~(a), particles in the lowest part of the layer touch the bottom wall of the container. In~(b), particles in the central small region adjacent to the bottom wall are pushed up and solidified, while the rest of the layer continues falling under gravity and looks loosely connected. Here, arcs in black show the schematic boundary between the solidified and the fluidized domains. In~(c), the compressed region develops to the other (upper) side of the layer, which further extends along the layer as described in~(d). At the end of the upward motion of the container~(e), almost all particles in the lowest layer touch the bottom wall, but those falling down (particles placed at both sides in~(a)) undergo stronger collision, so that the latter is compressed and solidified~(f). Meanwhile, the collision of the layer with the container {bottom} causes dilatation, which is resolved by the undulations of the next cycle~(g). As the container moves downward, the layer itself is in almost free flight motion, but {inward motion of the compressed region is recognized, which is exhibited as} the compression wave toward the central region ((g) and (h)). Figure~\ref{DensityWaveGlass}(i) is the configuration of half-period after~(a). Here, the central mound is solidified, which corresponds to the solid domain previously present at both sides in~(a). It may also be interpreted that the configuration of~(i) is the one spatially shifting~(a) by a half wavelength.

\begin{figure}[h]
\includegraphics[width=\textwidth]{DensityWaveGlass}
\caption{Density wave propagation in undulating monolayer ($W/d\approx 1$) of glass beads oscillating with frequency $f = \qty{40}{[Hz]}$ and amplitude $a = \qty{1.2}{[mm]}$. Positions of the particles relative to the container at one moment are shown in red open circles, while those \qty{1}{[ms]} later are shown in blue open circles. In each figure, the latter is superposed to the former. The clear blue circles show that they are immobilized with respect to the container, and the domain of such connected circles imply the solidified part.}\label{DensityWaveGlass}
\end{figure}

The previous observation was made on the monolayer of glass beads ($W/d \approx 1$), aiming at accurately capturing the position of the particle. To add an evidence of the density wave propagation along the undulating layer, ten-layer lead beads of diameter $d = \qty{1.0}{[mm]}$ in the container of $W/d = 5$ vibrated vertically with frequency $f = \qty{30.0}{[Hz]}$ and amplitude $a = \qty{2.10}{[mm]}$ were observed by a high-speed video camera with sampling time 1/2000\,\unit{[s]} (Figure~\ref{DensityWaveLead}). Every two consecutive images $I(\mathrm{frame number})$ were subtracted (and hence $\Delta I(2n)=I(2n)-I(2n-1)$, $n=1,2,\dotsc$), where common baseline of the images was chosen at the bottom wall of the container. The pattern period of undulations in this case is $2T$ ($= 1/15\,\unit{[s]}$), which corresponds to about 133.3 frames of the original images, so that Figure~\ref{DensityWaveLead} covers about $1/4$ pattern period. In each figure, points in white show the particles with relative motion, while points in black show either (1)~the void space or (2)~particles with no relative displacement during $\Delta t =1/2000\,\unit{[s]}$. The black domains {between lower boundary of the layer and the bottom wall that correspond to arches (hence denoted by a letter ``a'') are examples of~(1), besides the trivial void space above the upper surface of the layer. On the other hand}, those domains denoted ``C'' describe particles so tightly compressed to each other that no relative motion is allowed, which is an example of~(2). Note that our granular material consists of dry non-cohesive spheres, so that the spheres in ``C'' compose a temporally solidified domain. As shown in frame 78 on, the domain C initially a small part in contact with the bottom wall develops upward to the other side of the layer at frame 94. The latter domain grows horizontally in both directions along the layer until frame 98, whose propagation speed $V$ is about \qty{40}{[m/s]} in the present example. Then, the neighboring density waves of opposite directions collide at the middle part of the compacted domain, and exhibits burst-like behavior near the upper part of the layer (frames 98--100). After this stage, the compacted domain gradually loosens or dilates (as denoted ``L''), and becomes almost flat with more or less uniform distribution of particles.

\begin{figure}[ht]
\includegraphics[width=90mm]{UndDensityWave}\
\caption{Density wave propagation in undulating layer of lead spheres. {High-speed video image of sampling time 1/2000\,\unit{[s]} relative to the container is obtained. Each sampled frame is given by subtracting the previous frame. The black domain marked ``C'' shows the part that is immobilized with respect to the container, which implies the solidified part. The domain with} mark ``L'' shows the loosened (dilated) domain, and ``a'' shows the void space of arch ($f=\qty{30}{[Hz]}$, $a=\qty{2.1}{[mm]}$).}\label{DensityWaveLead}
\end{figure}

According to the theory of elasticity, the propagation velocity of longitudinal wave $V_{\solid}$ in 1D elastic material of density $\rho$ and Young's modulus $Y$ is given by $V_{\solid}=\sqrt{Y/\rho}$, which is about \qty{4000}{[m/s]} in solid state lead. Compared with the density wave propagation in an ordinary solid material, above-mentioned wave speed in a granular material is very slow, which implies that the granular material is regarded as an extraordinarily soft material. Measurement of the density wave may provide a fundamental data on dynamical properties, such as effective Young's modulus and effective Poisson's ratio, which will provide a clue to develop a new physics taking account of the configuration change of the constituent particles. A role of density wave in a granular layer to the arch formation will be considered in the later section (Section~\ref{kinematic}).

\section{Theory on undulations of vertically vibrated thin granular layer} \label{theory}

\subsection{Configuration change and dilatancy}

Throughout our experiments on undulation, we observed that the lowest boundary of the layer is rather clearly observed, while that of the upper layer becomes less clear. Mutual intrusion of the grains in the vertically vibrating container is limited to a few layers adjacent to the bottom wall, which forms a stronger stress network and contributes to the dilatancy. On the other hand, grains in the upper layers look simply piled up by the excluded volume effect~\cite{Duran} (see Figure~\ref{universal}), which explains rather insensitive dependence of the layer number $N$ ($\gtrsim 10$) on wavelength $\lambda$. These evidences lead us to regard the thin granular layer as a thin elastic plate~\cite{Sano2005, Sano2016}, which, however, needs some caution in the granular layer, because (1)~elastic constants such as Young's modulus will be associated with the configuration change of particles, (2)~constituent particles, here assumed uniform spheres of unit radius $r$ ($=1$), are almost in contact with each other in the granular material, so that their structure is maintained by the forces normal to their contact plane to retain the excluded volume effect, and (3)~in contrast to the elastic deformation (which may be energy conserving), particles will slide each other, which is associated with energy dissipation due to friction, etc.

In spite of these complexities, we first examine the formation of arch structure taking account of the analogy to the deformation of elastic bodies (Section~\ref{elastic plate}). Then, we shall consider a numerical model to understand the essential processes on the formation of undulations (Section~\ref{numerical}). In the latter, we focus our attention mainly up to three layers of particles adjacent to the bottom wall, which was shown to be the minimum number of layers that is dynamically important to create ripples and/or undulations~\cite{Kanai}. Spheres in the middle layer are squeezed between the bottom wall and the upper layers. For sufficiently strong external forcing, spheres in the lowest layer will intrude vertically into the middle ones, and the resulting dilatation in the horizontal direction will necessitate the bending of the layer to compromise the presence of the side walls (Figure~\ref{dilatation}).

\begin{figure}[h!]
\includegraphics[height=70mm]{UndArch}
\caption{Mutual intrusion of particles, dilatancy and arch formation (schematic).}\label{dilatation}
\end{figure}

\subsection{Buckling of a thin elastic plate}\label{elastic plate}

In general, granular material is in a quasi-equilibrium state under the force balance between external force, normal force supported by the excluded volume effect and the tangential force due to the frictional force. If the particles are crammed and frictional force is large enough, the configuration of the constituent particles will be unchanged, so that the granular material can be regarded as a continuum. Let us assume our thin granular layer as a thin elastic plate with horizontal length $L$, width $W$ and height $h$, and consider its deformation due to the longitudinal force originating in the dilatancy. We introduce the Cartesian coordinate system, with the $x$ and $z$ axes in the horizontal and vertical directions, respectively, assuming a constant width $W$ in the third direction~($y$ axis). In the case $h, W \ll L$, position $P$ in the layer (or arc) is well specified as $(x,z)$, which is also described by the length $s$ along the layer from the origin and the angle $\theta$ tangent to the arc at that position.

According to the theory of elasticity (see e.g.,~\cite{LL7}), deformation of the elastic plate mentioned above is described by the equation
\begin{equation}\label{eq_3_1}
\frac{\dif^2 \theta}{\dif s^2}= -\beta^2 \sin \theta,
\end{equation}
where $\beta = \sqrt{F_x/YI}$. Here, $F_x$ is the force applied horizontally at both ends of the layer, $Y$ is the Young's modulus of the material and $I$ ($\equiv Wh^3/12$) is the geometric moment of inertia of the cross-sectional area. After some calculation~\cite{Sano2005}, we obtain
\begin{equation}\label{eq_3_2}
x=\frac{2}{\beta}\bracks[\Bigg]{Z(\beta s, k)+\parens[\Bigg]{\frac{E}{K}-\frac{1}{2}}\, \beta s},
\qquad
z=\frac{2k}{\beta}\parens[\big]{1-cn(\beta s, k)},
\end{equation}
where $Z(\phi,k)$ is the Jacobi's zeta elliptic function, $cn(\beta s, k)$ is the Jacobi's elliptic function, $K$ ($=K(k)$) and $E$ ($=E(k)$) are the complete elliptic integrals of the first and second kind, respectively, and $k$ ($0\le k <1$) is called the elliptic modulus. As will be shown later, $k$ implies the amount of aspect ratio of undulation (see the insets of Figure~\ref{ElasticPlate}). At this moment, we have imposed the boundary condition $z = 0$ at $s=0$. In terms of the $(x,z)$ coordinate system, the boundary conditions on the granular layer require the perpendicular contact to the side walls at $x=0$ ($s=0$) and $x=L/M^*$ (i.e., $\beta s =2M^* K(k)$), where odd and even numbers of $M^*$ correspond to the antisymmetric modes $A_n$ and symmetric modes $S_n$, respectively. Considering the periodicity of $Z(\phi,k)$ function with period $2K$, the above boundary conditions yield
\begin{equation}\label{eq_3_3}
\beta L =2M^*\bracks[\big]{2E(k)-K(k)}.
\end{equation}
For $k \ll 1$, eq.~\eqref{eq_3_3} becomes
\begin{equation}\label{eq_3_4}
\beta L =\pi M^* \parens[\Bigg]{1 -\frac{3}{4} k^2 +O(k^4)},
\qquad \text{or} \qquad
k^2 \approx \frac{4}{3}\parens[\Bigg]{1 -\frac{\beta L}{\pi M^*}}.
\end{equation}
As far as the small $k$ is concerned, which is decisive for the onset of buckling, the condition that $k$ has a real number is $ \beta \le \pi M^* / L $. Then, the critical value of $\beta_c$ is given for the lowest mode $M^*=1$ (i.e., $A_1$ mode), from which the critical force to generate undulation is
\begin{equation}\label{eq_3_5}
{(F_x)}_c = \beta_c^2 YI=\frac{\pi^2}{L^2}YI=\frac{\pi^2 Wh^3}{12 L^2}Y,
\qquad \text{or} \qquad
{(f'_x)}_c \equiv \frac{{(F_x)}_c}{Wh}
=\frac{\pi^2}{12}\parens[\Bigg]{\frac{h}{L}}^2 Y.
\end{equation}
Above this critical external forcing, another branch of solution emerges. By applying the boundary condition that the layer meets perpendicularly to the side walls of the container, the solution of eq.~\eqref{eq_3_1} admits particular eigen-modes of undulation, whose behavior is described in terms of the elliptic functions~\cite{Sano2005}. Based on that theory, the aspect ratio of the wave height $\delta$ to wavelength $\lambda$ is
\begin{equation}\label{eq_3_6}
\frac{\delta}{\lambda}=\frac{2k}{2E(k)-K(k)} \quad
\, \bracks[\Big]{\,=\frac{4}{\pi}\, k\, \parens[\Bigg]{1+\frac{3}{4}k^2+\dotsb}
\quad \text{for} \quad
k \ll 1\,},
\end{equation}
where $k$ is determined by adjusting the boundary condition at the peak or foot of the undulations of the layer (see Figure~\ref{ElasticPlate}).

When the initial flat layer of length $L$ is undulated to form an arch, the amount of elongated layer $\Delta L$ is calculated to be
\begin{equation}\label{eq_3_7}
\frac{\Delta L}{L}=\frac{2\bracks[\big]{K(k)-E(k)}}{2E(k)-K(k)}\equiv \Delta,
\end{equation}
from which we obtain (see also Figure~\ref{ElasticPlate})
\begin{equation}\label{eq_3_8}
\Delta\approx k^2 \parens[\Bigg]{1+\frac{9}{8} k^2+\dotsb}
\quad \text{or} \quad k \approx\sqrt{\Delta}\, \parens[\Big]{1-\frac{9}{16}\Delta+\dotsb},
\quad \text{for} \quad k \ll 1.
\end{equation}
Given the information on $\Delta L/L$, and hence on $k$, aspect ratio of undulation $\delta/\lambda$ will be determined by eq.~\eqref{eq_3_6}:
\begin{equation}\label{eq_3_9}
\frac{\delta}{\lambda} \approx \frac{4}{\pi}\sqrt{\Delta}\, \parens[\Bigg]{1+\frac{3}{16}\Delta+\dotsb}.
\end{equation}

\begin{figure}[h]
\includegraphics[height=50mm]{AR-ElasticPlate3}\quad
\includegraphics[height=50mm]{dellamb-Delta}
\caption{Aspect ratio $\delta/\lambda$ and elongation of an elastic plate $\Delta L/L$. Insets show the image of the undulation with $k=0.1$ and $k=0.2$.}\label{ElasticPlate}
\end{figure}

\subsection{Stationary state undulations in the granular layer}\label{kinematic}

As has been shown in our previous section, vertical momentum applied to the granular layer by the vibration of container $z=a\parens[\big]{1-\cos(2\pi f t)}$ is transmitted to the side walls as a density wave if $fa$ exceeds a certain critical value $(fa)_c$. In the subcritical case, relative motion of spheres will not be induced because of larger gravity force and larger static friction between spheres in contact. As $fa$ exceeds $(fa)_c$, however, a certain period of oscillation appears that allows configuration change among spheres, whose extent will depend on $fa-(fa)_c$. Then, the velocity of the density wave $V$ is a function of $fa-(fa)_c$, and will be described as
\begin{equation}\label{eq_3_10}
V = 2\pi C_1 \parens[\big]{fa - (fa)_c}+ o\parens[\big]{fa - (fa)_c},
\end{equation}
near $fa = (fa)_c $, where $C_1$ is a constant. Once constituent spheres are fully mobilized, they will create a ``dynamical pressure''. In the stationary state undulations, a quasi-equilibrium state will be maintained between the balance of energy input by the external forcing and the energy dissipation due to a certain effective dynamic friction. To avoid complexities, however, we may apply the analogy to the kinetic theory of gasses, and estimate the pressure on the side walls (or the nodes of undulation waves) in the horizontal direction in the macroscopic sense. We consider a test sphere of mass $m$ moving with velocity $\upsilon_0$ toward side walls of the container (or nodes of waves) of length $L^* \equiv L/M^*$, where $M^*$ characterizes the eigen-mode defined by eq.~\eqref{eq_3_3}. The particle will collide the wall with momentum $m\upsilon_0$ and reflected with velocity $-e\upsilon_0$, where $e$ is the restitution coefficient and $\upsilon_0 \sim 2\pi \parens[\big]{fa - (fa)_c}$, so that the side wall undergoes a momentum $(1+e)m\upsilon_0$ at each collisions. The particles considered here are more or less connected to each other through stress network, so that the momentum transfer will be effectively repeated on both sides of the container. Note that our typical experiment (see e.g.~Figure~\ref{DensityWaveLead}) shows $V\sim \qty{40}{[m/s]}$, while $\upsilon_0 = 0.3\sim \qty{0.6}{[m/s]}$, which implies that the number of collisions $n$ to one of the walls, estimated to be $n \sim V/(2L^*)$, is sufficiently large. Then, the impulse that one particle exerts to the side wall will be $(C_2 m\upsilon_0) n $ per unit time, where $C_2$ is a certain constant of order unity reflecting the structure of the stress network. Accordingly, total number $N_{\tot}$ of the particles in the granular layer exerts a force $F_x$ on the side wall, which yields the dynamical pressure
\begin{equation}\label{eq_3_11}
\begin{split}
p
& = \frac{F_x}{h^*W}
\\	& = \frac{C_2}{h^*W}m\upsilon_0 \parens[\Bigg]{\frac{V}{2L^*}}N_{\tot}
\\	& = \frac{1}{2} C_2 \parens[\Bigg]{\frac{mN_{\tot}}{L^*h^*W}}\upsilon_0 V
\\	& = \frac{1}{2} C_2 M^* \rho^* \upsilon_0 V
\\	& = 2\pi^2 C_1 C_2 M^* \rho^* \parens[\big]{fa - (fa)_c}^2+\dotsb,
\end{split}
\end{equation}
where we estimate the effective density $\rho^* = mN_{\tot} / \bracks[\big]{L\times (h/2)\times W}$ by assuming effective layer height $h^*\approx h/2$. As a consequence, we have
\begin{equation}\label{eq_3_12}
p = \frac{F_x}{S_{yz}} =Y^* \frac{\Delta L}{L}= C_p \rho^* \parens[\big]{(fa-(fa)_c)}^2,
\end{equation}
where $C_p=2\pi^2 M^* C_1\, C_2$ is a constant and $Y^*$ is the effective Young's modulus taking account of the configuration change of particles. From eqs.~\eqref{eq_3_6}--\eqref{eq_3_8} and~\eqref{eq_3_12}, we have
\begin{equation}\label{eq_3_13}
\frac{\delta}{\lambda} \approx \frac{4}{\pi}\sqrt{\frac{\Delta L}{L}}
\approx \frac{4}{\pi}\sqrt{\frac{C_p \rho^*}{Y^*}}\, \parens[\big]{fa-(fa)_c},
\end{equation}
which explains the linear dependence of $\delta /\lambda$ on $fa-(fa)_c$ shown in Figure~\ref{AspectRatio}(b).

Our granular spheres in undulation is nearly in closest packing, which implies $\rho^* \approx 0.7 \rho$. On the other hand, we have no available data for $Y^*$. However, we have some knowledge on the velocity of density wave (see Section~\ref{secDensityWave}). The one shown in Figure~\ref{DensityWaveLead} is an example of an $S_2$ mode undulation of lead spheres with $f=\qty{30.0}{[Hz]}$ and $a=\qty{2.10}{[mm]}$ ($fa=\qty{63}{[mm/s]}$), where observed aspect ratio is $\delta/\lambda=0.14$ and the velocity of density wave $V \sim \sqrt{Y^*/\rho^*}$ is about \qty{40}{[m/s]}. Then, eq.~\eqref{eq_3_10} is reproduced by putting $C_1\approx100$. If we put $C_2=2$ (as is for an elastic collision), we have $C_p\approx 4000$, and
\begin{equation}\label{eq_3_14}
\frac{\delta}{\lambda} ={\alpha^*} (fa-(fa)_c),
\quad \alpha^*=\qty{0.002}{[s/mm]},
\end{equation}
which is a reasonable agreement with our experimental results, although $(fa)_c$ is not determined in the above treatment.

\subsection{Onset of undulations and effect of friction}\label{numerical}

In this subsection we shall consider the effect of friction on the configuration change of spheres due to the vibration of container $z=a\parens[\big]{1-\cos(2\pi f t)}$. In addition to the frame of reference $(x,z)$ ($S$-system) fixed to the observer, we introduce a frame of reference $(x', z')$ ($S'$-system) moving with a marked sphere~$A$, as well as the polar coordinate system $(R,\theta)$ with its origin at the center of sphere~$A$ (Figure~\ref{Intrusion}(a)). Here, the angle $\theta$ is measured in anti-clockwise direction. Figures~\ref{Intrusion}(a) and~(b) depict the relative motion of spheres near the bottom wall of the container that is in upward acceleration, i.e.~$t=t_0 \sim t_1$, where $t_0$ ($\approx 0$) is the time when the sphere~$A$ in the lowest layer touches the wall and $t_1$ is the time it is released from the wall (see also Figure~\ref{universal}(e)). The sphere~$A$ is forced to move upward with the velocity $\upsilon(t)=2\pi f a \sin(2\pi f t)$ during the period $t=t_0 \sim t_1$. In the upward phase, sphere~$A$ is in contact with the bottom wall, so that $z'=z-a\parens[\big]{1-\cos(2\pi f t)}$, and the acceleration of inertial force $\alpha(t)=(2\pi f)^2 a \cos(2\pi f t)$ is added to the acceleration of gravity $g$ in the $S'$-system.

\begin{figure}[h!]
\includegraphics[height=55mm]{Intrusion-3m}

\includegraphics[height=48mm]{OrbitRM1mu1-faLarge}
~
\includegraphics[height=48mm]{ContactPointR}

\includegraphics[height=48mm]{OrbitR_fa10-mu}
~
\includegraphics[height=48mm]{OrbitR-fa10-m}
\caption{Configuration change by intrusion. (a)--(b)~Schematics of the configuration of spheres. Development of the angular position of contact point in (c)~larger $fa$ case and (d)~smaller $fa$ case, where dotted lines show the one after the upward acceleration phase of the container (after the time $t=T/4$). Comparison of the movement of contact point: (e)~dependence on the friction constant $\mu$ under given $m^*$ ($=1$); (f)~dependence on the effective mass $m^*$ under given $\mu$ ($=1$). One example of the combination of $f=40$ and $a=0.25$ ($fa\, =10$) is shown in~(e) and~(f).}\label{Intrusion}
\end{figure}

Figure~\ref{Intrusion}(a) describes an example sketch of the relative motion of the spheres around a sphere~$A$. In general, spheres may distribute more irregularly. From our experiment, however, spheres are observed to take nearly regular hexagonal positions near the bottom wall of the container (as shown by dotted circles in Figure~\ref{Intrusion}(a)). At the impingement of the lowest-positioned sphere on the container bottom, sphere~$A$ will intrude into the second row of spheres in the upper layer (spheres~$B$ and~$D$), while sphere~$C$ (right above~$A$) is likely to move vertically in contact with the sphere~$A$. As far as the very early stage of the impingement is concerned, other spheres change positions according as they contact with about three points among them as illustrated in Figure~\ref{Intrusion}(a). At this moment on, friction between the spheres is regarded as the kinematic one that is characterized by the dynamic friction coefficient $\mu$. This process induces dilatation in the horizontal direction, creates a certain space among all spheres and settles in a new quasi-equilibrium configuration. The amount of the latter adjustment seems to attenuate quickly toward upper layers~\cite{Sano2010, Sano2011}, which explains the rather insensitive movement of grains near the upper layers (see Figure~\ref{universal}). These experimental findings suggest the importance of the dynamical processes in lower few layers near the bottom wall. In the nearly condensed system of spheres, spinning motion of respective spheres will be suppressed by ``frustration''. Centrifugal force may also be negligible because of the small free space for rotational or circulatory motion. Consequently, the sphere~$B$ will move in contact with~$A$, and the contact point~$P$ is well described by the angle $\theta$. Without information on detailed structure of the stress network among spheres and with an analogy to the fluid mechanics, we choose the sphere~$B$ as a test particle and introduce an effective mass ${\tilde m} \equiv m+\Delta m$ to include the averaged effect from other spheres, where $\Delta m$ is the added mass that any body in accelerating motion undergoes in the continuum. The validity of our approach will be discussed in the last section. Then, the equation of motion of sphere~$B$ in the tangential direction in the $S'$-system is given as
\begin{equation}\label{eq_3_15}
2{\tilde m} \frac{d \upsilon_\theta}{dt}=mg^*(\mu\sin\theta-\cos\theta),
\qquad r \frac{d \theta}{dt}=\upsilon_\theta -{\upsilon}(t) \cos\theta,
\end{equation}
where $g^*=g+\alpha(t)$ in upward acceleration phase of the container (phase~I: $t= 0\sim T/4$) and $g^*=g$ in downward acceleration phase (phase~II: $t=T/4 \sim 3T/4$). The latter reflects the non-adhesive character of dry granular material. Note that the motion in the plane passing the centers of sphere is considered here for brevity, and the balance of force in the radial direction at the contact point,
\begin{equation}\label{eq_3_16}
N-mg\sin\theta-m\alpha \sin\theta=0,
\end{equation}
is taken into account. Note also that we consider $\mu$ in the range from zero to $\sqrt{3}$ because relevant angular position of sphere~$B$ is limited between $\theta =\ang{30}$ to zero.

\subsubsection{Trajectory of the contact point}\label{trajectory}

We perform a numerical simulation of eq.~\eqref{eq_3_15} under the initial condition $\theta=\ang{30}$ and $v_\theta =0$, and examine the displacement of sphere~$B$ for a given $f$ and $a$. Figures~\ref{Intrusion}(c)--(f) are examples of the motion of contact point~$P$ along the surface of sphere~$A$. The abscissa is the time $t$, and calculation is made from $0$ to $T/2$. Behaviors of~$P$ for different $m^*$ ($\equiv {\tilde m}/m)$, $\mu$, $fa$ and~$f$ (and hence $a$) are compared. For a specified values of $fa$, contact point angle $\theta$ decreases monotonically from $\ang{30}$ to $\ang{0}$ in lower frequency cases (or in larger amplitudes $a$) through phase~I (Figure~\ref{Intrusion}(c)). Sphere~$B$ intrudes up to side-by-side positions ($\theta = \ang{0}$) of the sphere~$A$ on the bottom wall, which yields the maximum dilatation. This happens for frequencies less than $f_c$, which is determined by the combination of $m^*, \mu$ and $fa$. The amount of dilatation saturates for $f \le f_c$. On the other hand, sphere~$B$ ceases to be accelerated at an intermediate position during phase~I in higher frequency (or in smaller amplitudes), but continues moving under gravity (phase~II), so that it may stop at a certain angle $\theta_s $ (Figures~\ref{Intrusion}(d)--(f)). Figure~\ref{Intrusion}(e) compares the behavior of sphere~$B$ with different friction coefficient $\mu$ under the same $fa$ ($=10$) and $m^*$ ($=1$), which shows the decrease of intrusion as $\mu$ increases. Figure~\ref{Intrusion}(f) compares the behavior of sphere~$B$ with different effective mass $m^*$ under the same $fa$ ($=10$) and $\mu$ ($=1$). Generally speaking, development of $\theta$ asymptotes to the one as $m^*$ increases, while the difference between $m^*=1$ and $m^*=2$ seems larger. In the following, we shall mainly show $m^*=3$ cases. The quantitative difference associated with different $\mu$ and $m^*$ will be discussed in later sections.

\subsubsection{Mutual intrusion and dilatancy}

By the mutual intrusion of spheres, dilatation in the horizontal direction is induced because of the excluded volume effect. Under given $f$ and $a$, we first calculate the terminal position $\theta_s$ of the contact point (Figure~\ref{Intrusion}(c)--(f)). In our simplified model, the rate of dilatation is given by $\Delta \equiv \Delta L/L=2\cos(\theta_s)/\sqrt{3}-1$. Figure~\ref{Dilatancy}(a) shows the relation between the terminal position $\theta_s$ and $f$ in the $m^* =3$ and $\mu = 1$ case. Near the $\theta_s = \ang{30}$ region, the amount of dilatation will not likely be large enough to recognize undulations. The latter corresponds to smaller $fa$ case with larger $f$ and smaller $a$. Given the amount of intrusion, dilatancy $\Delta$ is calculated as shown in Figure~\ref{Dilatancy}(b). Here, $fa$ and $\Delta$ are chosen for abscissa and ordinate, respectively. The choice of $fa$ as the ordinate axis has an intention to examine whether it is appropriate to describes the onset of undulation. Smallest $\theta_s$ ($\approx 0$) corresponds to the largest intrusion $\Delta$ ($=0.1547$), which is attainable in a wide range of $f$ by the choice of $a$ under the constraint $fa=\mathrm{constant}$. On the other hand, $\theta_s \approx \ang{30}$ corresponds to smaller displacement $\Delta \approx 0$, which is realized in smaller $fa$ area. The latter gives a rather strong restriction on the choice of $f$ and $a$. In any case, the same $\Delta$ is realizable by suitable combinations of $f$ and $a$ provided that the above conditions are satisfied. The red open circles show the lower boundary of the dilatation for a given $fa$, which saturates to the maximum possible value $2/\sqrt{3}-1$ ($= 0.15470$) as $fa$ increases. Figure~\ref{Dilatancy}(c) is the comparison of $\Delta$ for different friction coefficient $\mu$ (simulated with $m^*=3$), and Figure~\ref{Dilatancy}(d) is the cumulative ones for various combination of $m^*$ and $\mu$. Throughout our numerical simulation, we found qualitatively the same dependence of $\Delta$ on $\mu$ and $m^*$. Namely, $\Delta$ asymptotes monotonically to the one with the increase of $\mu$ and/or $m^*$, although the behavior in $m^*=1$ and $\mu \lesssim 1$ cases differ considerably from others.

\begin{figure}[h]
\includegraphics[height=50mm]{ContactPoint}
\includegraphics[height=50mm]{DilatationM3mu1}

\includegraphics[height=50mm]{Fig-m3}
\includegraphics[height=50mm]{Fig-All-M-mu}

\includegraphics[height=50mm]{fac-mu}
\includegraphics[height=50mm]{fac-m}
\caption{(a)~Dependence of the terminal position $\theta_s$ of the contact point on the frequency $f$ under given $fa$. (b)~Dependence of the dilatancy $\Delta$ on $fa$ ($m^*=3$, $\mu=1$). (c)~Comparison of $\Delta$ for different friction coefficient $\mu$ ($m^*=3$), and (d)~comparison of $\Delta$ for various combination of $m^*$ and $\mu$. Onset of undulations: (e)~dependence of $(fa)_c$ on friction $\mu$ and (f)~that on effective mass $m^*$.}\label{Dilatancy}
\end{figure}

\subsubsection{Extrapolation to \texorpdfstring{$\Delta=0$}{Delta=0}}\label{extrapolation}

The lower boundary curve of the dilatation mentioned above is well fitted by
\begin{equation}\label{eq_3_17}
\Delta = \Delta_0 + {\tilde \alpha} \exp(-{\tilde \beta} fa),
\end{equation}
where $\Delta_0=2/\sqrt{3}-1$ ($=0.1547\dots$). Fitting constants ${\tilde \alpha}$ and ${\tilde \beta}$ are calculated for various combinations of $m^*$ and $\mu$. By the extrapolation of the curve to $\Delta=0$, the critical value $(fa)_c$ that describes the onset of undulations should be obtained. Figures~\ref{Dilatancy}(e) and~(f) show the dependence of $(fa)_c$ on $\mu$ and $m^*$, respectively, which reveals rather larger difference between $m*=1$ and $m*>1$. Except for the $m*=1$ data, critical values of the onset of undulation obtained are $(fa)_c = 2.5$, $2.6$, $2.7$, $2.9$ within the accuracy of $\pm 0.1$ for $\mu= 0.5$, $1$, $1.5$, $1.73$, respectively. We have shown that experimentally determined $(fa)_c$ ranges from~16 to~44 (Figure~\ref{AspectRatio}(b)), which are one order larger than those obtained in the above estimation. The disagreement may come from (1)~the accuracy of calculation and the range of applicability of our numerical simulation, and (2)~the realizability in experimental circumstances. As for the first point, eq.~\eqref{eq_3_17} yields a linear relation of $fa$, $\Delta=-{\tilde \alpha}{\tilde \beta}\parens[\big]{fa-(fa)_c)}+\dotsb$ near the critical point, where $(fa)_c=(\Delta_0+{\tilde \alpha})/{\tilde \alpha}{\tilde \beta}$. This expression certainly describes the linear relation of $fa$, but needs accurate calculation to assure the slope and intercept of linear fitting function in this region (see also Figure~\ref{ElasticPlate}(b) for the case of elastic plate). Thus, the trajectory of the contact point for very small $fa=\mathrm{constant}$ becomes less accurate for smaller $a$ and larger $f$ case, because the latter requires much smaller time-step, which implies that the extrapolation of $\Delta$ vs.\ $fa$ curve down to $\Delta=0$ owes much to the behavior in larger $fa$ region, and lessens the accuracy of estimation of $(fa)_c$. Second, experimentally observable boundary of the granular material is not smooth, which implies that the marginally distinguishable displacement of the latter is an averaged one over a few diameter size. The amplitude of external oscillation $a$ also needs the accuracy of the same or higher order of the grain size $d$, i.e., $ a = C_1 d$, where $C_1$ is of the order of~1. Meanwhile, the characteristic time scale of undulations is $2T$, during which the grains adjacent to the container bottom displace up to the height of $\delta = C_2 d$ (see Figure~\ref{universal}(e)), where $C_2$ will be of the order of~10. Assuming the free flight under the gravity $g$, the latter gives $f\approx \sqrt{g/(2C_2 d)}$, and hence $fa \approx (C_1/\sqrt{2C_2}) \sqrt{gd}\sim\sqrt{gd}$. For the typical experiment on the granular material of grain size $d\lesssim \qty{1}{[mm]}$, the latter gives $fa \lesssim \qty{100}{[mm/s]}$. In the customary experiment on vibrating granular material, reliable spatial resolution of the amplitude of external oscillation are $a \sim \qty{0.1}{[mm]}$. With this amplitude, the necessary frequency $f$ to achieve the condition $fa=10$ amounts to $f=\qty{100}{[Hz]}$, which is quite high in the experiment of undulations. On the other hand $f=\qty{1}{[Hz]}$ is realizable, but $a=\qty{10}{[mm]}$ (for $fa=10$) will be too large to maintain the connectedness of the granular material. Our numerical simulation also shows that the cases (large $f \times$ small $a$) and (small $f \times$ large $a$) under the same $fa$ are not equivalent because the type of acceleration $g^*$ and period applied to the test particle could be different between the two (Figure~\ref{Intrusion}(c)--(f)). These contemplations imply that $fa$ description may not be appropriate as far as the onset of undulation is concerned.

\subsubsection{Contour map approach}%\label{contour}

In this subsection, we shall examine the onset of undulation in the $f$-$a$ plane. Figure~\ref{contour}\eqref{contour_a} is the contour lines of $\Delta$ for $\mu=1.0$ case obtained in our simulation with $m^*=5$. The abscissa and ordinate are the frequency $f$ and amplitude $a$ of the external forcing, respectively, and the values of $\Delta$ are attached to each lines. We also superpose the lines $fa=10$, $20$, $30$, etc.\ in black dotted lines. Along the line $fa=\mathrm{constant}$, $\Delta$ is smaller for larger $f$, while it is larger for larger $a$. If we adopt the realizable experimental condition on the amplitude of oscillation as $a_0=\qty{0.2}{[mm]}$, we only need to pay attention to the line $a=0.2$ in Figure~\ref{contour}\eqref{contour_a}. Furthermore, if we assume the level of recognition of undulations $\Delta_0$ as 0.10, 0.11, 0.13, and 0.15, then the critical frequency $f_c$ becomes 39.8, 28.4, 20.3, and \qty{15.9}{[Hz]}, respectively (shown by triangles on the relevant positions in the figure). This implies that e.g.~small dilatation $\Delta=0.10$ can be recognized by choosing $f$ as large as \qty{39.8}{[Hz]} even if $a$ is as small as 0.2. The products of $f_c$ and $a_0$ ($=0.2$) yield $(fa)_c= 7.96$, $5.68$, $4.06$, and $3.18$, respectively, corresponding to the above-mentioned $\Delta_0$ in this order. The value $(fa)_c=7.96$ looks a few times smaller than the one obtained in the experiment (Figure~\ref{AspectRatio}), although which choice of $\mu$ in our numerical model is appropriate to compare with experiment is not clear at this moment.

Figure~\ref{contour}\eqref{contour_b} is the superposition of contour lines for $\mu = 0.1$ (green), $1.0$ (blue) and $1.73$ (red), where contour lines of $\Delta=0.09$, $0.11$, $0.13$ and $0.15$ are shown in solid lines, while others are in broken lines. Values of $\Delta$ are attached to respective lines. We also superpose the lines $fa=10$, $20$, $30$, etc.\ in black dotted lines. An example of the positions of recognizable undulation under the assumed observation level $\Delta_0=0.11$ is marked by triangles on $a=0.2$ line. In this case, critical frequencies $f_c$ are 11.2, 28.4 and 36.1 for $\mu=1.73$, $1.0$ and $0.1$ cases, respectively. The corresponding $(fa)_c$ values are 2.24, 5.68 and 7.22, respectively. The frequency range of undulations $f$ examined in our simulation covers the experimental ones adopted by many researchers (Figure~\ref{f-a-diagram}), but $(fa)_c$ obtained in the present simulation seems a few times smaller (Figure~\ref{AspectRatio}), which is left for our future study. Nevertheless, we can expect that the onset of undulations $(fa)_c$ is smaller for materials with larger $\mu$, and vice versa, and the latter ordering does not change as far as we checked down to $a=0.02$ with $\Delta_0=0.02$ in the present numerical simulation. Although we have no available data on the effective friction coefficient $\mu$ of the given material in a granular state, we can speculate $\mu_{\Glass}<\mu_{\Aluminum} \lesssim \mu_{\Lead}$ by comparing Figure~\ref{AspectRatio}(b) with Figure~\ref{contour}\eqref{contour_b}.

\begin{figure}[h]
\begin{subfigure}{.58\textwidth}
\includegraphics[width=\textwidth]{ContourA}
\caption{}\label{contour_a}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{.58\textwidth}
\includegraphics[width=\textwidth]{ContourB}
\caption{}\label{contour_b}
\end{subfigure}
\caption{Contour map of dilatation on the $f$-$a$ plane ($m^*=5$). \eqref{contour_a}~Example of the contour lines for $\mu=1.0$. Triangular marks show the critical point where undulations are recognized under given level of observation. \eqref{contour_b}~Superposition of contour lines for $\mu=0.1$ (green), 1.0 (blue) and 1.73 (red). Triangular marks on $a=0.2$ line show the positions of recognizable undulations under the observation level $\Delta=0.11$. For details, see the text.}\label{contour}
\end{figure}

In our present model, the difference of material constants are not incorporated except for the friction coefficient. The different properties of material are important and acquire a huge amount of studies in the solid state physics, which does not seem the case for their properties in granular state. For example, material properties of ($\text{lead}:\text{aluminum}:\text{glass}$) are $(1:0.24:0.22)$ for the density $\rho$, $(1:4.4:4.4)$ for the Young's modulus $Y$, and $(1:4.7:5.2)$ for the shear modulus~$G$, respectively. In contrast to rather closer material properties between aluminum and glass spheres, lead and aluminum spheres show similar granular properties in our experiment, which implies that the material constants in solid state are not important for the macroscopic properties of the granular material. One possible cause of discrepancy between glass beads and others, however, may be the influence of electricity, because the container used in our experiment is made of transparent acrylic plates for optical observation, where frictional electricity between the container and glass beads, among others, may not be negligible. The latter introduces factors other than mechanical ones, which is left for future studies.

In our numerical model, only a few layers near the bottom wall are included. In the experiment, about ten layers of grains are piled, so that the movement of the spheres in the lower layers must be suppressed by the above ones, the extent of which will depend on the stress network. Without detailed information on the structure of the latter, we include the interaction of the test sphere by taking account of the averaged effect from the other spheres. In an analogy to the virtual mass employed in fluid mechanics, we have introduced an effective mass $m^*$ in eq.~\eqref{eq_3_15}. The validity of the choice of $m^*$, however, is left for future studies.

\goodbreak

In Figure~\ref{universal}(a), we have shown rather clear undulation arch observed in a layer of sesame seeds, whose aspect ratio also shows a dependence on $fa$ with a certain $(fa)_c$ (see Figure~\ref{AspectRatio}(b)), although the material has non-spherical shape. In addition to the mutual intrusion of grains similar to the spherical ones, the former are highly associated with the change of orientation, and the tangential stress between flat sides of the grains seems to play an important role. This requires another mechanism of momentum transfer along the layer, which is left for future investigation.

\section{Conclusions}

Vertically vibrated granular material confined in thin vertical parallel walls is investigated both experimentally and theoretically. We obtain the following.
\begin{enumerate}[leftmargin=*]
\item In contrast to the vertically vibrated thin \emph{horizontal} granular layer, where typical planer patterns are satisfactorily characterized by the non-dimensional acceleration $\Gamma$ ($\equiv (2\pi f)^2 a/g$), vertically vibrated granular material confined in thin \emph{vertical} parallel walls shows ``ripples'' and ``undulations'' depending on the frequency $f$, amplitude $a$ and the timing of sinusoidal forcing on the granular layer. In the undulation, we found that the aspect ratio of arch height~$\delta$ to wavelength $\lambda$ is found to be well scaled by means of the momentum (proportional to $fa$).
\item Density wave along the layer is observed, whose propagation speed $V_{\granular}$ is much smaller than the one $V_{\solid}$ observed in the same solid material. For example, we obtained $V_{\granular} \sim \qty{40}{[m/s]}$ in the granular state of lead spheres, while $V_{\solid}\sim \qty{4000}{[m/s]}$ inside the block of the same material. According to the theory of continuum, $V_{\solid} =\sqrt{Y/\rho}$ in 1D elastic material, where $Y$ is the Young's modulus and $\rho$ is the density. If we regard the granular material as a continuum with effective Young's modulus $Y^*$ and effective density $\rho^*$, we can estimate $Y^*$ by measuring $V_{\granular}$ and $\rho^*$ (smaller than $\rho$ but of the same order of magnitude). Then, the granular material may be regarded as quite soft material.
\item Assuming an analogy to the kinetic theory of gasses, along with the buckling theory of 1D elastic body, we can derive a linear dependence of the aspect ratio $\delta/\lambda$ on the excess momentum (proportional to $fa - (fa)_c$) for well-developed undulations, where $(fa)_c$ is proportional to the critical momentum necessary for the onset of oscillation. Here, density wave plays an important role in transmitting the momentum of the oscillating container throughout the granular layer.
\item In the undulations, a few layers of particles near the bottom wall look decisively important. We propose a simple model (``intrusion model''), in which the spheres in the lowest positions are pushed up by the upward motion of the container bottom. In the very early stage, the former intrudes into the spheres in the upper second layer, which dilates the layer horizontally. The vertical momentum given by the oscillating container is transmitted along the layer as a density wave, which applies normal force to the side walls. Resulting reaction force from the side walls pushes the layer and creates an arch structure, if it exceeds a certain magnitude provided that the deformability of the grains is not too large.
\item Equations of motion of a test sphere is proposed to elucidate the onset of undulations, and influence of the dynamic friction coefficient $\mu$ on the critical value $(fa)_c$ is examined. To do this, direct interaction between the test sphere and the adjacent one in contact with the bottom wall is considered, where effect of the other spheres are incorporated macroscopically into the test sphere as an ``added mass'' $\Delta m$. Numerical simulation elucidates the influence of frequency $f$ and amplitude $a$ on the movement of test spheres under given $\mu$ and $\Delta m$.
\item Based on our numerical simulation, contour map of dilatation $\Delta$ in the $f$-$a$ plane is made, which suggests that undulation occurs in smaller $(fa)_c$ for materials with larger~$\mu$, and vice versa. The effective friction coefficients of the material in a granular state shown in Figures~\ref{f-a-diagram} and~\ref{AspectRatio} are estimated as $\mu_{\Glass}<\mu_{\Aluminum} \lesssim \mu_{\Lead}$, although possible influence of electricity on the glass spheres may not be removed.
\item We have approached the onset of undulation of granular layer from the side which is in a certain fluid-like state. However, the other side of the critical point is in solid-like state, where static friction and gravity forces prevail. At the critical point, there should be solid to fluid transition. Although we have good fitting on the linear $\delta/\lambda -fa$ relation in supercritical state of undulations (Section~\ref{kinematic}), it is not clear what happens at the critical point. In fact, an elastic plate shows nonlinear relation between $\delta/\lambda$ and dilatation $\Delta$ (Figure~\ref{ElasticPlate}). For better understanding of the critical sate, more accurate calculation and experiment will be required.
\end{enumerate}

\appendix

\section{}

In Section~\ref{linear}, we have briefly remarked an analogy between our vibrating undulation with the de Broglie wave. Here, we shall consider a few more.\ It is well-known that a beam of light shows a classical wave-like characters that exhibit reflection, refraction and interference, etc., but unable to explain photoelectric emission. To overcome this difficulty, Einstein developed Planck's ideas further and put forward the quantum theory of light (see e.g.,~\cite{Beiser}). According to him, light is composed of a small unit (wave packet, now known as photon) carrying a fixed amount of energy
\begin{equation}\label{eq_A_1}
E=hf,
\end{equation}
where $h$ is the Planck's constant, and a letter $f$ denotes the frequency as is used in this paper, although a letter $\nu$ is the conventional notation of frequency in quantum mechanics. Based on this hypothesis, Einstein successfully explained the features of photoelectric emission. In his theory, the energy of photoelectron $E$ due to incident light of frequency $f$ on the surface of the metal is given by
\begin{equation}\label{eq_A_2}
E = hf- W,
\qquad W=hf_c,
\end{equation}
where $W$ is called the work function equal to the minimum energy that is necessary for an electron to escape from the surface of metal, i.e., $f_c$ is the threshold frequency for the photoelectric effect to be recognized. After long dispute on the duality of light, a conclusion is reached such that the wave theory and the quantum theory complement each other.

If each photon has energy $hf$ which is equivalent to $mc^2$ ($m$: mass, $c$: velocity of light), it follows
\begin{equation}\label{eq_A_3}
\frac{hf}{c}=\frac{h}{\lambda}=mc,
\end{equation}
where $\lambda$ is the wavelength, and the last term of eq.~\eqref{eq_A_3} implies the momentum of a photon. De Broglie developed this idea, and proposed a theory that all particles have both wave-like and particle-like natures, such that the particle with momentum $p$ is linked to the wavelength by the relation
\begin{equation}\label{eq_A_4}
\lambda=\frac{h}{m \upsilon}=\frac{h}{p},
\end{equation}
where $m$ and $\upsilon$ are respectively the mass and velocity of the particle.

Further development of the idea of de Broglie leads Schr\"{o}dinger to postulate a basic principle. In the one-dimensional time-dependent case, it is described as
\begin{equation}\label{eq_A_5}
i \hbar \frac{\partial \Psi}{\partial t}
= - \frac{\hbar^2}{2m}\frac{\partial^2\Psi}{\partial x^2}+U \Psi,
\end{equation}
where $\Psi$ is called the wave function of the particle, $U$ is the potential energy and $\hbar=h/(2\pi)$. If we assume $\Psi=A \exp\bracks[\big]{i(kx-\omega t)}$, we have
\begin{equation}\label{eq_A_6a}
\omega= \frac{\hbar}{2m}k^2+\frac{U}{\hbar},
\end{equation}
or
\begin{equation}\label{eq_A_6b}
k^2=\frac{2m}{\hbar} (\omega-\omega_c),
\end{equation}
where $\omega_c=U / \hbar$. In terms of $\omega=2\pi f$ and $k=2\pi / \lambda$, eq.~\eqref{eq_A_6b} is rewritten as
\begin{equation}\label{eq_A_7a}
\lambda=\sqrt{\frac{\pi \hbar}{m}}\, \frac{1}{\sqrt{f - f_c}},
\end{equation}
or
\begin{equation}\label{eq_A_7b}
\frac{\delta}{\lambda}=\delta \sqrt{\frac{m}{\pi \hbar a}}\, \sqrt{fa - f_c a},
\end{equation}
where $f_c=2\pi U /\hbar$. Here $\delta$ and $a$ are introduced only for the purpose of comparison with our result eq.~\eqref{eq_2_1}. Different type of dependence on $fa$ between eq.~\eqref{eq_A_7b} and eq.~\eqref{eq_2_1} implies that our granular layer is not governed by the Schr\"{o}dinger equation.

If light beam is regarded as an aggregate of photons along the ray that carries energy and momentum, does vibrated thin granular layer (aggregate of grains) exhibit similar relation? In order to compare the microscopic view with the macroscopic one, we consider a thin elastic rod placed along the $x$ axis (on average) and examine its bending motion, where we denote the displacement in the $z$ direction by $z(x,t)$. We assume a gentle bending of the rod, whose cross-sectional area $S$ and geometric moment of inertia $I$ are constant. The tangential stress $f$ on the cross-section at the position $x$ is given by
\begin{equation}\label{eq_A_8}
f=\frac{1}{S} \frac{\partial M}{\partial x},
\qquad
M=\frac{YI}{R} \approx YI \frac{\partial^2 z}{\partial x^2},
\end{equation}
where $R$ is the radius of curvature, $Y$ is the Young's modulus, and $M$ is the moment to bend the rod (Bernoulli--Euler's law). Then the equation of motion of the small element $\Delta x$ of the rod with density $\rho$ is
\begin{equation}\label{eq_A_9}
(\rho S \Delta x)\frac{\partial^2 z}{\partial t^2}=-f(x+\Delta x,t)S+f(x,t)S
\approx - \frac{\partial f}{\partial x} S\Delta x = - \frac{\partial^2 M}{\partial x^2} \Delta x,
\end{equation}
which yields
\begin{equation}\label{eq_A_10}
\rho S \frac{\partial^2 z}{\partial t^2}
= - \frac{\partial^2}{\partial x^2} \parens[\Bigg]{YI \frac{\partial^2 z}{\partial x^2}}.
\end{equation}
If the material is uniform, eq.~\eqref{eq_A_10} becomes~\cite{LL7}
\begin{equation}\label{eq_A_11}
\frac{\partial^2 z}{\partial t^2} + \kappa^2 \upsilon^2 \frac{\partial^4 z}{\partial x^4}
=\parens[\Bigg]{\frac{\partial}{\partial t}-i \kappa \upsilon \frac{\partial^2}{\partial x^2}}
\parens[\Bigg]{\frac{\partial}{\partial t}+i \kappa \upsilon \frac{\partial^2}{\partial x^2}}\, z=0
\qquad \parens[\Big]{\upsilon \equiv \sqrt{\frac{Y}{\rho}}, \quad \kappa \equiv \sqrt{\frac{I}{S}}}.
\end{equation}
Equation~\eqref{eq_A_11} is reduced to a pair of Schr\"{o}dinger equations, and the dispersion relation assuming $z=z_0 \exp\bracks[\big]{i(kx-\omega t)}$ yields $\kappa \upsilon k^2=\omega$ or $\lambda=\sqrt{2\pi \kappa \upsilon / f}$, the latter being the same form as eqs.~\eqref{eq_A_6a}--\eqref{eq_A_6b} or~\eqref{eq_A_7a}--\eqref{eq_A_7b} with $U=0$.

If we further assume a wave-like dependence $\rho=\rho_0 \exp\bracks[\big]{i(\xi x-\sigma t)}$ and $Y=Y_0 \exp\bracks[\big]{i(\xi x-\sigma t)}$ that reflects the density wave along the rod (while $S$ and $I$ are kept constant as before), eq.~\eqref{eq_A_10} becomes
\begin{equation}\label{eq_A_12}
\frac{\partial^2 z}{\partial t^2} + \kappa^2 \upsilon^2 \parens[\Bigg]{\frac{\partial^4 z}{\partial x^4}
+2 i \xi \frac{\partial^3 z}{\partial x^3}-\xi^2 \frac{\partial^2 z}{\partial x^2}}=0,
\end{equation}
which yields a dispersion relation $\omega^2=\kappa^2 \upsilon^2 (k^4+2\xi k^3 +\xi^2 k^2)$. For $k \gg 1$ the dispersion relation is the same as the one mentioned above. However, for $k \ll 1$, we have
\begin{equation}\label{eq_A_13a}
\omega \sim \kappa \upsilon \xi k,
\end{equation}
or
\begin{equation}\label{eq_A_13b}
\lambda \sim \frac{\kappa \upsilon \xi}{f},
\end{equation}
which is of the same form as the de Broglie wave eq.~\eqref{eq_A_4}. Except for the proportional constant, eq.~\eqref{eq_A_13b} yields
\begin{equation}\label{eq_A_14}
\frac{\delta}{\lambda} \sim \parens[\Bigg]{\frac{\delta}{\kappa \upsilon \xi a}}\, fa,
\end{equation}
which is of the same form as eq.~\eqref{eq_2_1} with $(fa)_c=0$, as should be because elastic rod conserves energy. The appearance of the threshold term $(fa)_c$ may be reminiscent of the work function in the photoelectric effect. However, the most profound difference between the micro and macro worlds will be the presence of energy dissipation. In our vibrated granular layer, inclusion of the maximum static friction will presumably be the key point left for future investigation to elucidate the mechanism on the onset of undulation.

\section*{Declaration of interests}

The author does not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and has declared no affiliations other than their research organizations.

\printbibliography

\end{document}