Undulating compression and multi-stage relaxation in a granular column consisting of dust particles or glass beads

For fundamentally characterizing the effect of hierarchical structure in granular matter, a set of compression-relaxation tests for dust particles and glass beads confined in a cylindrical cell was performed. Typical diameter of both grains is approximately 1~mm. However, dust particles are produced by binding tiny ($\sim 5$~{\textmu}m) glass beads. The granular columns were compressed with a piston until reaching a maximum load force of 20~N with a constant compression rate $v$ ($0.17 \leq v \leq 2000$~{\textmu}m~s$^{-1}$). After that, the piston was stopped and the relaxation process was quantified. From the experimental results, we found that the compression force $F$ nonlinearly increases with the increase of compression stroke $z$ depending on particles; $F\propto z^{\alpha}$ for dust particles and $F \propto \exp(z/z_G)$ for glass beads ($\alpha=2.4$ and $z_G=70$~{\textmu}m). Besides, periodic undulation and sudden force drops were observed on $F(z)$ in dust particles and glass beads, respectively. Growing manners of periodic undulation and force drops were identical to those of mean compression forces, i.e., power law in dust particles and exponential in glass beads. In addition, the undulation amplitude and wavelength decreased as $v$ increased in dust-particles compression. On the basis of experimental results, we also discuss the physical meaning of granular-compression models used in engineering fields. The relaxation process was characterized by an exponential decay of the stress followed by a logarithmic dependence one in both kinds of particles. Physical meanings of these experimental findings were discussed based on the difference between dust particles and glass beads.


I. INTRODUCTION
Ubiquity and importance of granular matter in various natural phenomena have been recognized by recent numerous studies on granular behaviors [1]. Granular matter has been usually defined by a collection of rigid dissipative particles. The principal sources of dissipation are inelastic collisions and friction. While this type of typical granular matter exhibits various counterintuitive behaviors, a different type of granular matter, hierarchical granular matter, has also been studied recently.
Hierarchical granular matter consists of granules of tiny particles. When we treat tiny particles, various cohesive forces (e.g., capillary force, electrostatic force, van der Waals force) work as a binder to form granules. By collecting the formed granules, hierarchical granular matter is comprised. Hierarchical granular matter is characterized by two distinct length scales. One is the length scale of tiny monomers that construct porous dust particles. The other length scale is the size of dust particles. Due to this complexity, behaviors of hierarchical granular matter could become quite different from conventional granular matter. In general, understanding of the hierarchical structure is a crucial key to reveal the complex nature of various phenomena. Therefore, we focus on the role of the hierarchical structure in the physics of granular matter.
Large dust particles, usually called dust aggregates, have attracted planetary physicists because they are plausible candidates for the material of planetesimal formation [2,3].
To reveal the physical properties, mechanical characteristics on dust aggregates have been studied [4][5][6][7][8]. Recently, the existence of such porous particles/aggregates on small celestial bodies (comets/asteroids) has been reported by planetary explorations [9,10]. If comets consist of hierarchical dust particles, mechanical and thermal properties must be influenced by the hierarchical structure [11,12]. Furthermore, a granular porous projectile impacting onto a granular target produces particular crater shapes as experimentally demonstrated in [13]. Undoubtedly, the fundamental understanding of dust particles/aggregates is crucial to elucidate the planetary-scale phenomena as well as laboratory-scale phenomena.
The impact response of the collection of dust aggregates (hierarchical granular matter) was also studied recently [14]. In the study, impact-induced expansion of granular cluster consisting of dust particles or glass beads was measured and the similarity between dust particles and glass beads was reported. However, whether this similarity is truly universal or not should be checked by other experimental conditions. By revealing fundamental features of hierarchical granular matter through the systematic investigation, we might be able to establish a novel type of soft matter. Actually, one can easily find various entities of hierarchical granular matter even in our everyday life, e.g., snow balls, granules of sugar, medicines, etc. Thus, the fundamental study of hierar-chical granular matter affects various kinds of science and industries.
As a basic methodology, uniaxial compression test of granular matter has been frequently performed both in soil mechanics [15] and powder engineering [16]. Because these fields have been independently developed, different models have been applied to the analysis of similar compression data. In soil mechanics field, compression index has been used to characterize the pressure-dependent void shrinkage in granular (soil) material [15]. In powder engineering field, on the other hand, the compressibility of void space in granular matter has been characterized by using bulk modulus [16]. While these quantities have been independently used in each research field, their similarity and difference have not been sufficiently analyzed so far. Fundamental physical consideration is necessary to properly compare and unify these models. Of course, uniaxial compression test is one of the most fundamental methods to characterize soft materials. To characterize mechanical properties of novel soft materials such as hierarchical granular matter, uniaxial compression tests should be performed at first.
In the case of brittle materials, propagation of upward compaction bands was reported during vertical (uniaxial) compression of snow [17] and cereals [18]. Numerical simulations are able to reproduce this propagating band and have been helpful to determine the conditions required to observe such dynamics [19]. This propagating band is characterized by sudden drops in the compression stress produced by the breakage of particles at the bottom of the bed that allows its relaxation and further rearrangement.
In this study, we carried out a set of simple compression-relaxation tests of granular columns that consist of dust particles or rigid glass beads. We found nonlinearly growing compression force both for dustparticles and glass-beads. Physics of nonlinear growth in compression force is examined by comparing the models used in engineering fields. Furthermore, details of compression-force behaviors depend on constituent particles. Dust-particles compression shows periodic undulation whose amplitude grows as compression proceeds. However, glass-beads compression exhibits quasi-periodic force drops due to the stick-slip rearrangement of grains network structure. To understand the physical nature of the observed dust-particles force undulation, fracturing of dust particles during the compression was also analyzed based on the image analysis. Concerning the relaxation process, the final equilibrium state is reached relatively faster in the dust particles than the glass beads because the former can easily be deformed and many of them are pulverized during the compression process. In this paper, we report on the analysis results of these compression-relaxation behaviors. Then, the physical meaning of the observed results is considered.

