## Abstract

Upon compression, some soft granular crystals undergo pattern transformation. Recent studies have unveiled that the underlying mechanism of this transformation is closely tied to microscopic instability, resulting in symmetry breaking. This intriguing phenomenon gives rise to unconventional mechanical properties in the granular crystals, paving the way for potential metamaterial application. However, no consistent approach has been reported for studying other unexplored transformable granular crystals. In this study, we present a systematic approach to identify a new set of pattern-transformable diatomic granular crystals induced by microscopic instability. After identifying the kinematic constraints for diatomic soft granular crystals, we have generated a list of feasible particle arrangements for instability-induced pattern transformation under compression. Instead of computationally intensive finite element models (FEMs) with continuum elements, we adopt a simplified mass-spring model derived from granular contact networks to efficiently evaluate these feasible particle arrangements for pattern transformation. Our numerical analysis encompasses quasi-static analysis and microscopic/macroscopic instability analyses within the framework of linear perturbation. Subsequently, the pattern transformation of the identified particle arrangements is confirmed through quasi-static analyses employing detailed finite element (FE) simulations with continuum elements. Additional numerical simulations with continuum elements reveal that the pattern transformations of particle arrangements are significantly influenced by the initial void volume and some transformed granular crystals may exhibit strong low-frequency directional phononic band-gaps, which were not observed in the initial granular crystals.

## 1 Introduction

Granular crystals are composed of a tightly packed collection of particles which interact by exchanging forces and moments through their geometric contacts [1]. Highly nonlinear evolution in their contacts enables granular crystals to exhibit unusual wave phenomenons including solitary waves [2], shock waves [3], and discrete breathers [4]. Their uniqueness has promoted granular crystals as promising metamaterials opening up a wide range of applications such as impact mitigation [5], energy absorption [6], and energy harvesting [7]. Furthermore, researchers have begun to investigate the tunability of their wave characteristics, reporting potential applications in acoustics [8], optics [9], and electronics [10]. Their utility extends to non-destructive evaluation [11], seismic protection [12], and structural health monitoring [13]. Additionally, a strategy has been proposed in Ref. [14] to create three-dimensional (3D) architectured materials with exceptional properties and recyclability using granular crystals, thereby expanding their potential applications.

It is well known that phononic band structures are affected by particle arrangements of granular crystals [15] and can be tuned by the transformation in their patterns [16–18]. Some soft diatomic two-dimensional (2D) granular crystals are reported to exhibit compression-driven pattern transformations [19,20], which can be viewed as instability-induced symmetry breaking. For instance, compression rearranges an initial square granular crystal into a new pattern with assembled small particles as shown in Fig. 1. During this transition, the unit-cell size increases from its primitive one (containing one large particle and one small particle, delineated by the solid line in Fig. 1(a)) to an enlarged one (containing two large particles and two small particles, delineated by the solid line in Fig. 1(b)), so the number of particles in the unit-cell becomes doubled.^{2} Recently, employing a simplified spring-mass model derived from granular contact networks [21] have identified that the underlying mechanism of this pattern transformation in soft granular crystals is closely related to microscopic instability.

The objective of this study is to propose a systematic approach for constructing a new set of instability-induced, pattern-transformable diatomic granular crystals by leveraging the recently unveiled mechanism of pattern transformation in soft granular crystals. Considering the instability arising from kinematic constraints between particles, we first present a catalog of viable particle arrangements for pattern transformation. Subsequently, these particle arrangements are examined for pattern transformations through a series of numerical analyses. Simple mass-spring models are employed for instability analyses in the linear perturbation framework and quasi-static analyses with infinitely periodic configurations. Furthermore, we confirm their pattern transformations by conducting quasi-static analyses using detailed finite element (FE) simulations with continuum elements. Additionally, we discuss how the radius ratio affects the pattern transformation in the considered diatomic granular crystals. As an example of the potential applications of pattern-transformable granular crystals, lastly, we explore their phononic band structures, which exhibit tunability during deformation.

## 2 Method

This section outlines the principles governing particle arrangements that result in instability-induced pattern transformations under compression. Additionally, it reports on a series of numerical analyses employed for systematic exploration.

### 2.1 Kinematic Constraints.

There could be various granular particle configurations exhibiting instability-induced symmetry breaking behavior, which subsequently results in pattern transformations under compression. In this study, however, we introduce a systematic procedure to identify a new set of instability-induced *pattern-transformable* granular crystals. The approach involves constraining kinematics within diatomic configurations, but it is important to note that a similar procedure can be applied to triatomic or tetratomic configurations as well.

