Vibrational properties of two-dimensional dimer packings near the jamming transition

TL;DR

This study investigates the vibrational and contact properties of two-dimensional dimer packings near the jamming transition, revealing isostaticity, characteristic vibrational modes, and scaling laws similar to spherical systems.

Contribution

It provides the first detailed analysis of vibrational density of states and contact number scaling in 2D dimer packings, including rotational rattlers and mode characterizations.

Findings

Dimers become isostatic at jamming after correcting for rotational rattlers.

Vibrational density of states shows two plateaus separated by a peak.

Characteristic frequencies follow critical scaling laws near the transition.

Abstract

Jammed particulate systems composed of various shapes of particles undergo the jamming transition as they are compressed or decompressed. To date, sphere packings have been extensively studied in many previous works, where isostaticity at the transition and scaling laws with the pressure of various quantities, including the contact number and the vibrational density of states, have been established. Additionally, much attention has been paid to nonspherical packings, and particularly recent work has made progress in understanding ellipsoidal packings. In the present work, we study the dimer packings in two dimensions, which have been much less understood than systems of spheres and ellipsoids. We first study the contact number of dimers near the jamming transition. It turns out that packings of dimers have "rotational rattlers", each of which still has a free rotational motion. After correcting this effect, we show that dimers become isostatic at the jamming, and the excess contact number obeys the same critical law and finite size scaling law as those of spheres. We next study the vibrational properties of dimers near the transition. We find that the vibrational density of states of dimers exhibits two characteristic plateaus that are separated by a peak. The high-frequency plateau is dominated by the translational degree of freedom, while the low-frequency plateau is dominated by the rotational degree of freedom. We establish the critical scaling laws of the characteristic frequencies of the plateaus and the peak near the transition. In addition, we present detailed characterizations of the real space displacement fields of vibrational modes in the translational and rotational plateaus.

Figures

Equations

$$V=\sum _{i>j}\sum _{n_i{n}_j}\frac{ϵ}{2}\text{pqty}{1-\frac{r_{k_i{k}_j}}{{\sigma }_{ij}}}^{2}H\text{pqty}1-\frac{r_{k_i{k}_j}}{{\sigma }_{ij}},$$

$$\varphi =\frac{N_s{b}_s^2+N_l{b}_l^2}{2L^2}\text{bqty}\pi -{\cos }^{-1}\text{pqty}\alpha -1+\text{pqty}\alpha -1\sqrt{\alpha \text{pqty}2-\alpha },$$

$$\mathcal{L}=\sum _{i=1}^N\text{pqty}\frac{1}{2}{{\dot{x}}_i}^2+\frac{1}{2}{{\dot{y}}_i}^2+\frac{I_i}{2}{{\dot{\theta }}_i}^2-V,$$

$$\text{derivative}[2]\mathbf{r}t=-\text{partialderivative}V\mathbf{r}.$$

$$\text{derivative}[2]\mathbf{u}t=-\mathcal{M}\mathbf{u},$$

$${\mathcal{M}}_{kl}=\text{partialderivative}V{u}_k{u}_l,$$

$$z_{\text{iso}}^N$$

$$z_B^N$$

$${\tilde{z}}_{\text{iso}}^N$$

$${\tilde{z}}_B^N$$

$$g\text{pqty}\omega =\frac{1}{d_{f}N-N^{\prime }-d}\sum _{k=1}^{d_{f}N-N^{\prime }-d}\delta \text{pqty}\omega -{\omega }^k,$$

$$T_k=\sum _{i=1}^N\text{bqty}\text{pqty}{e_{i,x}^k}^2+\text{pqty}{e_{i,y}^k}^2,$$

$$R_k=1-T_k.$$

$$p_k=\frac{\text{pqty}{{\sum }_{i=1}^N\text{absolutevalue}{{\mathbf{e}}_i^k}^2}^2}{N{\sum }_{i=1}^N\text{absolutevalue}{{\mathbf{e}}_i^k}^4}.$$

$$\text{derivative}[2]{\phi }_{l}t=-\sum _{k=1}^N{\mathcal{M}}_{{\phi }_l,{\phi }_k}\text{pqty}{\phi }_k-{\phi }_{0k}.$$

$${\tilde{e}}_{\text{pqty}i,\pm 1,x}^k$$

$${\tilde{e}}_{\text{pqty}i,\pm 1,y}^k$$

$${\tilde{p}}_k=\frac{\text{pqty}{{\sum }_{i=1}^{2N}\text{absolutevalue}{{\tilde{\mathbf{e}}}_i^k}^2}^2}{2N{\sum }_{i=1}^{2N}\text{absolutevalue}{{\tilde{\mathbf{e}}}_i^k}^4}.$$

$${\tilde{R}}_k=\frac{\text{absolutevalue}{\tilde{\mathbf{e}}}^k-{\tilde{\mathbf{e}}}_{\text{trans}}^k}{\text{absolutevalue}{\tilde{\mathbf{e}}}^k}.$$

$${\tilde{e}}_{\text{trans}\text{pqty}i,\pm 1,x}^k$$

$${\tilde{e}}_{\text{trans}\text{pqty}i,\pm 1,y}^k$$

$$I_i$$

$$={\rho }_i{\int }_{-a_i/2}^{a_i/2}\text{differential}x2x^{2}y\text{pqty}x+{\rho }_i{\int }_{-a_i/2}^{a_i/2}\text{differential}x\frac{2}{3}y\text{pqty}{x}^3,$$

$${\rho }_i=\frac{m}{\frac{{a_i}^2}{2{\alpha }^2}\text{bqty}\pi -{\cos }^{-1}\text{pqty}\alpha -1+\text{pqty}\alpha -1\sqrt{\alpha \text{pqty}2-\alpha }},$$