II. EXPERIMENT
Simple compression-relaxation tests of a column consisting of hierarchical granular matter or glass beads were performed. To form hierarchical granular matter, dust particles were generated using a pan-type granulator (AS-ONE, PZ-01R). Tiny glass beads of diameter 5 µm (Potters Ballotini, EMB-10, size range is 2-10 µm) were poured on a pan and the pan was rotated. During the rotation, very small amount of water was sprayed on the rotated beads. Then, dust particles were formed by granulation of monomers (tiny glass beads). Monomers are presumably bound by humidity-induced micro capillary bridges and van der Waals force among monomer particles. The dust particles were sieved to collect particles with diameter d in the range of d = 0.6-1.4 mm. In this paper, we denote the typical particle diameter d = 1 mm. Dust particles made by this procedure are stable once they are formed under the usual atmospheric conditions. The resulting particles from the above protocol approximately have grain-scale packing fraction (in the particles) φ g 0.27 (similar to [14]) but their stiffness is sufficient to allow manipulation during the experiments.

A. Compression process
A cylindrical cell made of resin with inner diameter 20 mm and height 57 mm was mounted on a table of the universal testing machine (Shimadzu AG-X), see Fig.  1(a). The cell was filled with dust particles. Specifically, fixed-mass dust particles (8 g) are gently poured in the cylindrical cell. Because macroscopic dust particles are large enough, these particles are densely packed just like usual macroscopic granular matter as shown in Fig. 1(c). Then, the dust-particles column was verti-cally compressed by an aluminium piston with a constant rate. The compression rate v was precisely controlled by the testing machine and varied from 0.17 to 2000 µm s −1 (over 4 orders of magnitude). Because a 50 N-capacity load-cell sensor was used in the experiment, the column was compressed until the compression force became F max = 20 N. The vertical compression stroke and force were acquired by 10 or 100 Samples s −1 sampling rate. Glass beads of diameter 0.7-0.9 mm (AS-ONE, BZ08) were also used for comparison. Mass of glass beads poured in the cell was fixed at 15 g to form a nearly fixed height column. Basically, the identical experimental apparatus and protocol were used for the glass-beads experiments as well.
Besides, we performed compression tests of a single dust particle. Pictures exemplifying the compression of an individual dust particle and a granular column of these grains (hierarchical granular matter) contained in the cylinder are shown in Figs. 1(b) and (c), respectively. Note in (b) how a vertical fracture visible to the eye appears due to the compression. This measurement allowed us to determine the failure stress (tensile strength) and the corresponding yield strain of dust particles, which are approximately 15 kPa and 0.05, respectively. Therefore, the apparent elastic modulus of dust particles is about 300 kPa. These values are quite smaller than typical glass strength ∼ 7 MPa and its elastic modulus ∼ 60 GPa. See Appendix A for details of single dust particle compression test.
As can be seen in Fig. 1(b), dust particles have approximately spherical shape with rough surfaces. Thus, the macroscopic packing state should be similar between dust particles and glass beads. However, the initial macroscopic packing fraction φ m is roughly 0.6 for glass beads and 0.8 for dust particles; namely, initial bulk packing fraction of dust particles φ bulk = φ g × φ m is about 0.22. The relatively large φ m value for dust particles comes from polydispersity and deformability of dust particles.
The diameter of the cylindrical cell used in this experiment (20 mm) might not be large enough, particularly for d = 1 mm glass beads. In general, the wall friction could affect the compression force behavior. To observe the effects of wall friction and hierarchical structure, a compression test of a granular column simply consisting of tiny glass beads (diameter 5 µm) was also performed.
As can be observed in Fig. 1(c), the dust-particles column was significantly compressed and all the dust particles seem to be compressed more or less homogeneously. Thus, collective breakage of dust particles plays an essential role in compression force behavior. The collective particle breakage at the scale of the granular column allowed a visible compaction that can be recorded with a video camera [see snapshots in Fig. 1(c)], which was not possible for the glass beads. Digital Image Correlation analysis (DIC) [20] was used to determine the displacement of grains between consecutive frames (velocity field) during the compression process.