For a given granular crystal, first, we establish the corresponding contact network where a nodal joint is positioned at the center of the particle, and a link is assigned between two contiguous particles. Examples of contact networks are depicted in Fig. 2. Now, motivated by the previous study revealing the mechanism of the instability-induced pattern transformations in soft granular crystals [22], we start by specifying a set of kinematic constraints essential for this investigation: (*KC1*) the granular contact network must include connected links oriented in the direction of compression, and (*KC2*) each nodal joint must be connected to at least three connecting particle links. Notably, the kinematic constraint *KC1* forms a column configuration crucial for inducing instability upon compression, forming the pivotal constraint for exploring instability-induced pattern transformation. Furthermore, the kinematic constraint *KC2* is imperative for initiating symmetry breaking along the compressive loading direction, leading to instability-induced pattern transformation. Otherwise, nodal joints with only two connecting particle links could trigger symmetry breaking, leading to an arbitrary buckled shape under compression. These two kinematic constraints are vital for inducing the buckling of vertical column chains, which are relatively unaffected by geometric imperfections. Note that while these two kinematic constraints are necessary conditions for pattern transformations of granular crystals, they are not sufficient conditions.

### 2.2 Numerical Models.

In addition to the aforementioned kinematic constraints, the presence of voids among particles plays a critical role in triggering the buckling of vertical column chains. Consequently, particle arrangements that conform to the specified kinematic constraints may not necessarily guarantee a pattern transformation in granular crystals. The void ratio in the granular crystals is directly correlated with the radius ratio of constituent particles. Therefore, it becomes imperative to numerically investigate the nonlinear behavior of each particle arrangement under compression across a broad spectrum of radius ratios.

In order to capture the pattern transformation of soft granular crystals, the prior studies adopted the discrete element method (DEM) [19]. However, DEM was not adequate for unveiling the mechanism of pattern transformation in granular crystals or understanding why these transformations occur in different configurations [20]. Because investigating the underlying mechanisms requires adopting a sequential linear perturbation approach, such as eigenvalue analyses in an incremental fashion, DEM alone was insufficient. Instead, considering DEM simulation results, we identify the particle contact network incrementally. Until a new contact network forms upon loading, we approximate the particle contact behavior using a simple mass–spring model within the finite element model (FEM) framework. Once the contact network changes, we introduce a new mass-spring model sequentially. For the identified primitive unit-cell from the contact network, we perform a series of numerical simulations, including quasi-static analysis and microscopic/macroscopic instability analyses. This incremental study reveals that microscopic instability may lead to the pattern transformation of the granular crystals under compression [21]. We adopt the same numerical procedure in this study to identify which granular configurations are affected by microscopic instability upon incremental loading. Furthermore, we confirm their pattern transformations by conducting quasi-static analyses using detailed FE simulations with continuum elements. Additionally, the continuum element models are also employed to investigate the effects of the radius ratio and the tunability of the phononic dispersion relations.

#### 2.2.1 Material Properties and Model for Continuum Elements.

Properties | Hard particle | Soft particle |
---|---|---|

Density, $\rho (kg/mm3)$ | $2.15\xd710\u22126$ | $1.05\xd710\u22126$ |

Young’s modulus, $E(MPa)$ | 1000 | 0.360 |

Shear modulus, $G(MPa)$ | 337 | 0.121 |

Bulk modulus, $K(MPa)$ | 16,700 | 6.00 |

Dilatational wave speed, $cd=G[E\u22124G]/[\rho [E\u22123G]](mm/s)$ | $2820\xd7103$ | $76.6\xd7103$ |

Properties | Hard particle | Soft particle |
---|---|---|

Density, $\rho (kg/mm3)$ | $2.15\xd710\u22126$ | $1.05\xd710\u22126$ |

Young’s modulus, $E(MPa)$ | 1000 | 0.360 |

Shear modulus, $G(MPa)$ | 337 | 0.121 |

Bulk modulus, $K(MPa)$ | 16,700 | 6.00 |

Dilatational wave speed, $cd=G[E\u22124G]/[\rho [E\u22123G]](mm/s)$ | $2820\xd7103$ | $76.6\xd7103$ |

Every FE simulation with continuum elements was performed under plane stress conditions using a 4-node bilinear element with reduced integration and hourglass control (CPS4R in abaqus). In addition, a mesh refinement study was performed to ensure the convergence of FE simulations, resulting in the final sweeping seed size of 0.6 mm or smaller. FE analysis employs the Coulomb friction model with the friction coefficient of $\mu =0.01$ to model contact conditions.

#### 2.2.2 Mass-Spring Model for Nonlinear Contact Forces.

It is worth noting that the exploration of pattern transformation of soft granular crystals is a computationally expensive process. To circumvent the high computational expenses of continuum element simulations, we employ a simple mass-spring model (see Fig. 3), as presented in Refs. [19,22], to assess the potential for pattern transformation. Jain and Shim [22] showed that the particle contact motions can be effectively represented by a simple mass-spring model. Here, we provide only a brief description of the mass-spring model adopted in this study. The employed mass-spring model is composed of three linear springs (i.e., a normal spring $kn$, a tangential spring $kt$, and a rotational spring $kr$), which capture the qualitative mechanical behavior and enlarged periodicity of granular crystals during deformation.

*Model Parameters*. The structural element consists of nodes at each particle center with the three degrees-of-freedom, $un$, $ut$, and $\theta $, denoting the normal displacement, the tangential displacement, and the rotational angle, respectively (see Fig. 3(b)). The element comprises three linear springs, i.e., a normal spring $kn$, a tangential spring $kt$, and a rotational spring $kr$. The normal spring $kn$ is modeled by