$$y\text{pqty}x=\{\begin{array}{ll}& \sqrt{\frac{b^2}{4}-\text{pqty}{x-\frac{a-b}{2}}^2},x\ge 0,\\ & \sqrt{\frac{b^2}{4}-\text{pqty}{x+\frac{a-b}{2}}^2},x<0.\end{array}$$

$$I_i={a_i}^{2}C\text{pqty}\alpha ,$$

$$C\text{pqty}\alpha =\frac{2{\alpha }^2}{\pi -{\cos }^{-1}\text{pqty}\alpha -1+\text{pqty}\alpha -1\sqrt{\alpha \text{pqty}2-\alpha }}\{{\int }_{-1/2}^0\text{differential}x\sqrt{\frac{1}{4{\alpha }^2}-\text{pqty}{x+\frac{\alpha -1}{2\alpha }}^2}\text{bqty}2x^2+\frac{2}{3}\text{pqty}\frac{1}{4{\alpha }^2}-\text{pqty}{x+\frac{\alpha -1}{2\alpha }}^2$$

$$+{\int }_0^{1/2}\text{differential}x\sqrt{\frac{1}{4{\alpha }^2}-\text{pqty}{x-\frac{\alpha -1}{2\alpha }}^2}\text{bqty}2x^2+\frac{2}{3}\text{pqty}\frac{1}{4{\alpha }^2}-\text{pqty}{x-\frac{\alpha -1}{2\alpha }}^2\}.$$

$$r_{k_i{k}_j}$$

$$X_{k_i{k}_j}$$

$$Y_{k_i{k}_j}$$

$$f_{k_i{k}_j}$$

$${\kappa }_{k_i{k}_j}$$

$${\mathcal{M}}_{x_a{x}_a}$$

$${\mathcal{M}}_{y_a{y}_a}$$

Full text

I Introduction

Jammed particulate systems composed of (frictionless) particles are known to undergo the jamming transition when they are compressed or decompressed van Hecke (2009). In particular, the system of spherical particles has been studied extensively by many previous works (see, e.g., Ref. O’Hern et al. (2003)). At the transition point where the pressure reaches zero, the system becomes isostatic. The isostaticity here means that the number of constraints of the system equals the number of degrees of freedom (DOF); that is, for sphere packings, the contact number equals $z_{\text{iso}}=2d$ where $d$ is the spatial dimension. Close to the transition, the excess contact number, $\Delta z\equiv z-z_{\text{iso}}$, scales as the square root of pressure $p$ or excess density $\Delta\phi\equiv\phi-\phi_{J}(>0)$: $\Delta z\propto p^{1/2}\propto\Delta\phi^{1/2}$. Mechanical properties such as elastic moduli also show the power-law scaling as a function of $p$ or $\Delta\phi$. Moreover, soft modes emerge associated with the jamming transition. The vibrational density of states (vDOS) shows the characteristic plateau, which extends towards the zero frequency at the transition point Silbert et al. (2005). The onset frequency of the plateau, which is denoted by $\omega^{*}$, scales linearly with the excess contact number: $\omega^{*}\propto\Delta z$. The variational argument Wyart et al. (2005a, b); Yan et al. (2016) and effective medium theory Wyart (2010); DeGiuli et al. (2014) have been proposed to explain this scaling behavior of vDOS.

Jammed systems composed of various shapes of nonspherical particles, such as ellipsoidal particles Donev et al. (2004a, b, 2007); Mailman et al. (2009); Zeravcic et al. (2009); Schreck et al. (2010); Delaney and Cleary (2010); Xia et al. (2014); Xu et al. (2015); Schaller et al. (2015); Yuan et al. (2018), sphero-cylinders Williams and Philipse (2003); Wouterse et al. (2009); Zhao et al. (2012); Meng et al. (2016); Marschall and Teitel (2018), and composites of spherical particles Baule et al. (2013); Gaines et al. (2016, 2017); Griffith and Hoy (2018), have also been studied numerically, theoretically, and experimentally. One of the most studied nonspherical particles is the ellipsoidal particle. It is found for ellipsoidal packings that the contact number becomes smaller than the isostatic value $z_{\text{iso}}=2d_{f}$ ($d_{f}$ is the number of DOFs that an ellipsoidal particle possesses) at the jamming transition, which is referred to as hypostaticity Donev et al. (2004a, 2007); Mailman et al. (2009); Zeravcic et al. (2009); Schreck et al. (2010, 2012). This hypostaticity has attracted attentions of previous works because it drastically alters the mechanical and vibrational properties near the jamming transition Mailman et al. (2009); Schreck et al. (2012); Zeravcic et al. (2009); Yunker et al. (2011); Brito et al. (2018). One might expect that the packings of ellipsoidal particles are not mechanically stable since the hypostaticity implies the shortage of constraints. However, they actually have positive elastic moduli and do not have zero-frequency modes in the jammed phase $\phi>\phi_{J}$. The hypostaticity indeed induces the floppy modes, which, however, are not exactly zero-frequency modes but are stabilized to become finite frequency modes by the repulsive interactions between ellipsoidal particles Mailman et al. (2009); Schreck et al. (2012). Recent theoretical work Brito et al. (2018) proposed a variational argument to describe the vibrational states in ellipsoidal packings and explain the scaling behavior of the vDOS, including that of these floppy modes, near the jamming transition. Here, we remark that ellipsoidal shape is not the only nonspherical shape to generate hypostatic packings. Sphero-cylinders also show hypostaticity in their contact numbers, and vibrational properties were reported to be similar to those of ellipsoidal packings VanderWerf et al. (2018); Marschall and Teitel (2018).

In this work, we focus on another kind of nonspherical particles, dimer particles. Compared to the ellipsoidal particles described above, dimer particles have been still much less understood. The ellipsoidal particles show two major differences from the spherical particles. One is asphericity of constituent particles, and the other is hypostaticity. As a result, the vDOS of the ellipsoids has a number of different characteristics compared to spheres Mailman et al. (2009); Schreck et al. (2012); Brito et al. (2018). However, the dimer particles are known to be isostatic at the jamming transition Schreck et al. (2010). Namely, the dimer particles possess asphericity, but they are not hypostatic in the jammed phase. In recent experimental and numerical works Han and Kim (2012); Han et al. (2013, 2015), structures of two-dimensional asymmetric dimer packings were studied. For vibrational properties, it has been reported that the dimer particles exhibit rotational vibrations in the low-frequency region below the translational band Schreck et al. (2010). However, the vDOS of the dimers has not been studied yet. As a result, it remains to be addressed whether the vDOS of the dimers exhibits the critical scaling behaviors near the jamming transition, such as the extension of the characteristic plateau, which is indeed our main topic in this work.

This article is organized as follows. In Section II, we describe the numerical protocol to generate dimer packings and also formulate the equation of motion and the dynamical matrix to study vibrational eigenmodes. Section III contains our numerical results for the dimer packings. We start from the numerical results on the jamming density. We next clarify that the present model has a limitation when the density is much higher than the jamming density. Then, we study the contact numbers of dimer packings near the jamming density. It turns out that when we analyze the contact numbers, we have to take care of the “rotational rattlers”, which are a new type of rattler. We finally present the results of vibrational properties of dimer systems, such as the vDOS, the contribution of rotational vibrations, and degree of localization. Importantly, we find two characteristic plateaus in the vDOS: the translational plateau at higher frequency and the rotational plateau at lower frequency. Accordingly, we define two characteristic frequencies: $\omega^{*}$ for the onset of the rotational plateau and $\omega_{R}$ for the rotational peak that separates the translational and rotational plateaus. We numerically establish the scaling behaviors of $\omega^{*}$ and $\omega_{R}$ near the jamming transition with the excess density $\Delta\phi$ and with the aspect ratio of dimer particles. Section IV summarizes our results and gives conclusions and remarks.

II Model and Methods

II.1 System description

We generated a numerical system composed of dimer particles in two dimensions (2D). The number of constituent dimers $N$ ranges from $60$ to $3600$ in the present work. Our model of a dimer particle is constructed by connecting two monomers that are modeled by harmonic disks. Figure 1 illustrates conformation of the dimer particle. We set the major axis of the dimer as $a$ and the minor axis as $b$, and define aspect ratio $\alpha$ as $\alpha=a/b$. The system is a binary mixture in order to avoid crystallization: $2N/3$ of dimer particles have minor axis $b=\sigma$ (major axis $a=\sigma\alpha$), and $N/3$ have $b=1.4\sigma$ ($a=1.4\sigma\alpha$), as in the previous work Schreck et al. (2010), where the value of $\alpha$ is fixed to be identical for all dimer particles. We also set the same mass for all dimer particles as $m$.

In our dimer model, each monomer of a dimer interacts with monomers belonging to the other dimers via purely repulsive harmonic potential Schreck et al. (2010). The total potential of the system is described as the summation of harmonic potential over all the monomers:

[TABLE]

where $i$ and $j$ represent the dimer $i$ and dimer $j$, respectively. $k_{i}$, which represents a monomer of the dimer $i$, is composed of a set of integers: $k_{i}\equiv\pqty{i,n_{i}}$, where $n_{i}=\pm 1$ designates the monomer that belongs to the dimer $i$. $\epsilon$ sets the energy scale, $r_{k_{i}k_{j}}$ is the distance between the monomer $k_{i}$ and monomer $k_{j}$, $\sigma_{ij}=\pqty{b_{i}+b_{j}}/2$ ($b_{i}$ and $b_{j}$ are the minor axes of dimers $i$ and $j$, i.e., the diameters of monomers $k_{i}$ and $k_{j}$, respectively), and $H\pqty{x}$ is the Heaviside step function: $H\pqty{x}=1$ for $x\geq 0$ and $H\pqty{x}=0$ for $x<0$. Throughout this paper, we use $m$, $\sigma$, and $\epsilon$ as units of mass, length, and energy, respectively. We note that frequency is measured by $\tau=\pqty{m\sigma^{2}/\epsilon}^{-1/2}$.

We first prepared a random configuration at sufficiently small packing fraction $\phi_{0}=$0.2$$ and minimized the potential energy using the FIRE algorithm Bitzek et al. (2006). Here, the periodic boundary conditions are implemented on the system. The packing fraction $\phi$ of 2D dimers is defined as Schreck et al. (2010)

[TABLE]

where $L$ is the system length, $N_{s}=2N/3$ ($N_{l}=N/3$) is the number of smaller (larger) particles, and $b_{s}=\sigma$ ($b_{l}=1.4\sigma$) is the length of the minor axis of smaller (larger) particles. We then gradually compressed the system by increments of the packing fraction of $\delta\phi=${10}^{-3}. After each compression, the energy minimization was executed. Here, we note that the minimization is terminated when either of the following two conditions is satisfied Gao *et al.* ([2006](#bib.bib37)): (i) The total potential per particles is sufficiently small $V/N<${10}^{-16}; or (ii) The total potential energies of two successive FIRE steps $t$ and $t+1$ are nearly the same: $\absolutevalue{V_{t+1}-V_{t}}/V_{t}<${10}^{-16}$. We finally obtained the jamming transition point as the point where the potential energy per particle is in a range of$1\text{\times}{10}^{-16}$<V/N<$2\text{\times}{10}^{-16}$$. To numerically achieve this aim, we repeated compression and decompression on the system around the transition point. Each time the compression/decompression is switched, $\delta\phi$ is multiplied by $0.5$.

Once we obtained the configuration at the transition point, we compressed the system to obtain the configurations as a function of the excess density $\Delta\phi=\phi-\phi_{J}$. In these obtained packings, we recursively removed “rattler” particles that have contacts fewer than the number of DOF, $d_{f}\equiv 3$ (two translational DOFs and one rotational DOF). In the following, the number of particles without rattler particles is denoted as $N$.

II.2 Vibrational mode analysis

The Lagrangian $\mathcal{L}$ of the system is

[TABLE]

where $x_{i}$ and $y_{i}$ are spatial positions of the center of mass of the dimer $i$, and $\theta_{i}$ is rotational angle of the dimer $i$ with respect to the $x$-axis. $V$ is the total potential energy defined by Eq. (1), and $I_{i}$ is the inertia of the dimer $i$. The detailed formulation of inertia is given in Appendix A. We employ the normal coordinate of dimer $i$ as $\pqty{x_{i},y_{i},\varphi_{i}}$ with $\varphi_{i}=\sqrt{I_{i}}\theta_{i}$ and that of whole the system as $\mathbf{r}=\pqty{x_{1},y_{1},\varphi_{1},\dots,x_{N},y_{N},\varphi_{N}}$. Then, the equation of motion (Euler-Lagrange equation) is given as

[TABLE]

When we suppose the harmonic vibrations around the equilibrium coordination $\mathbf{r}_{0}$, the equation of motion (4) is described as

[TABLE]

where $\mathbf{u}\equiv\mathbf{r}-\mathbf{r}_{0}$ is the displacement from $\mathbf{r}_{0}$. $\mathcal{M}$ is the so-called dynamical matrix, and its elements are

[TABLE]

where $k,l=1,\dots,d_{f}N$, and $u_{k}$ is the $k$th component of $\mathbf{u}$. Explicit formulations of $\mathcal{M}_{kl}$ are given in Appendix B. By diagonalizing the dynamical matrix, we obtained a set of eigenvalues, $\lambda^{k}$, and eigenvectors, $\mathbf{e}^{k}\equiv\pmqty{undef}$, with $\mathbf{e}^{k}_{i}=\pmqty{undef}$ ($i=1,\dots,N$). The eigenfrequencies are given as $\omega^{k}=\sqrt{\lambda^{k}}$, and the eigenvectors are orthonormalized as $\mathbf{e}^{k}\cdot\mathbf{e}^{l}\equiv\sum_{i=1}^{N}\mathbf{e}^{k}_{i}\cdot\mathbf{e}^{l}_{i}=\delta_{kl}$ where $\delta_{kl}$ is the Kronecker delta.

In the present work, we performed vibrational eigenmode analysis on the “unstressed system”. Since the monomers are modeled by harmonic disks in the present model, forces acting between the monomers are always repulsive. For this reason, we refer to this original state as the stressed system. In the unstressed system, we retain all the stiffness (i.e., the second derivative of potential) between monomers but drop the force in the analysis. For the case of sphere packings, the theoretical understanding was first constructed based on the unstressed system Wyart et al. (2005a); Wyart (2010), which was then extended to the stressed system by considering effects of the forces Wyart et al. (2005b); DeGiuli et al. (2014). The unstressed system of sphere packings exhibits the characteristic plateau in vDOS and the critical behavior of this plateau near the jamming transition. The forces in sphere packings do not affect this plateau and its critical behavior, whereas they make the system mechanically unstable Lerner et al. (2014) and alter the nature of the very low-frequency vibrational modes Mizuno et al. (2017). In particular, quasilocalized vibrational modes are induced by the repulsive forces Mizuno et al. (2017); Lerner and Bouchbinder (2018); Shimada et al. (2018). In the present work, we limit ourselves to analyzing the vibrational properties of unstressed dimer packings. Based on the results of unstressed systems, we will discuss the role of repulsive forces in dimer packings in the near future.

III Results

III.1 Jamming density

By following the procedure described in Section II.1, we generated the configurations at the jamming transition where the potential per particle $V/N$ is in 1\text{\times}{10}^{-16}$<V/N<$2\text{\times}{10}^{-16}. Figure 2 plots the jamming density $\phi_{J}$ as a function of the aspect ratio $\alpha$. When $\alpha\to 1$, the jamming density approaches that for the binary mixtures of disks, $\phi_{J}\approx$0.84 O’Hern *et al.* ([2003](#bib.bib2)). We also observe the nonmonotonic dependence of $\phi_{J}$ on $\alpha$, where the peak value is taken at approximately $\alpha=$1.4. These results are consistent with previous observation by Ref. Schreck et al. (2010). Starting from the configurations at the jamming transition, we compress the system to obtain the configurations in the jammed phase $\phi>\phi_{J}$ (the excess density $\Delta\phi>0$). We stress that we obtained $\phi_{J}$ for each configuration; then, $\Delta\phi$ was defined for each configuration separately.

III.2 Limitations of the present model of dimers

The rigorous way to simulate dimer particles is to consider the particle as one elastic element of the dimeric shape. In this case, the potential energy of the system is given as the summation of the elastic deformation energies of all the dimer particles. By contrast, the present model is an approximate model, where the potential energy is evaluated as the summation of the interparticle potentials between constituent monomers. This approach introduces some limitations of the present model, which we now clarify.

There are two types of contacts that cannot be properly addressed by the present model. One is “double” contact, which is illustrated in Fig. 3 (a). When a monomer of a dimer deeply overlaps with another dimer, the monomer makes contact with both constituent monomers of the other dimer. This is an artificial effect in the present model, and as such, the contact cannot be properly addressed. The other is the “cusp” contact, as illustrated in Fig. 3 (b). In this contact, a monomer of a dimer contacts another dimer covering whole the cusp line of the other dimer. The present model is reasonable only for contacts where the deformation of a dimer can be decomposed into the deformations of two monomers. This “cusp” contact clearly violates this condition, and such contacts cannot be properly addressed. Thus, as long as we employ the present model, we have to discard configurations that have any of the double and the cusp types of contacts. These contacts can occur at the higher packing fraction, and therefore, this condition gives an upper limit of the excess volume fraction $\Delta\phi$ that can be studied in the present model.

We count the number of double and cusp contacts in configurations at each excess density $\Delta\phi$. Figure 4 plots the number of these contacts as a function of $\Delta\phi$ for the case of aspect ratio $\alpha=$1.1. The system size is $N=$3600. The number of cusp contacts is exactly zero up to $\Delta\phi=${10}^{-3} and then becomes finite for $\Delta\phi\geq$e-2.5. The number of the double contacts is exactly zero up to $\Delta\phi=$e-1.5. Thus, configurations with $\Delta\phi\leq${10}^{-3} include neither of these two contacts, and the upper limit of the excess density is determined as $\Delta\phi=${10}^{-3} for $\alpha=$1.1. We perform the same analysis for several different aspect ratios and obtain the upper limits of the excess densities, as $\Delta\phi=$e-1.5 for $\alpha=$1.5, and $\Delta\phi=$e-4.5 for $\alpha=$1.03. We find that these values of the upper limits are independent of the system size. In the present work, we limit ourselves to the densities below these upper limits.

III.3 Contact number

We next study the contact number $z_{J}$ at the jamming transition $\phi_{J}$. We first review the known results for $z_{J}$ of spherical particles. Generally, when each particle has $d_{f}$ DOFs, the isostatic condition in the infinite system size is $z_{\text{iso}}=2d_{f}$. For spheres, the number of DOFs equals the spatial dimension, $d_{f}=d$; then, $z_{\text{iso}}=2d$. However for finite-size systems with $N$ particles, two types of finite-size effects appear Goodrich et al. (2012). First, because the systems are allowed to have $d$ global translations, the isostatic contact number is reduced to $z^{N}_{\text{iso}}=z_{\text{iso}}-2d/N$. Second, because the systems must have positive bulk modulus, the isostatic contact number is increased to $z^{N}_{B}=z^{N}_{\text{iso}}+2/N=2d-2d/N+2/N$. It is established for spheres that $z_{J}=z^{N}_{B}$ holds for any system size Goodrich et al. (2012). Note that both $z^{N}_{\text{iso}}$ and $z^{N}_{B}$ converge to $z_{\text{iso}}$ in the thermodynamic limit of $N\to\infty$.

We now discuss whether this argument holds for the dimer packings in 2D. Following the above arguments, we obtain the isostatic contact numbers for dimers:

[TABLE]

with the DOFs $d_{f}=3$ and the spatial dimension $d=2$ in our 2D dimer system. To test the relation of $z_{J}=z^{N}_{B}$ for dimers, we numerically measured $z_{J}$ of dimers for various aspect ratios and plot the difference $z_{J}$ from $z_{B}^{N}$ of Eq. (8) in Fig. 5. Clearly, $z_{J}$ takes a smaller value than $z^{N}_{B}$, and the difference grows larger as the aspect ratio increases.

To understand the origin of $z_{J}<z^{N}_{B}$, we look into the vibrational eigenmodes obtained from the dynamical matrix 111Here, we study vibrational eigenmodes in the original stressed system. On the other hand, all of the results presented in Section III.4 are obtained from the unstressed system. Note that the stressed and unstressed system are identical exactly at the jamming density.. Remarkably, we find zero-frequency modes, in addition to two trivial zero-frequency modes of global translations. Figure 6 compares the number of these nontrivial zero-frequency modes (squres) and the number of shortages of contacts, $N\pqty{z^{N}_{B}-z_{J}}/2$ (circles). These two numbers are the same for all the studied aspect ratios. Thus, the zero-frequency modes directly correspond to the shortage of $z_{J}$ from $z^{N}_{B}$.

Then, we look into the nature of these zero-frequency modes. We find that these modes are spatially localized to the dimers that have characteristic contacts as illustrated in Fig. 7: they have more than three contacts, i.e., they are not rattler particles, but all of the contacts are concentrated in one monomer while the other monomer has no contacts. These dimers can obviously rotate around one monomer without any energy cost; thus, they can produce zero-frequency modes Alexander (1998). These dimers can be considered as another type of rattler particle: their translational DOFs are constrained by their three or more contacts, but the rotational DOFs are not. In this paper, we refer to these dimer particles to as “rotational rattlers”. Figure 6 (triangles) shows that the number of rotational rattlers coincides with the number of nontrivial zero-frequency modes and the shortage of contacts for all the studied aspect ratios. Thus, we conclude that the dimer packings have rotational rattlers that cause the nontrivial zero-frequency vibrational modes in addition to the apparent shortage of contacts 222The number of rotational rattlers, the number of nontrivial zero-frequency modes, and the shortage of contacts are exactly the same for approximately $80\text{\,}\mathrm{\char 37\relax}90\text{\,}\mathrm{\char 37\relax}$ of the configurations at $\phi_{J}$. They are almost the same but not exactly the same for the rest $10\text{\,}\mathrm{\char 37\relax}20\text{\,}\mathrm{\char 37\relax}$ of the configurations, for which we do not have clear explanation..

For the proper definition of the isostaticity of dimers, the rotational DOFs of rotational rattlers should be excluded, as in the case of the translational DOFs of usual rattlers. This exclusion is achieved by replacing the total number of DOFs $Nd_{f}$ with $Nd_{f}-N^{\prime}$, where $N^{\prime}$ is the number of rotational rattlers. This leads to the modified definition of the isostatic contact numbers for dimers:

[TABLE]

Hereafter, we use these definitions as the isostatic contact numbers.

Finally, we focus on the contact number at $\phi>\phi_{J}$. Although we focus on $\alpha=$1.5$$ in the following, all the other aspect ratios exhibit similar behaviors. In the case of spherical particles, the excess contact number $\Delta z=z-z_{\text{iso}}$ follows the square root law $\Delta z\sim\Delta\phi^{1/2}$ O’Hern et al. (2003). More precisely, the excess contact number of $N$ spherical particles follows the finite-size scaling law $N\pqty{z-z^{N}_{\text{iso}}}=F\pqty{\Delta\phi^{1/2}N}$, where $F(x)\sim 1$ for small $x$ and $F(x)\sim x$ for large $x$ Goodrich et al. (2012) 333In Goodrich et al. (2012), the pressure $p$ is used instead of the excess packing fraction $\Delta\phi$.. The existence of the finite-size scaling law demonstrates that the jamming transition is a phase transition. Here, we test whether these laws hold for the jamming transition of dimers or not. Employing the value of $\tilde{z}^{N}_{\text{iso}}$ in Eq. (9), we plot the excess contact number of dimers $z-\tilde{z}^{N}_{\text{iso}}$ as a function of the excess density $\Delta\phi$ for several different system sizes $N$ in Fig. 8 (a). At the high packing fractions, the excess contact number follows the square root law $\Delta z\sim\Delta\phi^{1/2}$, as in the case of spheres O’Hern et al. (2003). We also observe that $\Delta z$ deviates from the square root law when approaching the jamming density, which is a signature of the finite-size effect. We then replot data in the form of $N\pqty{z-\tilde{z}^{N}_{\text{iso}}}$ versus $\Delta\phi^{1/2}N$ in Fig. 8 (b). Clearly, all the data points collapse onto a single master function $N\pqty{z-\tilde{z}^{N}_{\text{iso}}}=G\pqty{\Delta\phi^{1/2}N}$, where $G\pqty{x}\sim 1$ for small $x$ and $G\pqty{x}\sim x$ for large $x$.

Therefore, the results in this subsection establish that, when the rotational rattlers are properly excluded, the excess contact number of dimers near the jamming obeys the same critical law and finite-size scaling law as those of spheres. These results also allow us to determine the range of densities in which the finite-size effect is negligible, that is, $\Delta\phi\geq${10}^{-7} for $N=$3600 and $\Delta\phi\geq${10}^{-6} for $N=$1200. Note that this lower bound can be extended by increasing the system size, while the upper bound discussed in Section III.2 cannot.

III.4 Vibrational properties

In this subsection, we will discuss the nature of vibrational eigenstates, which are obtained by diagonalizing the dynamical matrices as explained in Section II.2. We emphasize that the results are from the unstressed system.

To characterize the statistics of the eigenfrequencies $\omega^{k}$, we calculate the vDOS:

[TABLE]

where $\delta\pqty{x}$ is the Dirac delta function. We also analyze the nature of the eigenvectors $\mathbf{e}^{k}$ by quantifying the contribution of translational and rotational DOFs, and the participation ratio. The contribution of translational DOFs to mode $k$ is defined as

[TABLE]

as in the case of ellipsoids Mailman et al. (2009). The contribution of rotational DOFs is the counterpart of $T_{k}$, and defined as

[TABLE]

The participation ratio $p_{k}$ is defined as

[TABLE]

as in the case of spheres Mazzacurati et al. (1996). The participation ratio $p^{k}$ measures the extent of vibrational localization in the mode $k$. We emphasize that $R_{k}$ and $p_{k}$ give the information of the eigenvectors in normal coordinates described in Section II.2.

III.4.1 General nature of vibrational eigenmodes

We first focus on the vibrational states of dimers with the aspect ratio $\alpha=$1.5, which are presented in Fig. [9](#S3.F9). The excess packing fraction is $\Delta\phi=${10}^{-2}, and the system size is $N=$3600. The top panel of Fig. [9](#S3.F9) shows the vDOS. The structure of the vDOS is very similar to that of sphere packings. In particular, it has the characteristic plateau at intermediate frequencies. The second top panel shows the contribution of rotational DOFs, $R_{k}$. We see that $R_{k}\approx$0 for the high-frequency edge and low-frequency regime below the plateau. On the other hand, $R_{k}$ takes approximately $0.3$ to $0.6$ in the intermediate, plateau regime. Therefore, the rotational DOFs are relevant only in the vibrational states in the plateau region. The third top panel shows the participation ratio, $p_{k}$. The plot demonstrates that the modes in the high-frequency edge are strongly localized and the modes in the plateau regime are extended, which are very similar to the case for spheres. We also observe that the modes in the low-frequency regime below the plateau are largely extended and their eigenfrequencies are distributed discretely, which suggests that these modes are phonons. We indeed confirm that the frequencies of these modes are consistent with the discretized levels of acoustic phonons of this system. Therefore, at $\alpha=$1.5$$, the vibrational states of dimers are composed of the localized high-frequency edge, the intermediate plateau, and low-frequency acoustic phonons. These features are very similar to those of spheres, apart from the fact that the rotational DOFs are relevant in the plateau modes for dimers.

We turn our attention to the vibrational states of dimers with the aspect ratio $\alpha=$1.1. The excess packing fraction is $\Delta\phi=${10}^{-3}, and the system size is $N=$3600. The vDOS, $R_{k}$, and $p_{k}$ are presented in Fig. [10](#S3.F10), in the same way as in Fig. [9](#S3.F9). The vDOS for $\alpha=$1.1 appears to be very different from that for $\alpha=$1.5; we observe two characteristic plateaus, and they are separated by a peak. $R_{k}\approx$0 in the plateau of high frequency, while $R_{k}\approx$0.8 in that of low frequency. This result indicates that the rotational DOFs do not play any role in the vibrations in the high-frequency plateau, while they are relevant in the low-frequency plateau. Hereafter, we refer to the high-frequency plateau as the “translational plateau” and the low-frequency plateau as the “rotational plateau”. The participation ratio of the translational plateau modes is $p_{k}\approx$0.3, while the rotational plateau modes take smaller values of $p_{k}\lesssim$0.1. This observation suggests that the rotational plateau modes can be localized in the space of normal coordinates (we will discuss this point based on the real space displacement field later in Section [III.4.3](#S3.SS4.SSS3)). Interestingly, $R_{k}$ takes the largest values in the peak at approximately $\omega\approx$0.2. We refer to this peak as the “rotational peak”. In the regime of 0.06$\lesssim\omega\lesssim$0.3 including this rotational peak, the distribution of eigenfrequencies are discretized. We speculate that these discretized modes are optical types of modes. To fully clarify the nature of these modes, however, we need to analyze the dynamic structure factors Wang et al. (2015), which is beyond the scope of this work. In summary, the vibrational states of dimers at $\alpha=$1.1$$ are characterized by the emergence of the translational and rotational plateaus that are separated by the rotational peak.

We repeated a similar analysis for various different aspect ratios (1.001$\leq\alpha\leq$1.9). We find that the results of larger aspect ratios are similar to that for $\alpha=$1.5, while those for smaller aspect ratios are similar to $\alpha=$1.1. Therefore, the vibrational states of dimers are categorized into these two types of behaviors, depending on the aspect ratio.

III.4.2 Critical behaviors near the jamming transition

We now study the packing fraction dependence of vDOS near the jamming transition. Figure 11 plots the vDOS $g\pqty{\omega}$ for several different values of packing fraction $\Delta\phi$ and for aspect ratio $\alpha=$1.5 in (a), $\alpha=$1.1 in (b), and $\alpha=$1.03 in (c). For (a) $\alpha=$1.5, the vDOS shows the characteristic plateau at intermediate frequencies. We recall that the vibrational states in the plateau have contributions from rotational DOFs. As we approach the jamming transition as $\Delta\phi\to 0$, this plateau extends towards zero frequency, which is similar to that observed in the vDOS of spheres. For (b) $\alpha=$1.1 and (c) $\alpha=$1.03, the vDOSs have the translational plateau at high frequency and the rotational plateau at low frequency, which are separated by the rotational peak. As we approach the transition, the rotational plateau extends towards zero frequency, while the translational plateau and the rotational peak are nearly unchanged.

From the above observations, we define two characteristic frequencies: $\omega_{R}$ for the frequency of the rotational peak, and $\omega^{*}$ for the onset of the rotational plateau. We first study the excess density $\Delta\phi$ and aspect ratio $\alpha$ dependencies of $\omega_{R}$. We used $N=$1200 systems for this analysis. Note that although the rotational peak is absent for $\alpha=$1.5, we still observe the sharp increase in $R_{k}$ as shown in Fig. 9. Therefore as a reference, we report the frequency at which $R_{k}$ takes the maximum as $\omega_{R}$ for $\alpha=$1.5$$. Figure 12 plots the $\omega_{R}$ against $\alpha-1$ for various $\Delta\phi$. $\omega_{R}$ is nearly independent from $\Delta\phi$ but increases with $\alpha$. This increase can be fitted into the linear function of $\alpha-1$, especially for the small $\alpha$. This $\alpha$ dependence is natural for the fundamental frequency scale of rotation, as explained below. When we consider purely rotational vibration, the equation of motion (Eq. (5)) reduces to

[TABLE]

As shown in Eq. (41) of Appendix B, the dynamical matrix elements $\mathcal{M}_{\varphi_{l},\varphi_{k}}$ are proportional to $\pqty{\alpha-1}^{2}$; therefore, the fundamental frequency scale of rotation is proportional to $\alpha-1$. This $\alpha$ dependence originates from the geometric fact that the displacement of monomers due to the pure rotation of dimers is proportional to $\alpha-1$ and causes an energy cost proportional to $\pqty{\alpha-1}^{2}$ (with constant stiffness). Therefore, $\omega_{R}\propto\pqty{\alpha-1}$ is reasonable geometrically.

In relation to the scaling behavior of $\omega_{R}\propto\pqty{\alpha-1}$, we can discuss the heights of the rotational and translational plateaus of the vDOS. Figures 11 (b) and (c) show that the rotational plateau increases with decreasing $\alpha$, whereas the translational plateau decreases with decreasing $\alpha$. These tendencies are reasonable as discussed below. We establish that the rotational peak at $\omega_{R}$ separates the rotational and translational plateaus. Suppose that most of the translational DOFs and the rotational DOFs are confined in the translational ($\omega>\omega_{R}$) and rotational ($\omega<\omega_{R}$) plateaus, respectively, which is reasonably supported by the values of $R_{k}$ in Fig. 10. Then, when $\omega_{R}$ becomes close to the zero frequency as $\alpha$ approaches unity, the frequency regime in which rotational DOFs emerge decreases, and the height of the rotational plateau increases. We can observe this behavior in Fig. 11.

We next focus on the excess density $\Delta\phi$ and aspect ratio $\alpha$ dependencies of $\omega^{*}$. Figure 13 plots $\omega^{*}$ as a function of the excess contact number $\Delta z$ for different values of aspect ratio $\alpha=$1.03$,\$1.1$,\$1.5$$ 444Practically, the value of $\omega^{*}$ is determined as the frequency where $g\pqty{\omega}/\omega^{1/4}$ takes a maximum.. Remarkably, $\omega^{*}$ scales linearly with $\Delta z$, as $\omega^{*}\propto\Delta z$. For spherical particles, the vDOS shows the characteristic plateau, and its onset frequency $\omega^{*}$ also scales linearly with $\Delta z$ Silbert et al. (2005). Therefore, our result establishes the common scaling law for the onset of the “rotational plateau” of dimers and the plateau of spheres. In addition, $\omega^{*}$ slightly decreases as the aspect ratio $\alpha$ goes to unity. However, this tendency is very subtle, and the precise determination of the $\alpha$ dependence of $\omega^{*}$ requires the data at $\alpha$ to be much closer to unity. This requires conducting simulations with a much larger system size, which are not available in this work.

III.4.3 Real space displacements in the plateau modes

In the previous subsections, we established that the vDOSs of the dimers with a small aspect ratio have the translational plateau and rotational plateau, which exhibit the critical behaviors near the jamming transition. These two plateaus can be distinguished by the rotational DOFs, i.e., they are relevant in the rotational plateau but not in the translational plateau. This result is true in the space of normal coordinates. This subsection focuses on the real space displacement fields in the translational and rotational plateau modes.

We analyze the displacement fields of the vibrational modes in real space. In normal coordinates, the vibrational eigenstates are expressed by the translational and rotational components, $\mathbf{e}^{k}_{i}=\pmqty{undef}$, for dimer $i$. However, we can also express the vibrational eigenstates as displacement fields of monomers, which we denote by $\tilde{\mathbf{e}}^{k}_{i}=\pmqty{undef}$ for dimer $i$. These two descriptions can be transformed from each other by

[TABLE]

We note that $\mathbf{e}^{k}=\pmqty{undef}$ describes the vibrations in the phase space of $(x,y,\varphi)$, where translational and rotational motions are treated as different objects. On the other hand, $\tilde{\mathbf{e}}^{k}=\pmqty{undef}$ describes the vibrations in the real space of $(x,y)$, where both translational and rotational motions are expressed as translational motions of monomers, namely, treated on equal footing. We can define the participation ratio for the real space displacement fields as

[TABLE]

We can also define the contribution of the rotational DOFs for the real space displacement fields as

[TABLE]

Here, $\tilde{\mathbf{e}}^{k}_{\text{trans}}$ is defined as

[TABLE]

which corresponds to the displacement fields of monomers that originated from only the translational DOFs in $\mathbf{e}^{k}$.

In Fig. 14, we visualize the displacement fields, $\tilde{\mathbf{e}}^{k}$, in four representative eigenstates for the aspect ratio $\alpha=$1.1$$, whose frequencies are indicated in the top panel of Fig. 10. The panel (A) shows a mode in the high-frequency edge. This vibrational state is highly localized onto a few particles. The panel (D) shows a vibrational state in the low frequency edge below the rotational plateau. This vibrational state is more ordered and spatially extended, which exhibits a feature of an acoustic phonon. These features are similar to those of the corresponding vibrational modes in sphere packings Silbert et al. (2009).

Now, we discuss the features of translational and rotational plateau modes. The panel (B) shows a translational plateau mode, while the panel (C) shows a rotational plateau mode. In both cases, the vibrational states are disordered and spatially extended. These features are very similar to those of the plateau modes in sphere packings Silbert et al. (2009). We note that the rotational plateau mode (C) appears to be more “swirly” compared to the translational plateau mode (B). However, this difference cannot be attributed to the rotational nature of the modes because sphere packings already show similar swirly motions in the vibrational modes in the lower-frequency part of the plateau van Hecke (2009); Silbert et al. (2009). Apart from the swirly motions, the visualization of real space displacements does not indicate any strong difference between these two types of modes. In contrast, the values of $R_{k}$ and $p_{k}$ indicate large differences between these two modes, as shown in Fig. 10. We discuss this point below.

First, the displacement field of the rotational plateau mode (C) does not appear to be very “rotational”, regardless of high values of $R_{k}\approx$0.8. In the displacement fields, the purely rotational motion of dimers appears as the pair of arrows (rooted at two monomers) with the opposite directions. However in the panel (C), such pairs of arrows cannot be found. This visual inspection is supported by the values of $\tilde{R}_{k}$ shown in Fig. [15](#S3.F15): $\tilde{R}_{k}$ indicates much lower values of $\tilde{R}_{k}\approx$0.2 than those of $R_{k}\approx$0.8$$ for the rotational plateau modes. The large difference between $R_{k}$ and $\tilde{R}_{k}$ can be understood when we go back to the transformation formula Eqs. (16) and (17); the factor of $\pqty{\alpha-1}$ suppresses the impact of the rotation in the displacements in the real space. This situation is analogous to the rotation of a rod. When the rod is very short, a large rotation does not cause large displacements of its edges. Therefore, these results show that the translational motions of dimers are more dominant than the rotational motions of dimers for the real space displacements, even for the rotational plateau modes.

In addition, the displacement field of the rotational plateau mode (C) does not appear to be very localized, regardless of low values of $p_{k}\lesssim$0.1. This visual inspection is again supported by the values of $\tilde{p}_{k}$ shown in Fig. [15](#S3.F15): $\tilde{p}_{k}$ indicates the higher values of $\tilde{p}_{k}\approx$0.3 than those of $p_{k}\lesssim$0.1$$ for the rotational plateau modes, which are almost the same as those of the translational plateau modes. This observation suggests that the vibrational localization takes place only in the rotational DOFs that do not strongly affect the real space displacements.

In summary, these analyses reveal that the real space displacements for both the translational and rotational plateau modes are disordered and spatially extended. Remarkably, the displacement field of the rotational plateau mode is still dominated by the translational motions of dimers, and the rotational motions only introduce slight modifications to it, as evidenced by the low values of $\tilde{R}_{k}\approx$0.2. We emphasize, however, that the distinction between the translational and rotational plateau modes is still clear and definite. For translational plateau modes, the rotational motions are fully frozen and never introduce even small modifications to the displacement field, as evidenced from nearly zero values of $\tilde{R}_{k}\approx$0.

IV Conclusion

In this work, we study critical behaviors of vibrational properties near the jamming transition for dimer packings in two dimensions. Our dimer particle is modeled by connecting two harmonic disks. However, we revealed that this usual modeling of dimers encounters problems that involve the unusual contacts between dimers, i.e., the double contact and the cusp contact. We then clarified that the unusual contacts appear only at high packing fractions and that the model is free from this problem near the jamming transition. Our main results are obtained from the density region without this problem.

We first studied the contact number of dimers near the jamming transition. We found that the contact number at the jamming transition, $z_{J}$, is smaller than the naive isostatic contact number of dimers. This apparent shortage of contact number originated from the presence of “rotational rattlers”, in which the dimer has three contacts or more, but these contacts are concentrated in one monomer of the dimer so that the dimer still has a free rotational motion. Therefore, we corrected this effect of rotational rattlers on the isostatic contact number and defined the excess contact number $\Delta z$ as that measured from this corrected isostatic contact number. We then established that $\Delta z$ defined in this way obeys the same critical law and finite size scaling law as those of sphere packings; namely, it follows the square root law of $\Delta z\propto\Delta\phi^{1/2}$ for large systems, as well as the finite size scaling law of $N\Delta z=G\pqty{\Delta\phi^{1/2}N}$, where $G\pqty{x}\sim 1$ for small $x$ and $G\pqty{x}\sim x$ for large $x$.

We next studied the vibrational properties of dimers near the jamming transition. The vDOS of dimers with a large aspect ratio is composed of a high-frequency edge, the plateau in the intermediate frequency, and the acoustic phonon in the low-frequency edge, which are very similar to the features of vDOS of spheres. We note, however, that the vibrational states in the intermediate plateau have a rotational nature. Interestingly, for the case of the small aspect ratio, the plateau at the intermediate frequency is separated into two plateaus by a peak. The vibrational states of the high- and low-frequency plateaus are dominated by the translational and rotational DOFs, respectively, which are thus called “translational plateau” and “rotational plateau”. We also call the peak that separates the two plateaus the rotational peak. However, we remark that the displacement fields in the real space are still dominated by the translational DOFs even for the modes in the rotational plateau.

In addition, we defined two frequencies to characterize the above features of vDOS: $\omega_{R}$ for the frequency of the rotational peak and $\omega^{*}$ for the onset of the rotational plateau. The frequency $\omega_{R}$ to characterize the rotational peak is independent of the packing fraction but depends on the aspect ratio as $\omega_{R}\propto\pqty{\alpha-1}$, which is reasonable as the fundamental frequency scale of rotational vibrations in this model. On the other hand, the onset of the rotational plateau $\omega^{*}$ depends on the packing fraction; it scales as $\omega^{*}\propto\Delta z$, which is the same scaling law as that of the plateau in spheres. These results establish a fundamental view on the vibrational properties of dimers near the jamming transition. In this work, we restricted ourselves to the vibrational properties of the unstressed system. Of course, these findings should be supplemented with studies of original, stressed systems. We are now working on this direction.

As a final remark, we would compare the vibrational properties of dimers with those of ellipsoids Schreck et al. (2012); Brito et al. (2018). The characteristic frequencies of the rotational bands scale linearly with $\alpha-1$ in ellipsoids, which is the same as the scaling behavior of $\omega_{R}$ in dimers. On the other hand, ellipsoids have the quartic vibrational modes in the lowest-frequency region, which are totally absent in dimers. Instead, dimers have the extended rotational plateau in the vDOS, which does not appear in the vDOS of ellipsoids. We tentatively attribute these differences in vibrational properties to the difference between “hypostaticity” of ellipsoids and “isostaticity” of dimers. However, for a deeper understanding of the rotational plateau modes in dimer packings, we need further numerical and theoretical studies.

Acknowledgements.

We thank Harukuni Ikeda and Daniele Coslovich for useful discussions and suggestions. This work was supported by JSPS KAKENHI Grant Numbers 17K14369, 17H04853, 16H04034, 18H05225, and 19K14670. This work was also partially supported by the Asahi Glass Foundation. The computations were partially performed using the Research Center for Computational Science, Okazaki, Japan.

Appendix A Moment of inertia of a dimer

In this appendix, we provide the formulation of the moment of inertia of a dimer $i$, $I_{i}$.

We set a coordinate system, $\pqty{x,y}$, on the dimer particle that is illustrated in Fig. 1. We set the origin of this coordinate system at the center of the dimer, and $x$ and $y$ axes along the major and the minor axes of the dimer, respectively. Then, the moment of inertia of the dimer $i$, $I_{i}$, is formulated as

[TABLE]

where $\rho_{i}$ is an areal density of the 2D dimer $i$ with:

[TABLE]

and $y\pqty{x}$ is

[TABLE]

By performing the integration in Eq. (23), we obtain the moment of inertia $I_{i}$ as

[TABLE]

where $C\pqty{\alpha}$ is a function of only the aspect ratio $\alpha$:

[TABLE]

The behavior of $C\pqty{\alpha}$ is shown in Fig. 16.

Appendix B Dynamical matrix

Here, we provide the explicit formulations of the elements of dynamical matrix, $\mathcal{M}_{kl}$, defined in Eq. (6) ($k,l=1,\dots,d_{f}N$).

The distance between the monomer $k_{i}$ of the dimer $i$ and the monomer $k_{j}$ of the dimer $j$ is formulated as

[TABLE]

where $n_{i}$ ($n_{j}$) takes values of $1$ or $-1$ to represent the monomer of the dimer $i$ ($j$). Additionally, the force $f_{k_{i}k_{j}}$ and the stiffness $\kappa_{k_{i}k_{j}}$ acting between these two monomers (which correspond to the first and second derivatives of the potential, Eq. (1), respectively) are

[TABLE]

where we recall that $H\pqty{x}$ is the Heaviside step function: $H\pqty{x}=1$ for $x\geq 0$ and $H\pqty{x}=0$ for $x<0$. Then, we can formulate the elements of dynamical matrix, $\mathcal{M}_{kl}$, as follows. Note that we use variables $\xi_{i}$ to designate elements of the dynamical matrix ($\xi=x,y,\varphi$) in the following. In the sense of the definition in Eq. (6), this should be converted as $k=d_{f}\pqty{i-1}+n_{\xi}$, where $n_{\xi}=1$ for $\xi=x$, $n_{\xi}=2$ for $\xi=y$, $n_{\xi}=3$ for $\xi=\varphi$. We recall here that $k_{i}$ is a set of integers composed of $\pqty{i,n_{i}}$, and represents a monomer of dimer $i$.

The elements for $i=j$ are

[TABLE]

and the elements for $i\neq j$ are

[TABLE]