B. Relaxation process
Once the maximum compression force F max was reached, the piston was stopped and kept in that position to study the force relaxation process. The relaxation of compression force was measured as a function of time. The experiments were performed using the dust particles and glass beads. The room temperature, T = 25 • C, and humidity, ∼ 50%, were controlled and kept constant during each run.  Fig. 2(b)] subjected to different compression rates. In the case of glass beads, only small values of compression rates were used because the stiffness of the particles produces a faster increase of the force load, in comparison with dust particles in which case we can use larger compression rates. In these plots, negative values of t correspond to the compression phase and the positive ones to the relaxation phase. First, the compression process was performed by increasing the stroke z from zero to a maximum value z max (t = 0 s) at which the maximum force F max = 20 N was reached. Then, the relaxation phase began and it was characterized by the continuous decrease of F (t). For all the studied cases, the compression phase can be characterized by fluctuations of F , although the form and size of these fluctuations depend clearly on the nature of the granular particles. On the other hand, the relaxation phase is characterized by an initial exponential decay of F (t) that depends on the compression rate, followed by a slow and smooth decrease of F (t). Note that the relaxation to the equilibrium state was achieved faster by the dust particles than the glass beads. Because we would like to focus on the difference between dust particles and glass beads in force fluctuation and relaxation, we fixed most of the parameters. The principal parameter varied in this experiment was the compression rate v. To check the reproducibility of the experiment, slow compression tests (17 µm s −1 for dust particles and 0.17 µm s −1 for glass beads) were repeated as shown in Fig. 2 Let us first focus on the compression of granular columns and then compare the relaxation process for both kinds of materials.

A. Dust particles compression
Compression data for dust particles are shown in Fig. 3. When we plot F (t) as shown in Fig. 2, it is difficult to see the details on data behavior in compression phase (t < 0) because of the wide variation of compression rate v. Nevertheless, v-independent nature of F can FIG. 2. Compression force F as a function of time t measured at different compression rates followed by system relaxation for a) dust particles and b) glass beads. In these plots, t = 0 corresponds to the moment at which the maximum compression load Fmax = 20 N is attained by increasing the stroke z from 0 to zmax. After t = 0, the stroke is fixed and the force is measured as a function of time t.
be realized when we plot F (z). Therefore, the compression force F is plotted as a function of compression stroke z in Fig. 3. The level z = 0 is manually determined by the contact of the piston and the column's top surface, and the downward (compressive) direction is set as positive direction of z. In Fig. 3(a), force curves with various compression rates are displayed. As can be seen, F grows nonlinearly with increasing z. Although F (t) curves for various v values look quite different in Fig. 2(a), the form of F (z) is almost independent of v. This indicates that the compression proceeds in a quasi-static manner.
Another prominent feature observed in the force curves is periodic undulation of F (z). The inset of Fig. 3(a) shows the magnified view of the force curves. As seen in Fig. 3(a), amplitude of the undulation depends on v. The slower compression yields the larger amplitude. Wavelength of the undulation is in the order of mm and slightly depends on v.
To analyze these characteristic features, we fit the mean growth of F (z) to power-law function and extracted the undulation component. A typical example of this procedure is shown in Fig. 3(b), in which the force curve with v = 83 µm s −1 is shown. The red dashed curve indicates the power-law fitting of the mean growth of the compression force, F ∝ z α with α = 2.7. The undulation component computed by F − F is shown in the inset of Fig. 3(b). The nonlinear growth of undulation amplitude can also be clearly confirmed. The black dotted curve in the inset of Fig. 3(b) indicates the power-law fitting of the envelope, |F − F | env ∝ z αp with α p = 2.7. This type of power-law growth found in F and F − F can be observed in all F (z) curves of dustparticles compression. Thus, we fitted all the data with F = F * (z/d) α and |F − F | env = F * p (z/d) αp , where F * , F * p , α, and α p are fitting parameters.
FIG. 3. Force F vs. stroke z for a dust-particles-column compression. Color in (a) indicates the compression rate v as denoted in the legend. The inset of (a) displays the magnified view of periodic undulations. A typical example of F (z) (v = 83 µm s −1 ) is presented with its power-law fit (red dashed curve) in (b). The inset of (b) shows the temporal development of the periodic undulation with a power-law envelope fit (black dotted curve).
The obtained fitting-parameter values as functions of v are displayed in Fig. 4(a-c). The power-law exponents α and α p do not show any clear trend; they are rather almost constant ( Fig. 4(a)). Thus, we consider a unified value for the exponents α α p = 2.4 ± 0.7. For F * , clear v dependence cannot be confirmed again ( Fig. 4(b)). However, the data scattering particularly in small v regime is significant. Thus, only the typical order of F * can be estimated as F * 10 −2 N. The coefficient of force undulation F * p shown in Fig. 4(c) exhibits a slightly negative correlation with v. The dotted line shown in Fig. 4(c) is the power-law fitting, F * p ∝ v −β with β = 0.4 ± 0.1. This negative power-law correlation corresponds to the decrease tendency of undulation amplitude in the fast compression.
Next, we computed the dominant wavelength of the periodic undulation. By identifying the peak wavenumber in the power spectrum of F − F (inset of Fig. 4(d)), the dominant wavelength λ can be easily computed. The obtained λ values are plotted in Fig. 4(d). The dot-ted line in Fig. 4(d) indicates the logarithmic decay of λ, λ = λ * ln(v up /v), where λ * = 0.09 ± 0.01 mm and v up = 24±6 mm s −1 were computed by the fitting. Since λ becomes 0 at v = v up , v up should correspond to the upper limit of the compression rate for the periodic undulation observation. However, we did not perform such a fast compression test due to the technical limitation. By combining the above-mentioned analysis results, we finally arrive at the form of compression force, where the value of v 1 ( 10 −2 mm s −1 ) can be estimated from the fitting results. This form includes all the characteristic features found in F (z) curves of dust-particlescolumn compression.