Type | Parameter | Hard–hard | Hard–soft | Soft–soft |
---|---|---|---|---|

Particle arrangement (Figs. 5(a)–5(d) and 5(f)) | $kn(N/mm)$ | 3468 | 2.52 | 1.35 |

$h(mm)$ | 2.50 | 3.33 | 5.00 | |

$\beta $ | $10\u22122$ | $10\u22122$ | $10\u22122$ | |

$\gamma (mm)$ | 1 | 1 | 1 | |

Particle arrangement (Fig. 5(e)) | $kn(N/mm)$ | 3468 | 2.52 | 1.35 |

$h(mm)$ | 3.00 | 3.75 | 5.00 | |

$\beta $ | $10\u22122$ | $10\u22122$ | $10\u22122$ | |

$\gamma (mm)$ | 1 | 1 | 1 |

Type | Parameter | Hard–hard | Hard–soft | Soft–soft |
---|---|---|---|---|

Particle arrangement (Figs. 5(a)–5(d) and 5(f)) | $kn(N/mm)$ | 3468 | 2.52 | 1.35 |

$h(mm)$ | 2.50 | 3.33 | 5.00 | |

$\beta $ | $10\u22122$ | $10\u22122$ | $10\u22122$ | |

$\gamma (mm)$ | 1 | 1 | 1 | |

Particle arrangement (Fig. 5(e)) | $kn(N/mm)$ | 3468 | 2.52 | 1.35 |

$h(mm)$ | 3.00 | 3.75 | 5.00 | |

$\beta $ | $10\u22122$ | $10\u22122$ | $10\u22122$ | |

$\gamma (mm)$ | 1 | 1 | 1 |

*Mass and Stiffness Matrices*. The mass and stiffness matrices of the considered structural element can be derived by applying Euler–Lagrange equations. The kinetic energy $K$ and potential energy $\Pi $ of the element can be written as

Type | Properties | Hard particle | Soft particle |
---|---|---|---|

Particle arrangement (Figs. 5(a)–5(d) and 5(f)) | Radius, $r(mm)$ | 2.50 | 5.00 |

Linear mass, $m(kg)$ | $9.32\xd710\u22123$ | $1.96\xd710\u22124$ | |

Angular mass, $\iota (kgmm2)$ | $1.94\xd710\u22124$ | $1.64\xd710\u22125$ | |

Particle arrangement (Fig. 5(e)) | Radius, $r(mm)$ | 3.00 | 5.00 |

Linear mass, $m(kg)$ | $1.34\xd710\u22122$ | $1.96\xd710\u22124$ | |

Angular mass, $\iota (kgmm2)$ | $4.02\xd710\u22124$ | $1.64\xd710\u22125$ |

Type | Properties | Hard particle | Soft particle |
---|---|---|---|

Particle arrangement (Figs. 5(a)–5(d) and 5(f)) | Radius, $r(mm)$ | 2.50 | 5.00 |

Linear mass, $m(kg)$ | $9.32\xd710\u22123$ | $1.96\xd710\u22124$ | |

Angular mass, $\iota (kgmm2)$ | $1.94\xd710\u22124$ | $1.64\xd710\u22125$ | |

Particle arrangement (Fig. 5(e)) | Radius, $r(mm)$ | 3.00 | 5.00 |

Linear mass, $m(kg)$ | $1.34\xd710\u22122$ | $1.96\xd710\u22124$ | |

Angular mass, $\iota (kgmm2)$ | $4.02\xd710\u22124$ | $1.64\xd710\u22125$ |

### 2.3 Quasi-static Analysis With Infinitely Periodic Structures.

### 2.4 Instability Analyses in Linear Perturbation Framework.

By employing the mass-spring model within particle contact networks in an infinitely periodic setting, we examine both the microscopic instabilities with wavelengths on the order of the unit-cell size and the macroscopic instabilities with wavelengths much larger than the unit-cell size [26]. As our study specifically addresses compression-driven symmetry breaking in granular crystals, we focus on mechanical instabilities that lead to a unit-cell change in the periodicity of a particle contact network during deformation (i.e., pattern transformation induced by microscopic instability). The examination of an infinitely periodic structure is effectively performed by monitoring changes in unit-cell periodicity induced by instability.

Operationally, for a given particle contact network, the mass-spring model is utilized to conduct quasi-static analysis under incremental compressive loads in an infinitely periodic setting. Then, two types of instability analyses are executed. The first *microscopic* instability of various unit-cell periodicities is investigated by detecting the strain at which the smallest eigenfrequency associated with a non-trivial eigenmode becomes zero [27]. The minimum strain at which this occurs is considered the critical strain for microscopic instability in the finite periodic structure, and the corresponding unit-cell size determines the unit-cell periodicity after the microscopic instability. On the other hand, the second *macroscopic* instability is examined by checking the loss of ellipticity condition of the homogenized tangent modulus [28]. The loss of ellipticity condition is determined using a single primitive unit-cell with the conventional periodic boundary conditions. The critical strain, violating the positive definite condition of the elasticity tensor, is then identified as the critical strain for macroscopic instability. The aforementioned three analyses are conducted within the linear perturbation framework.

