## 3 [PERSON] et al. [25] set up a protocol by conducting isotropic compression tests in DEM to investigate the jamming of granular materials and proved their theoretical phase diagram. The following protocol to prepare specimens is used (it will show later that this method \begin{table}
\begin{tabular}{c c c c c c} \hline \(\rho_{grain}\) (kg/m3) & \(E\) (GPa) & \(\kappa\) & \(\beta_{n}\) & \(\beta_{t}\) & \(\mu\) \\ \hline
2650 & 0.5 & 0.5 & 0.2 & 0.2 & 0.5 \\ \hline \end{tabular}
\end{table} 
Table 1: DEM parameters. can only generate dense specimens): All following simulations in this study are performed in 3D periodic domains. (1) Grains are randomly inserted within a cuboidal domain (each side is 4 mm long) with the possibility of overlap until an initial solid fraction \(\phi_{0}\) (simply called density in this study) is achieved. The solid fraction is the volume of grains over the total volume. i.e., \(\phi=\frac{V_{g}}{V}=\frac{V_{g}}{V_{g}+V_{v}}\). Depending on the initial density, the number of grains ranges from 13,800 to 15,000. (2) Set the friction coefficient to zero (i.e., \(\mu=0\)) and grains are left to reach equilibrium. (3) Set the friction coefficient back to 0.5 (i.e., \(\mu=0.5\)) until equilibrium is reached again. The jamming, shearing, and critical state of granular soil with surface friction of \(\mu=0.5\) are studied. (4) The periodic domain shrinks isotopically, mimicking the isotropic compression tests. The deformation rate is fixed at a small value of 1 s\({}^{-1}\), making sure that the quasi-static deformation condition is met. During simulations, the stress is estimated from grains with [PERSON]'s equation as in the literature [13]  
\[\sigma=-\frac{1}{V}\sum_{g\in V}\sum_{c\in g}l_{gc}\otimes f_{gc} \tag{3}\]  
Here \(V\) is the volume of the representative element volume (REV), the contact vector \(l_{gc}\) is the vector pointing from the centre of gravity of grain \(g\) to the contact point \(c\). The summation is performed over all contacts (\(c\in g\)) of every grain (\(g\in V\)) in RVE. In the present study, compressive stress and strain are defined as positive. Moreover, this study focuses on normalised stress, i.e., \(\sigma^{\prime}=\frac{g}{E}\). The is the average number of contacts per grain, i.e., \(Z=\frac{M}{N}\), where \(N\) is the total number of grains, and \(M\) is the total number of contacts (two grains in contact have two contacts belonging to each grain). For granular assembly with uniform spherical grains, some grains are rattlers (having less than 4 contacts), and they do not contribute to the contact force chain [26], therefore, a corrected coordination number \(Z^{4}=\frac{M^{4}}{N^{4}}\) is used, where \(N^{4}\) is the total number of non-rattler grains, and \(M^{4}\) is the total number of non-rattler contacts. During each isotropic compression test, the normalised pressure \(p^{\prime}\), corrected the coordination number \(Z^{4}\), and the density \(\phi\) are recorded. Figure 1 gives two examples. A first assembly starts with an initial density of \(\phi_{0}=0.5233\). Initially, the assembly is unjammed with \(p^{\prime}=0\) and \(Z^{4}=4\). After it is compressed until larger than a jamming density of \(\phi_{j}=0.5779\), it becomes a jammed matter and then the pressure increases continuously with the increase of density. Figure 1b shows that the corrected coordination number also increases with density. Figure 1: Isotropic compression test to determine the jamming density. To determine the jamming density, the normalised pressure \(p^{\prime}\) and the density \(\phi\) are assumed to follow an equation [26] as  
\[p^{\prime}=\frac{\phi C}{\phi_{j}}p_{0}\log\left(\frac{\phi}{\phi_{j}}\right) \left[1-\gamma_{p}\log\left(\frac{\phi}{\phi_{j}}\right)\right] \tag{4}\]  
where \(p_{0}\), \(\gamma_{p}\), and the jamming density \(\phi_{j}\) are fitting parameters. Figure 1a shows that the DEM data are well-fitted by this equation. The assembly starting with \(\phi_{0}=0.5233\) has a jamming density of \(\phi_{j}=0.5779\).