B. Glass beads compression
To ensure the intrinsic nature of dust-particles compression, we performed similar compression tests with glass beads which are rigid enough in the current loading condition, F ≤ 20 N. The results of glass-beadscolumn compression are shown in Fig. 5. F (z) curves for glass beads with various compression rates are plotted in Fig. 5(a). The nonlinearity of F (z) curves appears to be enhanced in the glass-beads compression. Figure 5(b) shows log-linear plot of the same data. One can confirm the clear exponential growth in the late stage of compression. Moreover, the range of z to reach F max = 20 N is ∼ 60 times smaller than that of dust-particles compression. This rapid increase of F results in its exponential growth rather than the power-law increase.
A typical example of glass-beads compression with v = 17 µm s −1 is shown in Fig. 5(c). The red dashed curve represents the exponential fit. Although the notable fluctuation of F can be observed in Fig. 5(c), its behavior is quite different from the periodic undulation seen in Fig. 3. In glass-beads compression, sudden force drops are repeated, whereas the periodic fluctuation can clearly be observed in dust-particles compression. The fluctuation component (force drops) in F (z) was extracted also for glass-beads compression, as shown in the inset of Fig. 5(c). The growth of fluctuation amplitude is confirmed also in the glass-beads compression. The black dotted curve indicates the exponential growth of the fluctuation amplitude. By computing power spectrum of the fluctuation component (inset of Fig. 6(c)), we can identify the dominant wavelength of the fluctuation λ G also in the glass-beads compression.
The measured fitting parameters for glass-beads compression are shown in Fig. 6. Here, we employ the exponential growth both for the mean force and its fluctuation as, F = F * G exp(z/z G ) and |F − F | env ∝ exp(z/z Gp ), where F * G , z G , and z Gp are fitting parameters. In Fig. 6(a), characteristic length scales, z G and z Gp , are plotted as functions of v. Similar to the dust-particles compression, the growing manners of bulk compression force and its fluctuation amplitude coincide; z G z Gp = 0.07 ± 0.01 mm. The compression rate dependence of F * G is presented in Fig. 6(b). The measured F * G values show significant data scattering. As for the characteristic wavelength λ G , the measured values do not show clear v dependence (Fig. 6(c)). The average wavelength is λ G = 0.04 ± 0.02 mm. Note that the value of λ G is about one order of magnitude smaller than λ (dust-particles case). From the above-mentioned analyses, the compression force for glass-beads column can be written as, The parameter F * G is actually very sensitive to the initial condition. The value of F * G shown in Fig. 6(b) varies over six orders of magnitude (see also the fitting lines shown in Fig. 5(b)). The large variation of F * G comes from the subtlety of initial response in the compressed glass-beads layer. By definition, compression begins when the piston reaches the particle on top. The subtle structure of grains configuration affects the behavior of early-stage compression of the glass-beads column since the glass beads cannot be deformed. Wide variation of the initial grains-network structure results in the large dispersion of the compression force in the early stage. To develop a stiff internal contact network, initial particlelevel surface roughness must be flattened by compression. Figure 5(b) indicates that this preparation-related variation results in about one-particle-diameter scale fluctuation. Contrastively, the dust-particles column shows relatively good reproducibility because the particles can be easily deformed and yielded by the compression. Due to this deformability of particles, local deformation without network rearrangements is allowed. Similarly, the value of C G in Eq. (2) does not exhibit the universal behavior.
Probably, these parameters strongly depend on the initial conditions (preparation dependence) of the column.