For a potential pattern-transformable granular crystal, we opt for particle contact networks where the critical strain for the microscopic instability is smaller than the critical strain for the macroscopic instability. Following the completion of these instability analyses in the linear perturbation framework, we perform the quasi-static analysis on the enlarged unit-cell periodicity associated with the microscopic instability to identify the complete pattern transformation. The detailed numerical procedure can be found in Refs. [26,27]. However, for the sake of completeness, we present a concise implementation scheme of the quasi-static and instability analyses.

#### 2.4.1 Microscopic Instability Analysis.

#### 2.4.2 Macroscopic Instability Analysis.

In practice, the loss of ellipticity of the homogenized tangent modulus $L\xaf$ is examined by adopting a primitive unit-cell with the conventional periodic boundary conditions (12). We apply four independent linear perturbations for the homogenized deformation gradient $F\xaf\u2218$ at every deformation increment in the 2D setting. Subsequently, we extract the corresponding homogenized first Piola–Kirchhoff stress components $P\xaf\u2218$. Finally, we check all components of $L\xaf$ from (17) by exploring $m$ and $N$ at every $\pi /360$ radian increment. For each combination of $m$ and $N$, we assess the loss of ellipticity condition (16).

### 2.5 Phononic Dispersion Analysis.

Subsequently, we perform comprehensive FE quasi-static analyses using continuum elements for particles to confirm that the anticipated instability-induced pattern transformations occur in all the identified particle arrangements. Additionally, we calculate the phononic dispersion relation of the examined particle arrangements using finite element models with continuum elements, as the mass-spring model for granular crystals may provide inaccurate estimates of the phononic band structure for wave motion [18]. Jain and Shim [18] outlined a comprehensive procedure for obtaining the phononic dispersion relation from granular crystals. Here, we offer a condensed summary of their approach for completeness.

Additionally, the wave motion within periodic structures must also satisfy the Bloch periodic condition for $u\u2218\kappa $ shown in (13). Then, one can obtain the phononic dispersion relation $\omega =\omega (\kappa )$ of the considered granular crystal by solving the eigenvalue problem (18) with (13) for a given wavevector $\kappa $.

Operationally, we adopt the same procedure described in Sec. 2.4.1 to handle the complex values in the incremental displacement $u\u2218\kappa $ in (18). Moreover, we only consider wavevectors located on the perimeter of the irreducible Brillouin zone along with its evolution during deformation to calculate the phononic dispersion relation [26]. The procedure to identify the irreducible Brillouin zone of some granular crystals and a 2D bilayered composite is shown in Fig. 4. Finally, the natural frequencies of the considered FEM unit-cell model are calculated by performing the eigenfrequency analysis using the abaqus keyword, *FREQUENCY for each considered wave vector $\kappa $.

## 3 Results

This section presents a catalog of viable particle arrangements that lead to pattern transformations based on the kinematic constraints *KC1* and *KC2*. Then, the confirmation of pattern-transformable granular crystals is achieved through a series of numerical analyses involving both simple mass-spring models and detailed FE models with continuum elements.

### 3.1 Particle Arrangements Using Kinematic Constraints.

Based on the diatomic compact packing [30,31], we provide a list of various circle arrangements satisfying the kinematic constraints *KC1* and *KC2*, as illustrated in Fig. 5. A vertical column chain satisfying *KC1* can consist of different combinations of large and small particles. However, the kinematic constraint associated with three adjacent particle contacts *KC2* restricts the layout of vertical column chains in the 2D space. Thus, we have only four types of vertical column chains: (1) large particles shown in Fig. 5(a), (2) alternating one large particle and one small particle in Fig. 5(d), (3) alternating one large particle and two small particles in Fig. 5(e), and (4) alternating two large particles and two small particles in Fig. 5(f). No column chains other than these four can satisfy both kinematic constraints *KC1* and *KC2* in the 2D space. The second column of Fig. 5 also displays the particle contact networks associated with the particle arrangements and the corresponding primitive unit-cell. These four 2D particle arrangements shown in Figs. 5(a), 5(d), 5(e), and 5(f) serve as the basic configurations for the subsequent particle arrangements.

Building upon those four basic configurations, we can further construct feasible particle arrangements for instability-induced pattern transformation by removing a set of small particles while maintaining the two kinematic constraints. Specifically, the particle arrangement in Fig. 5(b) is obtained by eliminating the alternating horizontal lines of small particles from the particle arrangement in Fig. 5(a). The elimination of the alternating diagonal lines (or vertical lines) of small particles from Fig. 5(a) results in the particle arrangement in Fig. 5(b′) (or Fig. 5(b)″). Furthermore, the elimination of the second set of small particles generates the arrangements shown in Figs. 5(c) and 5(c)′ in a similar fashion. This small particle elimination procedure can also be applied to the other three basic structures shown in Figs. 5(d), 5(e), and 5(f). The elimination of the first set of small particles from these three structures results in Figs. 5(d)′, 5(e)′, and 5(f)′, respectively. It is noteworthy that the elimination of the second set of small particles from Figs. 5(d), 5(e), and 5(f) violates the constraint *KC1*.

Essentially, by beginning with the kinematic constraints identified in this study and subsequently applying the small particle elimination, we generate a list of 13 particle arrangements in Fig. 5. These arrangements fulfill the essential kinematic conditions required for compression-induced pattern transformation in granular crystals.

### 3.2 Analyses With Mass-Spring Models.

For all the 13 particle arrangements in Fig. 5, we performed two types of effective numerical analyses using the simple mass-spring models: (1) instability analyses (i.e., microscopic and macroscopic instabilities) in the linear perturbation framework and (2) quasi-static analysis. From those series of intensive numerical analyses, we find that the most distinctive pattern transformations are observed in the configurations A, B, C, D, and E. In other words, as compression increases, other granular configurations have multiple competing unit-cells for pattern transformation in deformed shapes, illustrating a close race between multiple different unit-cells in microscopic instability. Consequently, their resulting patterns under compression are not as clear as those observed in the configurations A, B, C, D, and E. Accordingly, this section presents only the numerical results from the five granular configurations A, B, C, D, and E showing conspicuous pattern transformations.

#### 3.2.1 Instability Analyses in the Linear Perturbation Framework.

A primitive unit-cell in the initial (undeformed) particle contact network is provided in the first column of Fig. 6. The second column of Fig. 6 illustrates two pieces of information regarding vertical compression: (1) the evolution of the eigenfrequencies obtained from the microscopic instability analyses (lines with markers for the left vertical axis of the graph) and (2) the positive definiteness condition of the homogenized fourth-order mixed elasticity tensor $L\xaf$ related to macroscopic instability (black solid line without markers for the right vertical axis of the graph). Note that the third column of Fig. 6 presents a close-up graph of its second column.

Consider the particle arrangement shown in Fig. 6(a). Its eigenfrequency associated with the enlarged $(1,2)$-periodicity becomes zero at the engineering strain level of $\epsilon eng=0.47$. However, the corresponding positive definite condition of $L\xaf$ is not violated (see the graph in Fig. 6(a)). Consequently, this particle arrangement is prone to deform following microscopic instability rather than macroscopic instability. In addition, the fourth and fifth columns of Fig. 6(a) illustrate the unit-cell corresponding to the enlarged $(1,2)$-periodicity after microscopic instability and its eigenmode shape, respectively.

From the numerical analyses of 13 granular patterns shown in Fig. 5, we excluded the eight patterns that violated the loss of ellipticity condition. Ultimately, we identified and presented a potential list of five granular patterns (A–E) in Fig. 6 as pattern-transformable granular crystals. In every particle arrangement depicted in Fig. 6, microscopic instability precedes macroscopic instability upon vertical compression, as critical strain for microscopic instability is significantly smaller than for macroscopic instability. Thus, the macroscopic instability does not play a role in the particle arrangements identified in this study. The analyses of microscopic instability using the simple mass-spring models indicate that the unit-cell periodicity evolves into either $(1,2)$- or $(2,1)$-periodicity depending on the type of particle arrangement. Specifically, granular patterns A, B, and C exhibit the enlarged $(1,2)$-periodicity, whereas granular patterns D and E display the enlarged $(2,1)$-periodicity. On the other hand, we also observe that microscopic critical strains regarding $(1,2)$- and $(1,3)$-periodicities are very close in granular patterns B and C (see Figs. 6(b) and 6(c)). This proximity in critical strains suggests a potential competition between those two enlarged unit-cell microscopic instabilities. However, note that the mass-spring models have a tendency to overestimate the contact stiffness behavior between particles as deformation increases, with the initial contact network never breaking. Therefore, the actual microscopic instability may favor the smaller enlarged unit-cell periodicity, especially given the prevalence of imperfections within the granular crystals. Consequently, FE simulations with continuum elements are conducted subsequently to provide further insights.

#### 3.2.2 Quasi-static Analyses With Infinitely Periodic Structures.

After completing the series of instability analyses, we executed quasi-static analyses under vertical compression using the mass-spring models for the configurations A, B, C, D, and E from Fig. 5. The stress–strain relations of the particle contact network constructed using the mass-spring model based on two different boundary conditions are depicted in the first column of Fig. 7. The first red line with markers in the graphs denotes the result from the infinitely periodic structure simulated by enlarged unit-cells illustrated in the column of Fig. 7. Conversely, the second green dotted line in the graphs indicates the result from a $6\xd78$ finite size model structure partially depicted in the third column of Fig. 7. These two stress–strain curves exhibit similarity, suggesting that the unit-cell simulation for infinitely periodic structure effectively captures the nonlinear behavior of the finite size model structure. It is noteworthy that a sudden stiffness change in the stress–strain curves is observed at the strain level associated with the onset of microscopic instability shown in Fig. 6. This nonlinear behavior of stress–strain curves confirms that the identified five granular patterns are clearly induced by instability. Furthermore, the deformed shapes in the second and third columns of Fig. 7 illustrate a symmetry breaking phenomenon initiated by the buckling of the vertical column chains. Numerical analyses with the mass-spring models effectively identify potential particle arrangements leading to the pattern transformation induced by microscopic instability.

### 3.3 Analyses With Continuum Element Models.