C. Image analysis of dust particles compression
Next, we analyze the collective particle breakage that could relate to the periodic undulation observed in the dust-particles compression. Figure 7 shows snapshots obtained from DIC analysis of dust-particles compression at (a) v = 17 µm s −1 and (b) v = 2000 µm s −1 . Only 18 frames are shown in each case, but they exemplify the behavior observed during all the compression processes. In these images, the upper black part corresponds to the piston position and each vector in color indicates the displacement of pixel patches between two consecutive frames. The color scale bar at the left of the figure indicates the corresponding displacement in mm. Note in Fig. 7(a) that at slow compression, there are frames practically colored in red, which indicates that all the granular column is rearranged as a whole in that precise moment. Then, the column displacement becomes negligible in the following frame and starts to increase again until a new general rearrangement occurs. We consider this general rearrangement event corresponds to the collective breakage of compressed dust particles. Note that this type of rearrangement is measurable only in dustparticles compression. If this collective breakage is the primary origin of the observed periodic undulation, wavelength shown in Fig. 4(d) should agree with the stroke interval between consecutive general rearrangements (undulation wavelength λ). The number of frames between two consecutive general rearrangements obtained from analyzing the whole video is in average 3.9. Since the pictures were taken at 0.073 frames s −1 , the rearrangement occurs approximately every 53 s, or equivalently, when the piston stroke increases about 0.9 mm. This value is close to the wavelength reported in Fig. 4(d) for the corresponding compression rate. Therefore, the periodic undulation of F could relate to the general rearrangement of the granular column during the compression process. On the other hand, when the granular column was subjected to large compression rates, see Fig. 7(b), only a local rearrangement at the upper part of the granular column (close to the piston) was observed. This can be associated with the small amplitude observed for this compression speed in the experiments. Therefore, at slow compression, the applied stress transmitted to the grains force them to enter into existing vacancies in the whole column and there is enough time for general rearrangement, but at fast compression, only the grains close to the piston have time to enter into the local vacancies and the rearrangement occurs only in that zone. The particle rearrangement can also be induced by particles breakage, in all the column for slow compression or localized below the piston at fast compression rates. This difference is reflected in the rate-dependent amplitude decrease found in dust-particles compression (Fig. 4(c)).

D. Relaxation of dust particles vs glass beads
In the literature, different models have been proposed to describe the stress relaxation of a material previously subjected to a compression force. The generalized Maxwell model expresses the relaxation of stress σ as a function of time as: σ(t) = σ e + Σ n j=1 C j exp(−t/τ j ), where σ e is the equilibrium stress, C j are constant coefficients and τ j are relaxation times. For practical purposes, a three-term exponential equation has been found enough to describe the relaxation process in soft materials [21]. For a granular system subjected to oscillatory compaction that allows to achieve highly packed granular systems, the stress relaxation is found dependent on the strain rate. At slow strain rates, the relaxation is logarithmic in time, while at fast strain rates, the system follows a two-step relaxation [22] given by the expression: σ(t)/σ max = A + B exp(−t/τ ) + C ln(t), where τ , A, B, and C are the relaxation time and fitting parameters; an instantaneous relaxation that decays exponentially followed by a logarithmic relaxation.
Let us analyze the relaxation process found in our experiments. Figures 2(a) and (b) show F vs. t for the relaxation of dust particles and glass beads, respectively. Whereas relaxations in very slow dust-particlescompression data (v = 1.7 and 0.17 µm s −1 ) could not be measured due the technical limitation, all the other available relaxation data were analyzed. We found that our data are not well described by the two step relaxation mentioned above. Instead, it is necessary to consider an expression consisting of two exponential terms of considerably different time scale plus the logarithmic relaxation, in the form: where, τ 1 and τ 2 are characteristic relaxation times. The values of τ 1 , τ 2 and the coefficients A, B, and C obtained by fitting the experimental data shown in Fig. 8 depend on the compression rate, and are reported in Ta (single exponential with constant and logarithmic terms) for considerably large values of τ 1 . This happens for dust particles, in which τ 1 is three orders of magnitude larger than τ 2 . From Table I, we can also notice that the main difference in the relaxation dynamics for dust particles and glass beads resides in the characteristic timescales, in particular in the fast relaxation characterized by τ 2 for dust-particles compression. This means that the relaxation in a dust-particles column is faster than in a glass-beads column. Moreover, an instantaneous relaxation with τ 2 ∼ 1 s is only found for glass beads when v = 50 µm s −1 . This fast relaxation observed at the largest compression rate would probably correspond to the two step relaxation for fast strain rates reported in Ref. [22]. Furthermore, longer time measurement is also required to precisely measure and discuss the longer relaxation time τ 1 . However, at the limit of t → ∞, F/F max diverges due to the ln(t) term in Eq. (3). Thus, the relaxation model of Eq. (3) is just for the practical timescale phenomena. A deeper analysis focused on the dependence of the parameters τ 1 , τ 2 , A, B and C on the compression rate, packing fraction, dust-particles stiffness, particle size and the maximum load are left for further research.