In this section, we verify the pattern-transformable behavior of the five distinctive particle arrangements (i.e., configurations A, B, C, D, and E in Fig. 5) by conducting quasi-static analysis under compression using FE simulations with continuum elements. Additionally, to explore the effect of pattern transformation, we utilize the same FE models with continuum elements for phononic dispersion analyses in the linear perturbation framework.

#### 3.3.1 Quasi-static Analyses With Infinitely Periodic Structures.

The quasi-static FE analysis of the five particle arrangements is presented in Fig. 8. As mentioned in Sec. 2.2, the detailed mechanical properties of the materials under consideration can be found in Sec. 2.2.1, and the procedure for quasi-static analysis in infinitely periodic structures is outlined in Sec.2.3.

Here, we confirm the qualitative agreement between the deformed shapes obtained from the mass-spring model simulations (see Fig. 7) and the ones from the FE simulations (see Fig. 8). At the onset of vertical compression, particle contacts are firmly established within the FE model. As compression progresses, stress accumulates with its distribution closely resembling the initial contact network of the corresponding granular crystals depicted in Fig. 5. Specifically, high stress become prominent along the vertical particle columns as compression increases. With further compression, the vertical particle columns begin to buckle, breaking the symmetry of the initial particle arrangements (see the second column of Fig. 8). The buckling process opens up voids among particles, which are rapidly filled by neighboring particles, initiating pattern transformation. Subsequent compressive force only amplifies the newly formed pattern in granular crystals (see the third column of Fig. 8).

Specifically, the initial granular configuration A has a uniform distribution of smaller particles (see Fig. 8(a)). However, as compression increases, two small particles coalesce, giving rise to a new pattern identical to the square granular crystal previously reported in Ref. [19]. To the best of the author’s knowledge, the remaining four configurations B, C, D, and E have not been previously documented as pattern-transformable granular crystals. Notably, in the initial granular configuration D, particles of the same size are completely separated as shown in Fig. 8(d). In other words, each small particle is surrounded by four large particles only, and each large particle is surrounded by four small particles only. However, with increasing compression, all the large particles become contiguous in a zigzag vertical direction. A similar pattern transformation is also observed in granular configuration E (see Fig. 8(e)). The compressed granular configurations D and E can be approximately viewed as infinitely periodic bilayered composites since the large and the small particles possess different mechanical properties.

Essentially, these quasi-static FE analysis results demonstrate that the prescribed kinematic constraints (*KC1* and *KC2*) and the simple mass-spring models for the contact network offer a systematic approach to identifying a new set of pattern-transformable diatomic granular crystals. In this study, we select the simplest diatomic configuration connected to the previously identified pattern-transformable granular crystals [19,20]. However, the presented systematic approach can be further applied to other types of configurations such as triatomic or tetratomic granular configurations.

#### 3.3.2 Phononic Dispersion Analyses.

Based on this newly identified set of pattern-transformable granular crystals, we further investigate how pattern transformations affect the phononic band structures of granular crystals. Two representative particle arrangements (i.e., one in Fig. 5(d) (or Fig. 8(d)), and the other in Fig. 5(e) (or Fig. 8(e))) are selected for this purpose due to their approximate bilayered feature after the pattern transformation. For each particle arrangement, we explore two different particle radius ratios: one granular crystal not showing pattern transformation and another one exhibiting pattern transformation. Note that the pattern transformation may not occur even for the granular configurations A, B, C, D, and E if the initial granular configuration has very few voids between the particles. Further discussion relating to the radius ratio effect is presented in Sec. 4.1.

First, consider the granular crystal D shown in Fig. 8(d), with its phononic band-structure evolutions summarized in Fig. 9. We choose the radius ratio of $rsmall:rlarge=1:2.22$ for a non-pattern-transformable granular crystal (see Fig. 9(a)) and the radius ratio of $rsmall:rlarge=1:1.66$ for pattern-transformable granular crystal (see Fig. 9(b)). Their phononic band structures are presented at three different strain levels (i.e., $\epsilon eng=0.01,0.15,0.25$). At each strain level, we provide three sub-figures: (1) a unit-cell FE model, (2) the corresponding irreducible Brillouin zone, and (3) its phononic dispersion graph. Note that a phononic dispersion graph is presented along the perimeter of the irreducible Brillouin zone (i.e., enclosed area by O–M–X–G–O) in the reciprocal space, and the schematic procedure to determine the irreducible Brillouin zone is illustrated in Fig. 4. The orange-shaded region in the phononic dispersion graph indicates complete phononic frequency band-gaps, in which the propagation of mechanical waves is impermissible for any wavevector directions. In addition, the green-shaded region indicates directional phononic band-gaps along the X–G path in the reciprocal space. Now, at the strain level $\epsilon eng=0.01$, overall particle contacts start to develop in both non-pattern-transformable and pattern-transformable granular crystals (see the FE models in the lattice space in Figs. 9(a-1) and 9(b-1), respectively). At this small strain level, contact areas between adjacent particles are small, and the corresponding contact stiffness is also very low. Thus, this low contact stiffness results in broad complete phononic band-gaps for both granular crystals because their contact network is similar at these small strain levels. It is also well known that uncompressed granular crystals exhibit extremely slow wave speed and huge phononic band-gaps due to negligible linear contact stiffness, a phenomenon referred to as sonic vacuum [32,33]. On the other hand, as compression increases, the particle contact area significantly enlarges. Eventually, an initial loosely contiguous granular crystal with a high void volume becomes approximately a continuum solid with a low void volume. At the strain levels of $\epsilon eng=0.15$, the initially broad complete phononic band-gaps mostly disappear regardless of the radius ratio. However, the pattern-transformable granular crystal exhibits strong directional phononic band-gaps along the X–G path (see Fig. 9(b-3)), whereas the non-pattern-transformable granular crystal does not (see Fig. 9(a-3)). At the strain level of $\epsilon eng=0.25$, we observe these directional phononic band-gaps affected by pattern transformation.

Second, this interesting phenomenon regarding directional phononic band-gaps is also observed in the particle arrangement shown in Fig. 5(e). The radius ratio of $rsmall:rlarge=1:3.33$ is selected for a non-pattern-transformable granular crystal (see Fig. 10(a)) and $rsmall:rlarge=1:2.00$ is for pattern-transformable granular crystal (see Fig. 10(b)). As shown in Figs. 10(a-1) and 10(b-1), there is not much difference in the phononic band structures at $\epsilon eng=0.01$. However, as compression increases, the pattern-transformable granular crystal exhibits strong directional phononic band-gaps along the X–G path (see Fig. 10(b-3)), whereas the non-pattern-transformable granular crystals show negligible directional phononic band-gaps (see Fig. 10(a-3)).

We propose that these intriguing directional phononic band-gaps exhibited by the granular crystals are associated with their unique bilayer structure following the pattern transformation. Further discussion is presented in Sec. 4.2.

## 4 Discussion

The findings of this research affirm the effectiveness of the proposed systematic approach for constructing a new set of pattern-transformable diatomic granular crystals, which also offer tunability in their phononic band structures. This section further discusses the influence of radius ratio and the directional phononic band-gaps observed in pattern-transformable granular crystals.

### 4.1 Effect of Radius Ratio on Pattern Transformation.

Throughout this investigation, we have discovered that the radius ratio between large and small particles plays a significant role in constructing the list of the feasible particle arrangements for instability-induced pattern transformations, as depicted in Fig. 5. Depending on the range of viable radius ratios, we can categorize the particle arrangements into two groups. The first group comprises arrangements with only one feasible radius ratio for a given granular configuration. Specifically, the particle arrangements shown in Figs. 5*c*′, *d*′, *e*′, and *f*′, the only corresponding feasible radius ratios are $rsmall/rlarge=2$, $2$, $(17\u22123)/4$, and $(5\u22121)/2$, respectively. Conversely, the second group encompasses arrangements with infinitely many viable radius ratios for a given granular configuration. In other words, all other particle arrangements (i.e., particle arrangements shown in Figs. 5(a)–5(f), (*b*′), (*b*$\u2033$), and (*c*$\u2033$)) possess a range of viable radius ratios for their undeformed configuration setting.

The pattern transformations of the particle arrangements are significantly influenced by the initial void volume. When the void volume is too low, particles lack sufficient space to maneuver even after the microscopic instability of the vertical particle columns occurs. On the other hand, excessive void volume allows the particles to move in an uncontrolled manner following the initiation of microscopic instability in vertical particle columns. It is worth noting that, for the particle arrangements permitting only one radius ratio (i.e., granular patterns $C\u2032$, $D\u2032$, $E\u2032$, and $F\u2032$), the void volume ratio is fixed. Consequently, this constraint eliminates the control parameter necessary to explore pattern transformations for a given particle arrangement. Numerical analysis reveals that none of the particle arrangements restricted by a single radius ratio exhibit noticeable pattern transformations. In contrast, distinctive pattern transformations occur in the particle arrangements feasible within a range of radius ratios, particularly those with relatively high void volume. The results depicted in Figs. 8–10 are derived from the findings of these numerical studies.

### 4.2 Directional Phononic Band-Gaps in Pattern-Transformable Granular Crystals.

For the pattern-transformable granular crystals shown in Figs. 8(d) and 8(e), pronounced directional phononic band-gaps along the X–G path in the reciprocal lattice space are observed at high levels of compressive strain, as illustrated in Figs. 9(b) and 10(b). Here, we hypothesize that these significant directional phononic band-gaps under high strain conditions are linked to the deformed configuration of the granular crystals, which can be approximately viewed as infinitely periodic bilayered composite. Figures 9(b-3) and 10(b-3) at $\epsilon eng=0.25$ illustrate that the assembly of small particles completely separates the vertical chains of large particles from adjacent chains of large particles. At these large strain levels, the particle contact areas become substantially large, resulting in linear mechanical behavior akin to that of an infinitely periodic bilayered composite. Subsequently, the wave characteristics of the diatomic granular crystal may behave similarly to those of an equivalent bilayered composite composed.