IV. DISCUSSION
In the compression phase, both dust-particles and glass-beads columns exhibit the nonlinear resistance force against the constant compression. However, the functional forms are different. Dust-particles and glass-beads compressions obey power-law and exponential growths, respectively. This qualitative difference originates from deformability of the constituent particles. Dust particles can be easily deformed and broken while the glass beads particle cannot be significantly deformed. Thus, only the rearrangements of particle-network structure are permitted during the compression of glass-beads column. As a result, the rapid increase of bulk compression force with sudden force drops is observed in glass-beads compression. This implies that the fluctuation observed in glassbeads compression results from the stick-slip motion of contact network. Stick-slip behavior is a typical outcome of friction-induced slipping. Although we cannot distinguish wall-particle and particle-particle slips from the data, macroscopic sudden rearrangement of glassbeads network due to the frictional slip must be an origin of this behavior. In dust-particles compression, on the other hand, the local-compression-band formation and its propagation could occur as shown in Fig. 7.
Various models for the compression test of granular materials have been proposed. In soil mechanics, the compression index C c has been measured from the relation between void ratio e = (1/φ bulk ) − 1 and applied pressure p = F/S [15], where φ bulk and S are bulk packing fraction and area of compressing piston, respectively. We found that the exponential growth of F (z) is consistent with this model, as shown in Appendix B. In the field of powder engineering, however, Heckel model has  Table I. been used to characterize the granular (void) compression [16]. In the Heckel model, porosity ε = 1 − φ bulk is related to its initial value ε 0 and applied pressure p as − ln(ε) = κ h p − ln(ε 0 ). Here, κ h corresponds to a compressibility of voids. Although the Heckel model is partly consistent with our data, we consider the compression index model is better to explain the glass-beads compression behavior. Details on the physical meaning of these models and the comparison with our experimental data are described in Appendix B.
However, these models cannot explain power-law growth observed in dust-particles compression. In this study, we empirically employed the power-law and exponential forms to simply show the coincidence of the growing manners between mean growth F and fluctuation growth F − F . This correspondence is rather natural because both mean growth of bulk compression force and fluctuation amplitude depends on φ bulk through the compression stress. In other words, the scale of plastic deformation and/or slips should relate to the scale of applied force. Furthermore, in some previous studies relating to dust compression, power-law relation between bulk packing fraction and compressive pressure has been assumed [4,6,7]. Although the situations (packing fraction regime and hierarchical structure) in these studies are different from ours, these results suggest there might be a universal power-law feature in easily-deformable dust compression.
As for the periodic undulation, similar phenomenon (periodic fluctuation of granular drag force) was also found in plowing experiment [23]. In their experiment, the amplitude of the periodic fluctuation relates to the initial packing fraction of the granular layer. And the deformation of the free surface plays a crucial role to cause the periodic fluctuation. In this study, we confined particles in the cylindrical cell. Thus, any free-surface deformation is not allowed. However, there is still sufficient free space left for local deformation and compaction in dust-particles columns. As mentioned earlier, the initial bulk packing fraction of a dust-particles column is φ bulk 0.22. Since φ bulk increases as the column is compressed, growth of undulation amplitude can be related to φ bulk . To check the quantitative relation between the plowing and compression, much more careful observation in the compacted cell is needed. This problem is open to future.
In order to understand the physical behavior of the compressed dust particles, the observed periodic undulation should be related to the collective breakage of dust particles. As demonstrated in Sec. III C, period (wavelength) of F (z) undulation observed in dust-particles compression agrees with the interval of collective particle breakage. As shown in Figs. 7 and 1(c), the particle breakage propagates to the entire column. Similar compression wave propagation has been found in compressed crushable soft materials [17][18][19]24]. Furthermore, more or less similar diffusive propagation of a compression band has been found in various conditions [25][26][27][28][29]. Therefore, we consider that the particle breakage and its propagation are the main sources of the periodic undulation found in compressed dust-particles. Although its compression-rate dependence was analyzed and formulated (Eq. (1)), fully convincing physical rationale of the obtained form has not yet been obtained. In this sense, our understanding of periodic undulation is still insufficient. According to [19], competition of the timescales is the principal source of oscillation found in the compression of crushable materials. However, this idea cannot be applied to the current case. In our experiment, the compression rate v significantly varies over four orders of magnitude and the relaxation time τ 2 is almost independent of v (Table I). Nevertheless, we observed the undulation with almost a constant λ. Timescale-independent modeling is necessary for the reasonable explanation. We consider the physical mechanism causing periodic undulation in dust-particle compression is partly similar to but somewhat different from the phenomena reported in [17][18][19]24].
Here we discuss a phenomenological description of the force oscillation. First, we consider Hooke's law, dF/dz = −k(z), where k(z) represents z-dependent (effective) stiffness of the compressed column. Next, we assume that the work dissipated by compression of z can be written as W = (δd) 2 k(z), where δ is a dimensionless constant. Then, a simple form d 2 F/dz 2 + (1/δd) 2 F = 0, which yields to the spatial oscillation (undulation) of F with wavelength λ = 2πδd, can be obtained from a relation dW/dz = F . This implies that if the energy and stiffness are related only by the characteristic length scale d (macroscopic particle size), d-dependent undulation can be induced. Put differently, k(z) /d corresponds to the effective strength (energy required to break dust particle per unit volume) of dust particles at depth z. However, this qualitative explanation is still too rough to elucidate the observed mechanical behaviors. We need more much more detail modeling to explain the growth of undulation amplitude and parameter dependence of the wavelength λ. Systematic experiments under various conditions have to be performed to reveal the physical nature of this phenomenon.
The remarkable feature found in the periodic undulation of dust-particles compression is the logarithmic dependence of wavelength λ on the compression rate v. Perhaps, this logarithmic rate dependence may originate from the rate-state-dependent friction law [30,31]. Indeed, temporally logarithmic relaxation is observed as modeled in Eq. (3). Moreover, the logarithmic relaxation is crucial particularly for dust particles, as discussed later. However, the observed undulation cannot be controlled by the timescales introduced by compression rate and force relaxation. Almost constant (at least in the same order) wavelength suggests spatial or geometrical conditions govern the undulation. To properly understand the undulation mechanism, the internal structural development of the compressed dust-particle column should be directly observed and analyzed. Unfortunately, such observation is not feasible at this time. Details on these logarithmic properties are crucial question opened to future.
Actually, we also performed single dust-particle compression tests (see Appendix A) and a compression test of tiny monomers layer without granulation (see Appendix C). However, both tests did not show any periodic undulation. Recently, compression tests of aggregates formed by mixture of snow, water, and tephra particles were also performed to discuss the mechanical properties of complex aggregates [32]. In the experiment, any periodic undulation was not found. Only the complex aggregation structure is not sufficient to produce periodicity. It can only be said that the hierarchical granular structure is necessary to induce periodic undulation. To discuss the origin of periodic undulation and its precise rate dependence, further compression tests and their statistical analysis are necessary.
Regarding the relaxation process, the slower relaxation at long timescales (larger values of τ 1 ) in the case of dust particles can be explained considering that these conglomerates are deformed and many of them pulverized during the compression process. Then, the material is closely packed and gets stuck due to the cohesion effect and the further rearrangement of the contact network takes longer time to occur. This is similar to what happens with the jammed state obtained when the system of glass beads is vibrated in Ref. [22]. For that reason, the model given in Ref. [22] is also able to predict the relaxation of dust particles. On the other hand, the glass beads remain practically intact after compression, which is reflected in a larger decrease of F and shorter values of τ 1 . On the other hand, the short values of τ 2 < 10 2 s for dust particles imply that the logarithmic relaxation starts to be important first for this kind of material. This is contrastive with the glass beads where the logarithmic relaxation is only observed for the largest compression rate, which indicates that the glass beads continue rearranging even at very longer time. The role of friction with the container and between particles could also be playing an important role in this particle rearrangement.
In this experiment, we fixed most of the parameters: particle size (d 1 mm), cylindrical cell diameter (20 mm), and constituent particle material (glass). These factors could affect the behavior of compressionforce fluctuation and its relaxation. In this study, by fixing these parameters, we have concentrated on the effect of hierarchical structure on the fluctuation of compression force and its relaxation. To completely understand the physical nature of granular compression, other parameters should also be varied and the obtained forms (Eqs. (1),(2),(3)) should be improved.