To test the hypothesis, we compute the analytical phononic dispersion relation of an infinitely periodic bilayered composite within a 2D configuration. This composite structure closely resembles the rotated square granular crystal D presented in Fig. 9. For the approximate bilayered composite, Fig. 11 illustrates a 2D unit-cell with a finite dimension of $a2$ based on the deformed shape of rotated square granular crystal D, along with its corresponding irreducible Brillouin zone. It is worth noting that in the conventional 1D setting of infinitely periodic bilayered composites, the dimension $a2$ should be infinitely small [34]. Detailed discussions on the fictitious modes associated with the finite dimension of $a2$ can be found in Ref. [35]. In this 2D unit-cell setting for an infinitely periodic bilayered composite, the analytic expression for wave motion in the sagittal plane [36] yields the phononic dispersion relation presented in Fig. 11. As anticipated, the strong directional phononic band-gaps along the X–G path from the rotated square granular crystals D (see Fig. 9) also appear in the phononic dispersion graph of the approximate bilayered composite as shown in Fig. 11.

This numerical finding underscores the mechanical tunability of the proposed pattern-transformable granular crystals. These crystals exhibit changes in their wave properties as they transition from a sparse initial configuration to a newly ordered compact packing under mechanical loading. Note that Shan et al. [22] demonstrated the phononic tunability of *continuum specimens* made from the same elastomers used in this study and showed that the phononic dispersion relations of these materials can be effectively modeled with elastic material models within the finite element framework. However, additional experimental studies following this research are necessary to validate the estimated phononic tunability of soft granular crystals, especially considering energy loss due to friction between particles and the viscoelastic properties of certain elastic materials.

In this research, as an example of the potential applications of pattern-transformable granular crystals, we chose to focus on their phononic band structures rather than other applications such as light switching [37]. Interestingly, because the microscopic instability introduces symmetry breaking in the primitive unit-cell, the Willis coupling [38,39] may appear during the pattern transformation. Given the emergence of the asymmetric micro-structures is the key feature of the instability-induced pattern transformation, the coupled pressure–strain and momentum–velocity relation could be an interesting topic for future study. In particular, in order to investigate the Willis effect, we could further explore the complex band gap structures near $\omega ~\u223c0.12$, which are highly affected by the pattern transformation (see Figs. 9 and 10). In future research, besides simply comparing the phononic band structures of granular crystals and their corresponding bilayered composites, we could further achieve a clearer description of the spatial dispersion effect by extending beyond the typical quasi-static approximation. This can be done using dynamic homogenization methods, such as the Willis approach [40] and the two-scale approach [41]. This elegant route, utilizing dynamic homogenization, could provide an in-depth understanding of granular crystals and a procedure to design Willis-type metamaterials in future research.

## 5 Conclusion

We have provided a systematic approach to identify a new set of pattern-transformable granular crystals whose motions are triggered by microscopic instability. After identifying the kinematic constraints for diatomic soft granular crystals to undergo pattern transformation, we have compiled a list of feasible particle arrangements that can transform under compression. Instead of computationally intensive FE models with continuum elements, we introduce a simplified mass-spring model derived from granular contact networks to efficiently evaluate these feasible particle arrangements for pattern transformation. Our numerical analysis encompasses quasi-static analysis and microscopic/macroscopic instability analyses within the framework of linear perturbation. Subsequently, the pattern transformation of the identified particle arrangements is confirmed through quasi-static analyses employing detailed FE simulations with continuum elements. Additional numerical simulations with continuum elements reveal that the pattern transformations of particle arrangements are significantly influenced by the initial void volume. Moreover, as an example of the potential applications of pattern-transformable granular crystals, their phononic dispersion relations are studied, demonstrating the tunability of phononic band structures during deformation. This numerical study implies that external loading could be applied to tune the phononic dispersion relation of soft granular crystals, but follow-up experimental studies are needed to validate the estimated phononic tunability of soft granular crystals. Furthermore, since the emergence of instability-induced asymmetric micro-structures is a key feature of instability-induced pattern transformation, the Willis coupling effect could be further examined to gain a deeper understanding of these granular crystals. Lastly, this study is based on the simplest diatomic configuration, but the presented systematic approach can be further applied to other types of configurations such as triatomic or tetratomic granular configurations. We anticipate that this pattern transformation phenomenon could also be activated by different driving forces if the constituents of granular crystals respond to external stimuli such as changes in temperature, pH, or water content.

## Footnote

An enlarged unit-cell periodicity is denoted by $(n1,n2)$-periodicity (or $n1\xd7n2$) with $ni$ indicating the number of the primitive units in the $ai$ direction ($i=1,2$). Thus, for the square granular crystal in Fig. 1, the enlarged unit-cell is indicated by $(n1,n2)=(1,2)$ whereas the primitive unit-cell is denoted by $(n1,n2)=(1,1)$.

## Acknowledgment

The information, data, or work presented herein was funded in part by the National Science Foundation Project under Award Number CMMI-2140224 and the Advanced Research Projects Agency-Energy (ARPA-E), U.S. Department of Energy, under Award Number DE-AR0001771. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

## Conflict of Interest

All the authors declare that there is no conflict of interest with any financial organization regarding the material discussed in the paper.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.