V. CONCLUSION
We experimentally compared the compressive reaction and force relaxation of dust-particles column and glassbeads column. Both granular columns exhibit the nonlinear increase of compression force. The increasing behaviors can be fitted by power law and exponential functions for dust particles and glass beads, respectively. In dust-particles compression, periodic undulation can be observed in addition to the nonlinear growth. In glassbeads compression, however, sudden force drops rather than undulation were observed during the compression. Characteristic quantities characterizing the compressionforce behaviors were systematically analyzed using the experimental data. As a result, we found that most of the characteristic quantities are roughly independent of compression rate while the compression rate was varied over four orders of magnitude. The amplitude and wavelength of the periodic undulation found in dust-particles compression depend on the compression rate. By combining all the analyzed results, compression-force behaviors were modeled by Eqs. (1) and (2) for dust-particles and glass-beads compressions, respectively. In addition, force relaxation in the compressed granular layers were measured and modeled by Eq. (3). Physical meanings of the obtained result were discussed based on various previous studies in the fields of physics, planetary science, soil mechanics, and powder engineering. Particularly, we found that the exponential force growth observed in glass-beads compression is consistent with the compression index model used in soil mechanics. Although the physical understanding of these behaviors was not completely established, reasonable understanding was obtained by considering the physical difference between dust particles and glass beads.

ACKNOWLEDGMENTS
We thank JSPS KAKENHI Grant No. 18H03679 and Nagoya University for financial support. Appendix A: Single particle compression test and stress comparison

Single particle compression test
To quantitatively assess the mechanical strength of dust particles, we performed single-particle compression tests. Because dust particles are very weak, it is difficult to pickup and handle a single small particle. Due to this technical limitation in handling dust particles, we compressed relatively large dust particles produced by the protocol mentioned in Sec. II. Dust particles of diameter d p = 1.8-7.8 mm (measured by a vernier caliper) were compressed directly by the piston (without sidewall, see Fig. 1(b)). The compression rate was set v = 1.7 µm s −1 except one trial (v = 17 µm s −1 only for 6.0 mm diameter particle). Compression force and stroke must be appropriately normalized to obtain the scale-insensitive mechanical property. Therefore, here the relation between stress and strain is considered. To evaluate stress, the compression force F was divided by the cross-sectional area of particle πd 2 p /4. For the definition of strain, stroke z is simply divided by d p . The measured stress-strain relations are shown in Fig. 9. As expected, any systematic size dependence cannot be observed. As can also be seen, the linear stress-strain relation continues upto about 5% strain. Then, the fracturing could occur around 15 kPa stress. It should be noted that the dust particles are basically plastic than elastic. Thus, the plastic deformation is induced even in the initial linear compression regime. Anyway, we can roughly capture the mechanical properties of dust particles from this experiment.
FIG. 9. Stress-strain relations for single dust particle compression tests. Particle size (diameter) was measured with a vernier caliper. Color indicates the particle size as labeled in the legend. Compression rate is basically fixed to v = 1.7 µm s −1 except for 6.0 mm size particle case (v = 17 µm s −1 ).

Stress comparison with periodic undulation
In Fig. 10, stress of the periodic undulation part for dust-particles compression, (F − F )/S computed from the data shown in Fig. 3, is shown as a function of strain per grain diameter, z/d (with d = 1 mm). The undulation amplitude reaches about the mechanical strength of dust particles around 10d compression. We are not sure whether the upper limit of the amplitude of periodic undulation is restricted by the strength of particles, or not. Much more compression is necessary to reveal the behaviors of highly compressed dust-particles column. However, the early stage of compression cannot be resolved when we focus on the highly compressed stage by using a large capacity load cell. In this study, we have focussed on the early stage (lightly loaded state) of compression of a granular column. The heavy compression test is a possible future work.  Fig. 3(a). Definition of strain is same as the single particle compression.
be expressed as dF/dz = αF/z. Thus, K for power-law form (dust-particles compression) becomes For exponential growth of glass-beads compression, the corresponding differential equation is written as dF/dz = F/z G . Thus, K should be Here we consider the compression with ideal gas to develop the physical image of bulk modulus formulation. Adiabatic ideal gas satisfying pV γ = const. yields K gas = γp, where specific heat ratio γ is constant. By comparing this relation with Eqs. (B2) and (B3), it can be said that the granular bulk modulus is sensitive to its column height (volume) than gas.

Compression index
Next, we introduce the compression index which has been used to characterize soil-compression characteristics [15]. Void ratio e is defined by the ratio between void volume V void and solid volume V solid , e = V void /V solid = (1/φ bulk ) − 1. Then, the compression index C c is defined by the relation between e and p as Specifically, the slope in e vs. ln p corresponds to C c . The corresponding plot is shown in Fig. 11. One can confirm the wider linearity in Fig. 11(b) than (a). This is in fact because Eq. (B4) is consistent with Eq. (B3) as explained below. By considering the above-mentioned relations and φ bulk = M tot /ρ true V , where the total mass M tot and true density ρ true are constant, bulk modulus becomes This relation provides a physical meaning of the length scale z G . In other words, the exponential form of F (z) means that the compression index is a constant. For dust-particles compression, linear region is restricted as shown in Fig. 11(a). This is the reason why the power law fits better than the exponential form for dust particles. The crossover point (∼ 10 kPa) roughly agrees with the strength of a single particle. The substantial difference between glass beads and dust particles is deformability.
In the glass-beads column, contact-network rearrangement is the only way to compress the layer. That is, the glass-beads column is almost in the jammed state even in its initial state. In such a jammed state, exponential growth of compression force is observed. Indeed, soil is usually in such a jammed state. If there is enough space to be compressed, power-law form can practically be applied because the constant void compressibility is reasonable in this state. Therefore, we expect the dust-particles layer should also obey the exponential form once it is significantly compressed (at least around jamming density). In this experiment, however, we cannot compress the dust-particles column to this state due to the loadcell capacity. Such a high compression test is required to conclude this expectation.  Fig. 3(b) and Fig. 5(b). Dotted lines are displayed as guides to the eye. The steep angle in small p of (b) comes from the subtlety in the initial stage of glass-beads compression (Fig. 5).
FIG. 13. Force curve of the compression test of monomers. Tiny (5 µm) glass beads were used as monomer particles (φ bulk 0.3). Monomers are simply poured into the cylindrical cell without granulation process. Neither periodic undulation nor sudden force drops are observed during the compression. The inset shows the log-linear plot of the same data.