<EMAIL_ADDRESS>
# Search for the production of dark matter in the framework of Mono-Z′ portal
at the ILC simulated electron-positron collisions at $\sqrt{s}=500$ GeV

S. Elgammal  [  Centre for theoretical physics, The British University in
Egypt.

(August 30, 2024)

###### Abstract

In the present work, we study the possible production of the light neutral
gauge boson (Z′) candidates, which originated from a simplified-model scenario
based on the Mono-Z′ model, in association with dark matter. This study has
been performed by studying events with dimuon plus missing transverse energy
produced in the simulated electron-positron collisions at the foreseen
International Linear Collider (ILC), operating at 500 GeV center of mass
energy and integrated luminosity of 100 fb-1. In case no new physics has been
discovered, we set upper limits at 95% confidence level on the masses of
various particles in the model as, spin-1 (Z′), as well as the fermionic dark
matter.

††preprint: APS/123-QED

##  I Introduction

The Standard Model of particle physics (SM) has succeeded in describing many
features of nature that we observe in our experiments [1]. The most famous
example, arguably, is the agreement between the SM prediction and the
experimental measurement of the electroweak theory, which led to the old
discovery of $W^{\pm}$, $Z^{0}$ bosons and the recent discovery of the SM
Higgs boson by ATLAS experiment [2] and CMS experiment [3].

Although its success, SM is well known that it is not a complete theory. It
describes only the visible, baryonic, matter in the Universe, which represents
about 5% of the energy density distribution content of the Universe. The rest
is attributed to the mysterious dark energy and dark matter (DM). Many
cosmological observations, based on recent observations by [4, 5, 6], support
the existence of dark matter. These observations suggest that dark matter is
non-decaying and weakly interacting massive particles. So, SM could be
considered as a low-energy manifestation of other theories realized at high
energy, generically known as BSM (Beyond the Standard Model) theories [7].

Both the ATLAS and CMS Collaborations have previously searched for the massive
extra neutral gauge boson Z′, which is induced by the Grand Unified Theory
(GUT) and Supersymmetry [8, 9, 10, 11], with no evidence of their existence
using the full RUN II period of the LHC data [12, 13]. Results by the CMS have
excluded, at 95% confidence level (CL), the existence of Z′ in the mass range
between 0.6 - 5.15 TeV, while the ATLAS experiment has excluded the mass range
between 0.6 - 5.1 TeV.

In addition, many searches for DM have been performed by analyzing the data
collected by the CMS and ATLAS experiments during RUN II. These searches rely
on the production of a visible object ”X”, which recoils against the large
missing transverse energy from the dark matter particles leaving a signature
of ($\text{X}+E^{\text{miss}}_{T}$) in the detector [14]. The visible particle
could be an SM particle like W, Z bosons or jets [15, 16], photon [17, 18] or
SM Higgs boson [19, 20, 21].

The foreseen International Linear Collider (ILC) is an electron-positron
collider, which has been proposed to operate at 500 GeV center of mass energy
($\sqrt{s}$) as a startup with possible upgrade to $\sqrt{s}=1000$ GeV in run
II [22, 23, 24, 25]. The linear electron-positron colliders are characterized
by controllable initial particle energy, low QCD background compared to hadron
colliders, and controllable beam polarization.

In this analysis, we present a search for light neutral gauge bosons (Z′)
(i.e. $M_{Z^{\prime}}<90$ GeV), which are originated from a simplified-model
scenario, which is known as the light vector (LV), in the framework of Mono-Z′
model [26], at the ILC simulated electron-positron collisions with 500 GeV
center of mass energy. The topology of the studied simulated events is dimuon,
from the decay of Z′, plus large missing transverse energy which is attributed
to dark matter. A similar analysis has been presented for searching for dark
matter in association with the visible particle being a $Z$ boson decaying to
dimuon at the ILC [27].

ATLAS experiment has performed a search for dark matter, in the framework of
Mono-Z′ model, with the muonic decay channel of Z′ at the LHC [28]. This
search has excluded Z′ masses between 200 and 1000 GeV, and seated separate
limits on the coupling of the Z′ to the SM leptons, $\texttt{g}_{l}$. Thus,
$\texttt{g}_{l}$ ranges from 0.01–0.025 and 0.02–0.38 for Z′ masses between
200 and 1000 GeV are excluded for the light and heavy dark-sector in the
light-vector scenario. A search, which has been done by CMS experiment in
final states with dimuon, excluded Z′ boson in the mass region below 200 GeV
[29]. In addition, the LHCb has placed tight constraints on the Z′ mass and
Z′-dimuon coupling in the dimuon mass region that the current paper probes
[30, 31].

If the Z′ does not couple to quarks, the HL-LHC and future hadron colliders
would provide no limits on the existence of Z′. Then the future electron-
positron colliders as ILC will play an important role.

In the following section II, the theoretical formalism of the simplified model
based on the Mono-Z′ portal and its free parameters are presented. Then the
simulation techniques used for events generation for the signal and SM
background samples are displayed in section III. Afterward, the selection cuts
and the strategy of the analysis are explained in section IV. Finally, the
results and the summary of this analysis are discussed in sections V and VI
respectively.

Figure 1: Feynman diagram for the Light Vector (LV) scenario based on Mono-Z′
model; for the production of neutral light gauge boson (Z′) in association to
dark matter ($\chi_{1}$) pair [26].

##  II The simplified model in the framework of Mono-Z′ portal

Our target model is known as Mono-Z′, which has been discussed in [26],
assumes the production of dark matters from electron-positron collisions at
the ILC through a new light gauge boson Z′. The dark matter can be produced
through three different possible scenarios in the Mono-Z′ model. In the
current analysis, we focus on one of these scenarios, which is known as the
Light Vector (LV) scenario and is presented in figure 1. The proposed dark
matter can be produced through the process of pair annihilation of electron-
positron mediated by a light vector boson Z′, which then undergoes two dark
matters, a light dark matter $(\chi_{1})$ and a heavy one $(\chi_{2})$.
$\chi_{2}$ is heavy enough to decay to a Z′ and another light dark matter
$\chi_{1}$ (i.e. $\chi_{2}$ $\rightarrow$ Z′ $\chi_{1}$) as shown in figure 1.

The interaction term, in the Lagrangian, between the dark fermions and Z′ is
given by [26]

$\frac{\texttt{g}_{DM}}{2}Z^{\prime}_{\mu}\large(\bar{\chi_{2}}\gamma^{\mu}\gamma^{5}\chi_{1}+\bar{\chi_{1}}\gamma^{\mu}\gamma^{5}\chi_{2}\large),$
(1)

where $\texttt{g}_{DM}$ is the coupling of Z′ to the dark matter $\chi_{1}$
and $\chi_{2}$. In the rest of the paper, the notation $(\chi)$ refers to
$(\chi_{1})$. For future reference, the coupling of the Z′ to the visible
leptons will be represented by $\texttt{g}_{l}$, while the coupling of the Z′
to the dark matter will be represented by $\texttt{g}_{DM}$ in the rest of
this paper.

Scenario | Masses assumptions
---|---
|  $M_{\chi_{1}}=1,5,...,200$ GeV
Light dark sector |
|  $M_{\chi_{2}}=M_{\chi_{1}}+M_{Z^{\prime}}+25$ GeV
Table 1: The Light mass assumptions for the dark sector for the light vector
scenario [26].

The only allowed decays in the LV scenario are assumed to be the decay of
${Z}^{\prime}\rightarrow\chi_{1}\chi_{2}$,
$\chi_{2}\rightarrow{Z^{\prime}}\chi_{1}$ and
${Z}^{\prime}\rightarrow\mu\tilde{\mu}$. Where the total decay width of the
${Z}^{\prime}$ and ${\chi_{2}}$ can be calculated given the values of the
masses of ${Z}^{\prime}$ and dark matter and the coupling constants.

The free parameters in this scenario are the lightest dark matter mass
$M_{\chi_{1}}$, $M_{\chi_{2}}$, the Z′ mass $(M_{Z^{\prime}})$ and the
coupling of Z′ to both leptons and dark matter particles $\texttt{g}_{l}$ and
$\texttt{g}_{DM}$ respectively. Based on the LEPII results [32] for
$M_{Z^{\prime}}<209$ GeV; $\texttt{g}_{l}\leq 0.04$ is excluded, so that we
have chosen $\texttt{g}_{l}=$ 0.05. While $\texttt{g}_{DM}=1.0$ according to
[26]. Meanwhile the values of the masses ($M_{Z^{\prime}}$, $M_{\chi_{1}}$ and
$M_{\chi_{2}}$) are not fixed but scanned over.

The signature that these processes leave at the detector consists of two
oppositely charged leptons, which are produced from the decay of Z′, in
addition to a large missing transverse energy from the stable dark matter
$\chi_{1}$. So our studied events are characterized by the following topology
($\mu^{+}\mu^{-}+E^{miss}_{T}$).

The intensive searches for Z′ over many years by CMS and ATLAS detectors have
proven that the absence of heavy neutral gauge bosons (Z′) in the mass range
between 0.2 and 5.15 TeV so that we will focus our analysis on the search for
light neutral gauge bosons (Z′) in mass range below 90 GeV. For this reason,
we consider the LV scenario with the use of light-dark sector case for
acquiring mass to dark matter particles ($\chi_{1}$ and $\chi_{2}$), which is
presented in table 1. This is a particular choice to fix $M_{\chi_{1}}$ and
$M_{\chi_{2}}$ according to a prescription given in [26]. However, this is
only one of the possible choices.

$M_{\chi_{1}}$ $M_{Z^{\prime}}$ | 10 | 20 | 30 | 40 | 50 | 60 | 70 | 80
---|---|---|---|---|---|---|---|---
1 | $2.15\times 10^{2}$ | $1.95\times 10^{2}$ | $1.87\times 10^{2}$ | $1.80\times 10^{2}$ | $1.74\times 10^{2}$ | $1.69\times 10^{2}$ | $1.67\times 10^{2}$ | $1.63\times 10^{2}$
5 | $2.12\times 10^{2}$ | $1.91\times 10^{2}$ | $1.81\times 10^{2}$ | $1.74\times 10^{2}$ | $1.69\times 10^{2}$ | $1.65\times 10^{2}$ | $1.62\times 10^{2}$ | $1.59\times 10^{2}$
10 | $2.11\times 10^{2}$ | $1.85\times 10^{2}$ | $1.75\times 10^{2}$ | $1.68\times 10^{2}$ | $1.63\times 10^{2}$ | $1.59\times 10^{2}$ | $1.55\times 10^{2}$ | $1.52\times 10^{2}$
25 | $2.02\times 10^{2}$ | $1.72\times 10^{2}$ | $1.60\times 10^{2}$ | $1.52\times 10^{2}$ | $1.46\times 10^{2}$ | $1.42\times 10^{2}$ | $1.38\times 10^{2}$ | $1.35\times 10^{2}$
50 | $1.85\times 10^{2}$ | $1.55\times 10^{2}$ | $1.37\times 10^{2}$ | $1.27\times 10^{2}$ | $1.21\times 10^{2}$ | $1.17\times 10^{2}$ | $1.12\times 10^{2}$ | $1.09\times 10^{2}$
100 | $1.75\times 10^{2}$ | $1.06\times 10^{2}$ | $8.76\times 10^{1}$ | $7.83\times 10^{1}$ | $7.18\times 10^{1}$ | $6.66\times 10^{1}$ | $6.20\times 10^{1}$ | $5.80\times 10^{1}$
150 | $7.22\times 10^{1}$ | $4.51\times 10^{1}$ | $3.62\times 10^{1}$ | $3.08\times 10^{1}$ | $2.67\times 10^{1}$ | $2.30\times 10^{1}$ | $1.96\times 10^{1}$ | $1.66\times 10^{1}$
180 | $2.43\times 10^{1}$ | $1.67\times 10^{1}$ | $1.29\times 10^{1}$ | $1.01\times 10^{1}$ | $7.72\times 10^{0}$ | $5.69\times 10^{0}$ | $3.94\times 10^{0}$ | $2.51\times 10^{0}$
190 | $1.41\times 10^{1}$ | $1.00\times 10^{1}$ | $7.51\times 10^{0}$ | $5.48\times 10^{0}$ | $3.77\times 10^{0}$ | $2.38\times 10^{0}$ | $1.33\times 10^{0}$ | $5.94\times 10^{-1}$
199 | $7.44\times 10^{0}$ | $5.53\times 10^{0}$ | $3.88\times 10^{0}$ | $2.51\times 10^{0}$ | $1.44\times 10^{0}$ | $6.77\times 10^{-1}$ | $2.29\times 10^{-1}$ | $4.74\times 10^{-2}$
200 | $6.86\times 10^{0}$ | $5.11\times 10^{0}$ | $3.55\times 10^{0}$ | $2.25\times 10^{0}$ | $1.25\times 10^{0}$ | $5.62\times 10^{-1}$ | $1.74\times 10^{-1}$ | $3.22\times 10^{-2}$
Table 2: The light vector scenario production cross section times branching
ratios (in fb) at Leading Order (LO) for different choices of the DM mass
$M_{\chi_{1}}$ (in GeV) and Z′ mass $M_{Z^{\prime}}$ (in GeV); for the light-
dark sector mass assumption with the following couplings constants
$\texttt{g}_{l}=0.05,\leavevmode\nobreak\ \texttt{g}_{DM}=1.0$. Given that the
polarized degrees of electron and positron beams are $P_{e^{-}}=0.8$,
$P_{e^{+}}=-0.3$ at the ILC with $\sqrt{s}$ = 500 GeV.

Figure 2: Dependence of cross sections for the signal process
$e^{+}e^{-}\rightarrow\chi\bar{\chi}Z^{\prime}$($Z^{\prime}\rightarrow\mu^{+}\mu^{-})$
induced by the light vector (LV) on the DM mass with $\sqrt{s}$ = 500 GeV and
the electron-positron beams polarization are $P_{e^{-}}=0.8$,
$P_{e^{+}}=-0.3$. The blue dots refers to $M_{Z^{\prime}}=20$ GeV, while the
red dots for $M_{Z^{\prime}}=80$ GeV.

Figure 3: Dependence of cross sections for the signal process
$e^{+}e^{-}\rightarrow\chi\bar{\chi}Z^{\prime}$($Z^{\prime}\rightarrow\mu^{+}\mu^{-})$
induced by the light vector (LV) on the $Z^{\prime}$ mass with $\sqrt{s}$ =
500 GeV and the electron-positron beams polarization are $P_{e^{-}}=0.8$,
$P_{e^{+}}=-0.3$. The blue dots refers to $M_{\chi_{1}}=1$ GeV, the red dots
for $M_{\chi_{1}}=100$ GeV and the green dots for $M_{\chi_{1}}=200$ GeV.

##  III Simulation of signal samples and SM backgrounds

The SM background processes yielding muon pairs in the signal region are
Drell-Yan (DY $\rightarrow\mu^{+}\mu^{-}$) production, the production of top
quark pairs ($\text{t}\bar{\text{t}}\rightarrow\mu^{+}\mu^{-}+2b+2\nu$) and
production of diboson ($W^{+}W^{-}\rightarrow\mu^{+}\mu^{-}+2\nu$,
$ZZ\rightarrow\mu^{+}\mu^{-}+2\nu$ and $ZZ\rightarrow 4\mu$).

The signal samples, for the LV scenario, and the corresponding SM background
processes have been generated using WHIZARD event generator 3.1.1 [33], the
ISR effect was included, interfaced to Pythia 6.24 for parton shower model and
hadronization [34] and DELPHES [35] for a fast detector simulation of ILD
detector. They were generated from electron-positron collisions at the ILC
with 500 GeV center of mass energy, which corresponds to the circumstances of
RUN I, given that the polarized degrees of electron and positron beams are
$P_{e^{-}}=0.8$, $P_{e^{+}}=-0.3$ respectively.

For the light vector scenario with the use of dark matter particles
($\chi_{1}$ and $\chi_{2}$) mass assumptions, which are summarized in table 1,
and $\texttt{g}_{l}=0.05,\leavevmode\nobreak\ \texttt{g}_{DM}=1.0$: table 2
indicates the production cross section times branching ratios at Leading Order
(LO) corresponding to several mass points of Z′ and DM. Given that the
polarized degrees of electron and positron beams are $P_{e^{-}}=0.8$,
$P_{e^{+}}=-0.3$ at the ILC with $\sqrt{s}$ = 500 GeV.

The production cross sections times branching ratios are plotted, in figure 2,
versus the dark matter mass $M_{\chi_{1}}$ for two Z′ mass points
($M_{Z^{\prime}}$ = 20 and 80 GeV). They are generated at the ILC with
$\sqrt{s}$ = 500 GeV and the electron-positron beams polarization are
$P_{e^{-}}=0.8$, $P_{e^{+}}=-0.3$ respectively. The behavior of the cross
sections times branching ratios is almost flat up to dark matter mass equal to
100 GeV, then it decreases rapidly. In figure 3, the cross sections times
branching ratios are illustrated versus the neutral gauge boson (Z′) masses
for three choices of the dark matter mass $(M_{\chi_{1}})$ = 1, 100 and 200
GeV. A strong drop in the cross sections times branching ratios versus the
mass of Z′ has been noticed as the mass of dark matter increases.

All Monte Carlo of the SM background samples, used in this analysis and their
corresponding cross sections, were calculated at leading order and listed in
table 3. Thus, the contributions of the signal samples and the SM background
processes have been estimated from the Monte Carlo simulations, at which they
are normalized to their corresponding cross-section and integrated luminosity
of 100 fb-1.

An ad-hoc flat 10% uncertainty is applied to cover all possible systematic
effects.

Process | Deacy channel | Generator | $\sigma\times\text{BR}\leavevmode\nobreak\ (\text{fb})$ | Order
---|---|---|---|---
DY | $\mu^{+}\mu^{-}$ | Whizard | 1767 | LO
$\text{t}\bar{\text{t}}$ | $\mu^{+}\mu^{-}+2b+2\nu$ | Whizard | 10.4 | LO
WW | $\mu^{+}\mu^{-}+2\nu$ | Whizard | 232.8 | LO
ZZ | $\mu^{+}\mu^{-}+2\nu$ | Whizard | 3.7 | LO
ZZ | $4\mu$ | Whizard | 0.5 | LO
Table 3: The simulated SM backgrounds generated from electron-positron
collisions at the ILC with the polarized degrees of electron and positron
beams are $P_{e^{-}}=0.8$, $P_{e^{+}}=-0.3$ at $\sqrt{s}=500$ GeV. Their
corresponding cross-section times branching ratio for each process and the
generation order are presented. Names of these MC samples and the used
generators are stated as well.

##  IV Event selection

The selection of event, for the analysis, has been designed to reconstruct a
final state with two muons with low transverse momentum $(p_{T})$ in
association with missing transverse energy accounting for the DM candidate.
The selection is made in the form of cuts applied on different kinematic
parameters. Each of the two muons should pass the following preliminary
selection:
$\bullet$ $p^{\mu}_{T}$ (GeV) $>10$,
$\bullet$ $|\eta^{\mu}|$ (rad) $<$ 2.5,
$\bullet$ $\text{IsolationVar}<0.1$,
”IsolationVar” represents the isolation cut in DELPHES software to reject
muons produced inside jets. In this cut, it is required that the scalar
$p_{T}$ sum of all muon tracks within a cone of $\Delta R=0.5$ around the muon
candidate, excluding the muon candidate itself, should not exceed 10% of the
$p_{T}$ of the muon. Thus, each event has been selected with two opposite
charge muons.

Figure 4: The measured dimuon invariant mass spectrum, after applying
preliminary selection(i) summarized in table 4, for the estimated SM
backgrounds and different choices of neutral gauge boson (Z′) masses generated
based on the LV scenario, with dark matter mass ($M_{\chi_{1}}=1$ GeV).

Figure 5: The distribution of the missing transverse energy, after applying
selection(i) summarized in table 4; for the expected SM backgrounds, and mass
of $Z^{\prime}$ ($M_{Z^{\prime}}$ = 50 GeV) produced by the LV scenario, with
mass of dark matter ($M_{\chi_{1}}=1$ GeV).

Figure 4 shows the distribution of the dimuon invariant mass; the cyan
histogram represents the Drell-Yan background, the yellow histogram stands for
the vector boson pair backgrounds (WW and ZZ), and the $t\bar{t}$ background
is represented by the red histogram. These histograms are stacked. While the
signals of the LV scenario, which have been generated with different masses of
the Z′ boson and fixing the dark matter mass ($M_{\chi_{1}}=1$ GeV), are
represented by different colored lines, and are overlaid. The corresponding
distribution of the missing transverse energy is presented in figure 5. It is
seen in figure 5 that the signal sample is overwhelmed by the backgrounds for
missing transverse energy less than 100 GeV. Thus, it is necessary to apply a
tighter cut to discriminate signals from SM backgrounds as will be explained
in the next paragraph.

Step | Variable | Requirements
---|---|---
|  $p^{\mu}_{T}$ (GeV) |  $>$ 10
Selection(i) |  $|\eta^{\mu}|$ (rad) |  $<$ 2.5
| $\Sigma_{i}p^{i}_{T}/p^{\mu}_{T}$ |  $<$ 0.1
Selection(ii) | Mass window (GeV) | $(0.9\times M_{Z^{\prime}})<M_{\mu^{+}\mu^{-}}<(M_{Z^{\prime}}+25)$
Table 4: Summary of cut-based final event selection used in the analysis.

In addition to the preliminary selection, an extra tighter cut has been
applied. This tight cut is related to the invariant mass of the dimuon, at
which we restricted the invariant mass of the dimuon to a small range around
the mass of the light gauge boson Z′, such that $(0.9\times
M_{Z^{\prime}})<M_{\mu^{+}\mu^{-}}<(M_{Z^{\prime}}+25)$. This cut has been
suggested by [26], which aims to suppress the SM backgrounds by restricting to
a mass window around the Z′ mass point. The full selection of the analysis is
summarized in table 4.

##  V Results

The shape-based analysis has been used based on the missing transverse energy
distributions ($E^{\text{miss}}_{T}$), which are good discriminate variables
since the signals distributions are characterized by relatively large
$E^{\text{miss}}_{T}$ values compared to the SM backgrounds. The distribution
of the missing transverse energy, after the application of the final event
selection listed in table 4, is illustrated in figure 6.

Figure 6: The distribution of the missing transverse energy, after applying
the final analysis selection(i+ii) listed in table 4, for the expected SM
background and one signal benchmark corresponding to the LV with
$M_{Z^{\prime}}=50$ GeV is superimposed.

The event yields passing the analysis selection(i) and final selection(i+ii),
for each of the SM backgrounds and the signal of simplified LV scenario, which
was generated with masses of light gauge boson $M_{Z^{\prime}}$ = 50 GeV, and
dark matter mass $M_{\chi_{1}}=1$ GeV; corresponding to an integrated
luminosity of 100 fb-1 are presented in table 5. Uncertainties include both
statistical and systematic components, summed in quadrature.

The significance of the signal over the background was computed using the
formula given in [37]

$S=2\times\big{(}\sqrt{N_{s}+N_{b}}-\sqrt{N_{b}}\big{)},$ (2)

where $N_{s}$ and $N_{b}$ are the number of signals and a total of SM
background events passing the selections (i+ii) listed in table 4. In figure
7, we depict the $5\sigma$ detection region (defined as $S\geq 5$) on the
$M_{\chi_{1}}$ \- $M_{Z^{\prime}}$ plane by taking $\sqrt{s}=$ 500 GeV,
integrated luminosity of $100\leavevmode\nobreak\ \text{fb}^{-1}$ for the LV
scenario, and fixing the values of the coupling constants
$\texttt{g}_{l}=0.05,\leavevmode\nobreak\ \texttt{g}_{DM}=1.0$.

Figure 7: $5\sigma$ detection region on the $M_{\chi_{1}}$ \- $M_{Z^{\prime}}$ plane for the $\chi\bar{\chi}Z^{\prime}$($Z^{\prime}\rightarrow\mu^{+}\mu^{-})$ production induced by the light vector scenario at the $\sqrt{s}=$ 500 GeV ILC with integrated luminosity of $100\leavevmode\nobreak\ \text{fb}^{-1}$. Process | No. of events passing (i) | No. of events passing (i+ii)
---|---|---
$\text{DY}\rightarrow\mu^{+}\mu^{-}$ | $87575.6\pm 8762.6$ | $388.3\pm 43.5$
$\text{t}\bar{\text{t}}\rightarrow\mu^{+}\mu^{-}+2b+2\nu$ | $629.3\pm 67.7$ | 111.8$\pm 15.4$
$\text{WW}\rightarrow\mu^{+}\mu^{-}+2\nu$ | $13602.8\pm 1365.3$ | $118.4\pm 16.1$
$\text{ZZ}\rightarrow\mu^{+}\mu^{-}+2\nu$ | $251.5\pm 29.7$ | $0.4\pm 0.7$
$\text{ZZ}\rightarrow 4\mu$ | $46.2\pm 8.2$ | $0.5\pm 0.7$
Sum Bkgs | $102105.5\pm 10215.5$ | $619.4\pm 66.8$
Signal of LV scenario | $11501.6\pm 1155.2$ | $11366.5\pm 1141.6$
(at $M_{Z^{\prime}}$ = 50 GeV and $M_{\chi_{1}}=1$ GeV) |  |
Table 5: Cut flow table: number of events passing the pre-selection (middle
column) and full selection (right column), as obtained from simulation for
backgrounds and signal with a luminosity of 100 fb-1 at $\sqrt{s}=500$ GeV.
The signal sample corresponds to the simplified-model scenario LV with
$M_{Z^{\prime}}$ = 50 GeV, $M_{\chi_{1}}=1$ GeV, $g_{DM}=1.0$ and
$g_{l}=0.05$; the total uncertainties, including the statistical and
systematic components, are quoted for the simulated signal and background
samples.

To make a statistical interpretation for our results, we performed a
statistical test based on the profile likelihood method, with the use of the
modified frequentist construction CLs [38, 39] used in the asymptotic
approximation [40] to derive exclusion limits on the product of signal cross
sections and branching fraction Br($Z^{\prime}$ $\rightarrow\mu\mu$) at 95%
confidence level.

The 95% upper limit on the cross-section times the branching ratio versus the
mass of Z′, for the LV simplified scenario based on Mono-Z′ model, is
presented in figure LABEL:figure:fig7, with the muonic decay of the Z′ and
coupling constant values of $\texttt{g}_{l}=0.05$ and $\texttt{g}_{DM}=1.0$
and dark matter mass $M_{\chi_{1}}=$ 1 GeV (fig.LABEL:figure:fig7.a), 100 GeV
(fig.LABEL:figure:fig7.b) and 190 GeV (fig.LABEL:figure:fig7.c), which are
represented by the black solid curves.

For the LV scenario, based on the Mono-Z′ model, the cross section times the
branching ratio limit is presented in figure 10 as a function of the
mediator’s masses Z′ and the masses of the dark matter $M_{\chi_{1}}$. The
region enclosed by the contour is excluded using the muonic decay channel of
the Z′ boson for the following benchmark scenario: $\texttt{g}_{l}=0.05$ and
$\texttt{g}_{DM}=1.0$.

Figure 10: The 95% CL upper limits on the product of the cross-section and
branching fraction from the inclusive search, for variations of pairs of the
LV scenario parameters ($M_{Z^{\prime}}$ and $M_{\chi_{1}}$). The filled
region indicates the upper limit. The dotted black curve indicates the
expected exclusions for the nominal Z′ cross-section.

##  VI Summary

The ILC electron-positron collider is an optimal machine for the possible
detection of particles from BSM since it delivers a clean signature of unknown
particles as dark matter, extra neutral gauge bosons, and Kaluza-Klein
excitation concerning the QCD background.

In this view, we have studied the effects of a simplified-model scenario,
which is known as light vector (LV), via dark matter pair production
associated with a Z′ boson at the ILC. The LV signal samples have been
simulated based on electron-positron collisions corresponding to the foreseen
ILC RUN I with 500 GeV center of mass energy, for an integrated luminosity of
100 fb-1. Results from the muonic decay mode of Z′ are discussed, with fixing
the values of the coupling constants to be $g_{DM}=1.0$, $g_{l}=0.05$. Given
that the polarized degrees of electron and positron beams are $P_{e^{-}}=0.8$,
$P_{e^{+}}=-0.3$ respectively at the ILC.

At the beginning of the ILC data taking with an integrated luminosity of 100
fb-1 at 500 GeV center of mass energy; the $5\sigma$ discovery of the Z′
boson, originated from the LV scenario, can be reached in the dimuon mass
spectrum between 35 and 80 GeV. Nevertheless, we found that the studied
signal’s benchmarks with the choice of lighter dark matter masses
$M_{\chi_{1}}<50$ GeV are overproducing dark matter at the ILC.

If this signal is not observed at the ILC, we set upper limits on the masses
of Z′ and dark matter $(\chi_{1})$ at the 95% CL for the charged muonic
channel decay of Z′. The most stringent limit is set for light vector scenario
with $\texttt{g}_{l}=0.05$ and $\texttt{g}_{DM}=1.0$, excluding $M_{\chi_{1}}$
= 199 GeV at $M_{Z^{\prime}}$ = 20 GeV, and up to $M_{Z^{\prime}}$ = 69.8 GeV
for $M_{\chi_{1}}$ = 1 GeV.

###### Acknowledgements.

The author of this paper would like to thank Tongyan Lin, co-author of [26]
for providing us with the UFO model files, helping us with generating the
signal events, and cross-checking the results.

## References

  * [1] F. Halzen and A. D. Martin, (1984), “Quarks And Leptons: An Introductory Course In Modern Particle Physics.” isbn: 0471887412, 9780471887416.
  * [2] The ATLAS Collaboration, Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC, Phys.Lett. B716 (2012) 1-29.
  * [3] The CMS Collaboration, Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC, Phys. Lett. B 716 (2012) 30.
  * [4] Planck Collaboration, Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594 (2016) A13 [arXiv:1502.01589] [iNSPIRE-HEP].
  * [5] Planck 2018 results. VI. Cosmological parameters. A&A 641, A6 (2020) arXiv:1807.06209 [astro-ph.CO].
  * [6] Lage, C and Farrar, G, The bullet cluster is not a cosmological anomaly, JCAP, vol. 2015, no. 2, 038. https://doi.org/10.1088/1475-7516/2015/02/038.
  * [7] P. Langacker, “The standard model and beyond”, (2010). isbn: 9781420079067.
  * [8] M. Cvetic and S. Godfrey, “Discovery and identification of extra gauge bosons”, arXiv:hep-ph/9504216.
  * [9] A. Leike, “The Phenomenology of extra neutral gauge bosons”, Phy. Rep. 317, 143 (1999) [arXiv:hepph/9805494].
  * [10] M. Cvetic, P. Langacker, and B. Kayser, “Determination of g-R / g-L in left-right symmetric models at hadron colliders”, Phys. Rev. Lett. 68 (1992) 2871.
  * [11] S. Dimopoulos and H. Georgi, “Softly Broken Supersymmetry And SU(5)”, Nucl. Phys. B 193 (1981) 150.
  * [12] CMS Collaboration, Search for resonant and nonresonant new phenomena in high-mass dilepton final state at $\sqrt{s}$ = 13 TeV, JHEP 07 (2021) 208 [arXiv:2103.02708v2] [hep-ex].
  * [13] ATLAS Collaboration, Search for high-mass dilepton resonances using 139 fb-1 of $pp$ collision data collected at $\sqrt{s}$ = 13 TeV with the ATLAS detector, Phys. Lett. B 796 (2019) 68.
  * [14] Krovi AniruDF, Low Ian and Zhang Yue, Broadening dark matter searches at the LHC: mono-X versus darkonium channels. JHEP 10 (2018) 026 [arXiv:1807.07972] [hep-ph].
  * [15] CMS Collaboration, Search for new physics in final states with an energetic jet or a hadronically decaying W or Z boson and transverse momentum imbalance at $\sqrt{s}$ = 13 TeV, Phys. Rev. D 97 (2018) 092005. [arXiv:1712.02345] [hep-ex].
  * [16] ATLAS Collaboration, Search for dark matter in events with a hadronically decaying vector boson and missing transverse momentum in pp collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector, JHEP 10 (2018) 180 [arXiv:1807.11471] [hep-ex].
  * [17] CMS Collaboration, Search for new physics in the monophoton final state in proton-proton collisions at $\sqrt{s}$ = 13 TeV, JHEP. 10 (2017) 073, [arXiv:1706.03794v2] [hep-ex].
  * [18] ATLAS Collaboration, Search for dark matter in association with an energetic photon in pp collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector, JHEP 02 (2021) 226, [arXiv:2011.05259v2] [hep-ex].
  * [19] CMS Collaboration, Search for dark matter particles produced in association with a Higgs boson in proton-proton collisions at $\sqrt{s}$ = 13 TeV, JHEP 03 (2020) 025, [arXiv:1908.01713v2] [hep-ex].
  * [20] ATLAS Collaboration, Search for dark matter produced in association with a Standard Model Higgs boson decaying into b-quarks using the full Run 2 dataset from the ATLAS detector, JHEP 11 (2021) 209, [arXiv:2108.13391v2] [hep-ex].
  * [21] ATLAS Collaboration, Search for dark matter in events with missing transverse momentum and a Higgs boson decaying into two photons in pp collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector, JHEP 10 (2021) 13, [arXiv:2104.13240v2] [hep-ex].
  * [22] J. Brau, Y. Okada, N.Walker, arXiv:0712.1950 [physics.acc-ph].
  * [23] A. Djouadi et al., arXiv:0709.1893 [hep-ph]
  * [24] N. Phinney, N. Toge, N. Walker, arXiv:0712.2361 [physics.acc-ph]
  * [25] T. Behnke et al., ArXiv e-prints (2007). arXiv:0712.2356 [physics.ins-det]
  * [26] Marcelo Autran, Kevin Bauer, Tongyan Lin and Daniel Whiteson, Mono-Z′: searches for dark matter in events with a resonance and missing transverse energy. Physical Review D 92 (2015) 035007 [arXiv:1504.01386] [hep-ph].
  * [27] Wan Neng, Song Mao, Li Gang, Ma Wen-Gan, Zhang Ren-You and Guo Jian-You, Searching for dark matter via mono-Z boson production at the ILC, Eur. Phys. J. C (2014) 74:3219.
  * [28] ATLAS Collaboration, Search for a new leptonically decaying neutral vector boson in association with missing transverse energy in proton-proton collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector, ATLAS-CONF-2023-045, 18 August 2023.
  * [29] CMS Collaboration, Search for a narrow resonance lighter than 200 GeV decaying to a pair of muons in proton-proton collisions at $\sqrt{s}$ = 13 TeV, Phys. Rev. Lett. 124, 131802 (2020), [arXiv:1912.04776] [hep-ex].
  * [30] LHCb Collaboration, Searches for low-mass dimuon resonances, JHEP 10 (2020) 156, [arXiv:2007.03923] [hep-ex].
  * [31] LHCb Collaboration, Search for dark photons produced in 13 TeV pp collisions, Phys. Rev. Lett. 120, 061801 (2018), [arXiv:1710.02867 [hep-ex].
  * [32] Arnab Dasgupta, P. S. Bhupal Dev, Tao Han, Rojalin Padhan, Si Wang, Keping Xie, Searching for heavy leptophilic $Z^{\prime}$: from lepton colliders to gravitational waves, JHEP 12 (2023) 011.
  * [33] Wolfgang Kilian, Thorsten Ohl, and Jurgen Reuter. WHIZARD—simulating multi-particle processes at LHC and ILC. The European Physical Journal C, 71(9), Sep 2011.
  * [34] T. Sjostrand, S. Mrenna and P. Skands, PYTHIA 6.4 Physics and Manual [arXiv:hep-ph/0603175].
  * [35] J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, M. Selvaggi, DELPHES 3, A modular framework for fast simulation of a generic collider experiment, JHEP 1402 (2014).
  * [36] I. Bozovic-Jelisavcic, S. Lukic, G. Milutinovic Dumbelovic, M. Pandurovic, I. Smiljanic, Luminosity measurement at ILC, JINST 8 P08012.
  * [37] W. S. Emam and S. A. Moussa, Search for Stueckelberg Z′ at LHC, Eur. Phys. J. C 60, 441–447 (2009). https://doi.org/10.1140/epjc/s10052-009-0948-8.
  * [38] A. L. Read, Presentation of search results: the CLs technique, J. Phys. G: Nucl. Part.Phys. 28 (2002) 2693, doi:10.1088/0954-3899/28/10/313.
  * [39] T. Junk, Confidence level computation for combining searches with small statistics, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, Volume 434, Issues 2–3, 1999, Pages 435-443, ISSN 0168-9002, https://doi.org/10.1016/S0168-9002(99)00498-2.
  * [40] G. Cowan et al., Asymptotic formulae for likelihood-based tests of new physics, Eur. Phys. J. C 71 (2011), p. 1554, doi: 10.1140/epjc/s10052-011-1554-0, arXiv: 1007.1727 [physics.data-an], Erratum: Eur. Phys. J. C 73 (2013) 2501.


# Testing mean-field theory for jamming of non-spherical particles: Contact
number, gap distribution, and vibrational density of states

Harukuni Ikeda<EMAIL_ADDRESS>Graduate School of Arts and
Sciences, The University of Tokyo 153-8902, Japan

###### Abstract

We perform numerical simulations of the jamming transition of non-spherical
particles in two dimensions. In particular, we systematically investigate how
the physical quantities at the jamming transition point behave when the shapes
of the particle deviate slightly from the perfect disks. For efficient
numerical simulation, we first derive an analytical expression of the gap
function, using the perturbation theory around the reference disks. Starting
from disks, we observe the effects of the deformation of the shapes of
particles by the $n$-th order term of the Fourier series $\sin(n\theta)$. We
show that the several physical quantities, such as the number of contacts, gap
distribution, and characteristic frequencies of the vibrational density of
states, show the power-law behaviors with respect to the linear deviation from
the reference disks. The power-law behaviors do not depend on $n$ and are
fully consistent with the mean-field theory of the jamming of non-spherical
particles. This result suggests that the mean-field theory holds very
generally for nearly spherical particles whose shape can be expressed by the
Fourier series.

††preprint: APS/123-Qed

##  I Introduction

The jamming transition is a phenomenon that granular materials suddenly get
finite rigidity at a certain density referred to as the jamming transition
point $\varphi_{J}$ van Hecke (2009); Liu and Nagel (2010). Near
$\varphi_{J}$, several physical quantities, such as energy, mechanical
pressure, shear modulus, and contact number, exhibit the power-law behaviors
O’Hern  _et al._ (2003). This implies that the jamming transition is a
critical phenomenon van Hecke (2009); Liu and Nagel (2010), such as the
second-order phase transition of equilibrium systems Nishimori and Ortiz
(2010).

The simplest model to study the jamming transition is the system consisting of
frictionless spherical particles. The systematic numerical simulations
confirmed that the critical exponents of several physical quantities do not
depend on the spatial dimensions $d$ for $d\geq 2$ O’Hern  _et al._ (2003);
Vågberg  _et al._ (2011); Charbonneau  _et al._ (2014), while different
exponents appear in a quasi-one-dimensional system Ikeda (2020). The
dimensional independence of the critical exponents suggests that the exact
values of the critical exponents can be calculated by considering the large
dimensional limit $d\to\infty$, where the mean-field theory becomes exact
Charbonneau  _et al._ (2014); Parisi  _et al._ (2020). This calculation has
been done by using the replica method, which was originally developed for the
spin-glasses, but is now widely used for many disordered systems Mézard  _et
al._ (1987); Nishimori (2001). Resultant critical exponents well agree with
the numerical results in $d\geq 2$ Charbonneau  _et al._ (2014); Parisi  _et
al._ (2020). Several other mean-field theories, such as the variational
argument Wyart  _et al._ (2005); Yan  _et al._ (2016), effective medium theory
DeGiuli  _et al._ (2014a, b), and random matrix theory Beltukov (2015); Ikeda
and Shimada (2020), are also successfully applied to derive the scaling laws
of the jamming transition.

Motivated by the success of the mean-field theories for the jamming of
spherical particles, we recently developed a mean-field theory of the jamming
of nearly spherical particles Brito  _et al._ (2018); Ikeda  _et al._ (2019).
The theory predicts that several physical quantities, such as the contact
number and gap distribution, exhibit the singular behaviors in the limit of
perfect spheres. The theoretical conjectures, however, have been confirmed
mainly for the breathing particles, which are believed to have the same
universality class of non-spherical particles Brito  _et al._ (2018). It is,
of course, desirable to perform a more direct test for non-spherical
particles. That is the purpose of this work.

The problem is that for particles of general shape, it is numerically
demanding and technically involved to calculate the gap function, which is the
minimal distance between two particles and used to judge if two particles are
overlapped or not Lu  _et al._ (2015). However, since we are interested in the
case where the particle shapes are close to perfect spheres, we can derive an
analytic form of the gap function by using a perturbation expansion from the
reference spheres Ikeda  _et al._ (2020a, b). This allows us to perform an
efficient numerical simulation.

We consider a particle system in two dimensions $d=2$. Starting from perfect
disks, we deform the shapes of particles by the $n$-th order term of the
Fourier series $\propto\sin(n\theta)$. We observe how this deformation affects
the physical quantities, such as the contact number, gap distribution
function, and characteristic frequencies of the vibrational density of states,
at $\varphi_{J}$. We find that the qualitative behaviors do not depend on $n$
and fully consistent with the mean-field predictions.

The organization of the paper is as follows. In Sec. II, we summarize the
previous results for jamming of frictionless spherical and non-spherical
particles. In Sec. III, we introduce the model and several important physical
quantities. In Sec. IV, we discuss the approximation to calculate the gap
function and interaction potential. In Sec. V, we discuss the numerical
algorithm to generate configurations at the jamming transition point. In Sec.
VI, we check the validity of our approximation by comparing the result of our
model for $n=2$ with a previous result of ellipses. In Sec. VII, we discuss
the behavior of $\varphi_{J}$ and the fraction of rattles. In Sec. VIII, we
discuss the scaling of the contact number at the jamming transition point. In
Sec. IX, we discuss the scaling of the gap function. In Sec. X, we discuss the
behavior of the vibrational density of states. Finally, in Sec. XI, we
summarize the results and conclude the work.

##  II Previous results

We here summarize the previous numerical and theoretical results for the
jamming transition of frictionless spherical and non-spherical particles.

###  II.1 Spherical particles

Frictionless spherical particles are isostatic at the jamming transition point
in the thermodynamics limit. We first explain the concept of the isostaticity.
A system is said to be isostatic if it satisfies the following condition:

$\displaystyle N_{c}=N_{f},$ (1)

where $N_{c}$ denotes the total number of constraints of the system, and
$N_{f}$ denotes the total number of degrees of freedom. For $N$ frictionless
spherical particles in $d$ dimensions, the total number of degrees of freedom
is

$\displaystyle N_{f}=Nd.$ (2)

For a jammed configuration generated by the isotropic compression, $N_{c}$ is
written as

$\displaystyle N_{c}=\frac{Nz}{2}+d,$ (3)

where $z$ denotes the number of contacts per particle. In Eq. (3), the first
term $Nz/2$ represents the total number of contacts, and the second term $d$
comes from the requirement for the positive bulk modulus Goodrich  _et al._
(2014). Therefore, if the system is isostatic, the contact number per
particles is

$\displaystyle z_{\rm iso}=2d-\frac{2d}{N}.$ (4)

The systematic numerical simulations revealed that frictionless spherical
particles are almost isostatic at $\varphi_{J}$, more precisely, the contact
number at $\varphi_{J}$, $z_{J}\equiv z(\varphi_{J})$, is Goodrich  _et al._
(2012)

$\displaystyle z_{J}-z_{\rm iso}\propto\frac{1}{N}.$ (5)

In particular, $z_{J}\to z_{\rm iso}$ in the thermodynamic limit $N\to\infty$
Bernal and Mason (1960); O’Hern  _et al._ (2003).

$z$ increases on increasing the packing fraction $\varphi$. Numerical studies
near $\varphi_{J}$ revealed that $z$ exhibits the following power-law behavior
O’Hern  _et al._ (2003); Goodrich  _et al._ (2012):

$\displaystyle\delta z\equiv z-z_{\rm iso}\propto\delta\varphi^{1/2},$ (6)

where $\delta\varphi=\varphi-\varphi_{J}$. Another interesting power-law
appears if one observes the gap distribution,

$\displaystyle g(h)=\frac{1}{N}\sum_{i<j}\Theta(-h_{ij})\delta(h-h_{ij}),$ (7)

where $\Theta(x)$ denotes the Heaviside step function, and $h_{ij}$ denotes
the gap function defined as

$\displaystyle h_{ij}=\left|\bm{r}_{i}-\bm{r}_{j}\right|-R_{i}-R_{j}.$ (8)

Here $\bm{r}_{i}$ and $R_{i}$ denote the position and radius of the $i$-th
particle, respectively. At $\varphi_{J}$ in the thermodynamic limit, $g(h)$
exhibits the power law for $h\ll 1$ Donev  _et al._ (2005); Charbonneau  _et
al._ (2012):

$\displaystyle g(h)\sim h^{-\gamma},$ (9)

where $\gamma\approx 0.41$ Charbonneau  _et al._ (2014, 2020). For
$\varphi>\varphi_{J}$, on the contrary, the power-law is truncated at finite
$h$ Charbonneau  _et al._ (2012); Franz  _et al._ (2017):

$\displaystyle
g(h)\sim\begin{cases}\delta\varphi^{-\mu\gamma}&h\ll\delta\varphi^{\mu}\\\
h^{-\gamma}&h\gg\delta\varphi^{\mu}\end{cases},$ (10)

where

$\displaystyle\mu=\frac{1}{2(1-\gamma)}.$ (11)

Figure 1:  Schematic picture of $D(\omega)$. (a) $D(\omega)$ of spherical
particles. (b) $D(\omega)$ of non-spherical particles.

Finally, in Fig. 1 (a), we show the schematic figure of the vibrational
density of states near $\varphi_{J}$. $D(\omega)$ exhibits a plateau down to
the characteristic frequency $\omega_{*}$. $\omega_{*}$ exhibits the following
power-law behavior O’Hern  _et al._ (2003); Wyart  _et al._ (2005):

$\displaystyle\omega_{*}\sim\delta z\sim\delta\varphi^{1/2}.$ (12)

Several mean-field theories, such as the variational argument Yan  _et al._
(2016); Wyart  _et al._ (2005), effective medium theory DeGiuli  _et al._
(2014a, b), and replica method Charbonneau  _et al._ (2014); Franz  _et al._
(2017), have been successful in reproducing the above scaling behaviors of the
jamming of frictionless spherical particles. This motivates us to develop the
mean-field theories for the jamming of non-spherical particles, which we shall
discuss in the next subsection.

###  II.2 non-spherical particles

The vibrational argument predicts the same scaling laws as Eqs. (6) and (12),
for a system satisfying the isostaticity Eq. (1) at $\varphi_{J}$ Yan  _et
al._ (2016). The same argument holds even for non-spherical particles. For
instance, frictionless dimers become isostatic at $\varphi_{J}$ Schreck  _et
al._ (2010), and a systematic numerical simulation indeed confirmed the same
scaling laws as spherical particles Shiraishi  _et al._ (2019, 2020). More
generally, for non-spherical particles consisting of spherical particles, such
as dimers, trimers, …, $n$-mers, one can show that the systems are isostatic
at $\varphi_{J}$ (see Appendix. A for a theoretical argument).

On the contrary, the deformation from the perfect sphere (or dimer, trimer,
…), in general, breaks the isostaticity. For instance, frictionless ellipsoids
Donev  _et al._ (2004, 2007); Zeravcic  _et al._ (2009); Mailman  _et al._
(2009); Schreck  _et al._ (2010, 2012), spherocylinders Williams and Philipse
(2003); Blouwolff and Fraden (2006); Azéma and Radjaï (2010); Marschall and
Teitel (2018), superballs Jiao  _et al._ (2010), superellipsoids Delaney and
Cleary (2010), and circulo-polygons VanderWerf  _et al._ (2018) have been
known to be hypostatic $N_{c}<N_{f}$ at $\varphi_{J}$. In the previous works,
we have constructed the mean-field theories and derived the scaling laws for
non-spherical particles, which are hypostatic at $\varphi_{J}$ Brito  _et al._
(2018); Ikeda  _et al._ (2019, 2020a). Below we summarize the main results.

The diameter $R$ of a non-spherical particle depends on the direction
$\bm{n}$. We shall write $R$ as

$\displaystyle R(\bm{n},\Delta)=R_{0}\left(1+\Delta F(\bm{n})\right),$ (13)

where $R_{0}$ denotes the radius of the reference sphere, $\Delta$ represents
the linear deviation from the reference sphere, and $F(\bm{n})$ characterizes
the particle shape. The mean-field theory predicts that for $\Delta\ll 1$, the
contact number at $\varphi_{J}$ behaves as

$\displaystyle\delta z_{J}=z_{J}(\Delta)-z_{\rm iso}\propto\Delta^{1/2},$ (14)

and $g(h)$ behaves as

$\displaystyle g(h)\sim\begin{cases}\Delta^{-\mu\gamma}&h\ll\Delta^{\mu}\\\
h^{-\gamma}&h\gg\Delta^{\mu}\end{cases}.$ (15)

In Fig. 1, we also show the schematic picture of $D(\omega)$ of non-spherical
particles. As in the case of spherical particles, $D(\omega)$ exhibits a
plateau down to $\omega_{*}$. In addition, $D(\omega)$ of non-spherical
particles has two additional peaks at $\omega_{0},\omega_{1}\ll\omega_{*}$.
The characteristic frequencies exhibit the power-law behaviors:

$\displaystyle\omega_{0}\sim\Delta^{1/2}\delta\varphi^{1/2},\
\omega_{1}\sim\Delta,\ \omega_{*}\sim\Delta^{1/2}.$ (16)

In particular, at $\varphi_{J}$, we get

$\displaystyle\omega_{0}=0,\ \omega_{1}\sim\Delta,\
\omega_{*}\sim\Delta^{1/2}.$ (17)

This suggests that the modes in the lowest band become zero modes at
$\varphi_{J}$. These zero modes are consequence of the hypostaticity
$N_{c}<N_{f}$ and referred to as the quartic modes Mailman  _et al._ (2009).
Eqs. (14), (15), and (17) imply that increasing $\Delta$ at $\varphi_{J}$
causes the same results as increasing $\delta\varphi$ of spherical particles.

So far, Eq. (14) has been checked only for systems with small system size
$N=100$, where the power-law is barely visible Brito  _et al._ (2018), Eq.
(15) has not been checked for non-spherical particles, and Eq. (17) has been
checked only for ellipsoids Zeravcic  _et al._ (2009). The purpose of this
work is to test the mean-field predictions for general shapes of non-spherical
particles at $\varphi_{J}$.

##  III Settings

Here we introduce the model and several important physical quantities.

###  III.1 Model

Figure 2: Particle shapes for various $n$ and $\Delta$. Black solid, red
dotted, and blue dashed lines denote polar plots of
$R_{i}(\theta)/R_{i}^{0}=1+\Delta F_{n}(\theta)$ for $\Delta=0$, $0.1$, and
$0.3$, respectively.

We consider a two dimensional system consisting of $N$ non-spherical
particles. The radius of a non-spherical particle depends on its angle
$\theta\in[0,2\pi)$. We shall write the radius of the $i$-th particle as

$\displaystyle R_{i}(\theta)=R_{i}^{0}\left[1+\Delta F(\theta)\right],$ (18)

where $R_{i}^{0}$ denotes the radius of the reference disk defined by

$\displaystyle R_{i}^{0}=\frac{1}{2\pi}\int_{0}^{2\pi}d\theta R_{i}(\theta).$
(19)

$\Delta$ represents the linear deviation from the reference disk, and the
function $F(\theta)$ characterizes the particle shape. $F(\theta)$ satisfies

$\displaystyle\int_{0}^{2\pi}d\theta F(\theta)=0.$ (20)

Since $F(\theta)$ is a periodic function, it is natural to express $F(\theta)$
by the Fourier series as

$\displaystyle F(\theta)=\sum_{n=1}^{\infty}a_{n}\sin(n\theta+b_{n}),$ (21)

where the constant term does not appear due to Eq. (20). Starting from the
reference disk $R_{i}(\theta)=R_{i}^{0}$, we want to investigate how the
$n$-th term of the Fourier series perturbs the physical quantities at
$\varphi_{J}$. For this purpose, we shall consider the following functional
form of $R_{i}(\theta)$:

$\displaystyle R_{i}(\theta)=R_{i}^{0}\left[1+\Delta F_{n}(\theta)\right],$
(22) $\displaystyle F_{n}(\theta)=\frac{\sin(n\theta)}{n}.$ (23)

When $\Delta=0$, the $i$-th particle is, of course, a disk of the radius
$R_{i}(\theta)=R_{i}^{0}$. On increasing $\Delta$, the shape of the $i$-th
particle gradually deviates from the reference disk. In Fig. 2, we illustrate
the corresponding particle shape for various $n$ and $\Delta$.

###  III.2 Volume fraction

We define the volume fraction $\varphi$ as

$\displaystyle\varphi=\frac{\sum_{i=1}^{N}S_{i}}{L^{2}},$ (24)

where $L$ denotes the linear distance of the system, and $S_{i}$ denotes the
surface of the $i$-th particle calculated as

$\displaystyle
S_{i}=\int_{0}^{2\pi}d\theta\int_{0}^{R_{i}(\theta)}rdr=\frac{(R_{i}^{0})^{2}}{2}\int_{0}^{2\pi}d\theta\left[1+\Delta
F_{n}(\theta)\right]^{2}.$ (25)

###  III.3 Asphericity

For a particle of general shape, it is not always straightforward to find out
the parameter corresponding to the linear deviation from the reference disk
$\Delta$. It is sometimes convenient to use the asphericity:

$\displaystyle A_{i}=\frac{P_{i}^{2}}{4\pi S_{i}},$ (26)

where $P_{i}$ denotes the perimeter calculated as

$\displaystyle
P_{i}=\int_{0}^{2\pi}d\theta\sqrt{R_{i}(\theta)^{2}+\left(\frac{dR_{i}(\theta)}{d\theta}\right)^{2}}d\theta.$
(27)

$A_{i}$ is calculated straightforwardly for particles of any shape.
Furthermore, it is possible to derive the scaling relation between $A_{i}$ and
$\Delta$, as follows. $A_{i}$ takes the minimal value $A_{i}=1$ for a disk
$\Delta=0$, and $A_{i}>1$ for a non-disk $\Delta\neq 0$. Therefore, we get

$\displaystyle
A_{i}(\Delta)=1+\frac{\Delta^{2}}{2}A_{i}^{\prime\prime}(\Delta)$ (28)
$\displaystyle\to\Delta\propto(A_{i}-1)^{1/2},$ (29)

where the linear order term does not appear, as $A_{i}(\Delta)$ has a minimum
at $\Delta=0$. Eq. (29) allows to convert the scaling laws for $\Delta$ Eqs.
(14,15,16) to that for $A_{i}$ Ikeda  _et al._ (2020a). By substituting Eq.
(22) into Eq. (26), one can see that $A_{i}$ does not depend on $i$,
$A_{i}=A$. For $\Delta\ll 1$, we get

$\displaystyle A$ $\displaystyle\approx 1+\frac{\Delta^{2}}{2\pi}\int
d\theta\left(F_{n}^{\prime}(\theta)^{2}-F_{n}(\theta)^{2}\right)$ (30)
$\displaystyle=1+\frac{\Delta^{2}}{2}\left(1-\frac{1}{n^{2}}\right).$ (31)

Note that $O(\Delta^{2})$ order term vanishes for $n=1$. This means that the
$n=1$ Fourier component only causes the translation of a particle for the
lowest order correction Tarama  _et al._ (2013), as illustrated in Fig. 2 (a).
Hereafter we only consider $n>1$. Later, we use Eq. (31) to compare our result
with the previous work of ellipsoids.

##  IV Interaction potential

We consider the harmonic potential O’Hern  _et al._ (2003):

$\displaystyle V_{N}=\sum_{i<j}\frac{h_{ij}^{2}}{2}\Theta(-h_{ij}),$ (32)

where $\Theta(x)$ denotes the Heaviside step function, and $h_{ij}$ denotes
the gap function, which is the minimal distance between particles $i$ and $j$.
In general, it is a non-trivial task to calculate $h_{ij}$ for general shapes
of non-spherical particles. Since we are interested in the scaling behaviors
for $\Delta\ll 1$, we calculate $h_{ij}$ by using the first order expansion
w.r.t $\Delta$, see Appendix. B for details of the calculation. After some
manipulations, we get

$\displaystyle h_{ij}$ $\displaystyle=h_{ij}^{(1)}+O(\Delta^{2}),$ (33)
$\displaystyle h_{ij}^{(1)}$
$\displaystyle=\left|\bm{r}_{i}-\bm{r}_{j}\right|-R_{i}^{0}-R_{j}^{0}$
$\displaystyle-\Delta\left[R_{i}^{0}F_{n}(\theta_{i}-\theta_{ij})+R_{j}^{0}F_{n}(\theta_{j}-\theta_{ji})\right],$
(34)

where $\bm{r}_{i}=\\{x_{i},y_{i}\\}$ denotes the position of the $i$-th
particle, $\theta_{i}$ denotes the direction of the $i$-th particle, and
$\theta_{ij}$ denotes the relative angle of particles $i$ and $j$, namely,
$\theta_{ij}={\rm atan2}(y_{j}-y_{i},x_{j}-x_{i})$, see Fig. 3.

Figure 3: Definitions of $\bm{r}_{i}$, $\bm{r}_{j}$, $\theta_{i}$,
$\theta_{j}$, $\theta_{ij}$, and $h_{ij}^{(1)}$. Red arrows denote particles
directions, black arrow denotes the vector $\bm{r}_{j}-\bm{r}_{i}$, and dashed
lines are parallel to the $x$-axis.

Hereafter, we use Eq. (34) to calculate the interaction potential. Obviously,
the approximation holds only for $\Delta\ll 1$. Another limitation is that by
construction, our approximation allows that two particles have at most a
single contact, which is not true for particles of non-convex shapes, such as
dimers, even for small aspect ratios Schreck  _et al._ (2010).

##  V Numerics

We perform numerical simulations for $N$ particles confined in a $L\times L$
box. We impose the periodic boundary conditions for the both $x$ and $y$
directions. To avoid the crystallization, we consider an equimolar binary
mixture: $R_{i}^{0}=0.5$ for $i=1,\dots,N/2$ and $R_{i}^{0}=0.7$ for
$i=N/2+1,\dots,N$. $\varphi_{J}$ is the point at which $V_{N}$ begins to have
a finite value. In practical, we define $\varphi_{J}$ as a packing fraction
satisfying O’Hern  _et al._ (2003)

$\displaystyle 10^{-16}<V_{N}/N<2\times 10^{-16}.$ (35)

Here we explain how to generate configurations at $\varphi_{J}$. We first
generate a random initial configuration at small packing fraction
$\varphi=0.1$. Next, we slightly increase the density
$\varphi\to\varphi+\varepsilon$ with $\varepsilon=10^{-3}$, and then minimize
$V_{N}$ by using the FIRE algorithm, which combines the standard molecular
dynamics of the Newton equation

$\displaystyle m\frac{d^{2}\bm{r}_{i}}{dt^{2}}=-\frac{\partial
V}{\partial\bm{r}_{i}},$ $\displaystyle
I\frac{d^{2}\theta_{i}}{dt^{2}}=-\frac{\partial V}{\partial\theta_{i}},$ (36)

with adaptive damping of the velocity Bitzek  _et al._ (2006). We find that
the FIRE converges in a reasonable time if we set $m=1$ and $I=\Delta$ so that
$\ddot{\bm{r}}_{i}=O(\Delta^{0})$ and $\ddot{\theta}_{i}=O(\Delta^{0})$. We
stop the FIRE algorithm when $\left(\frac{\partial
V_{N}}{\partial\bm{r}_{i}}\right)^{2}+\left(\frac{\partial
V_{N}}{\partial\theta}_{i}\right)^{2}<10^{-25}$, or $V_{N}/N<10^{-16}$ O’Hern
_et al._ (2003). We repeat the above compression
$\delta\varphi\to\varphi+\varepsilon$ and minimization protocols as long as
$V_{N}/N\leq 10^{-16}$ after the minimization. On the contrary, if
$V_{N}/N>10^{-16}$ after the minimization, we then decompress the system by
changing the sign and amplitude of $\varepsilon$ as
$\varepsilon\to-\varepsilon/2$. We repeat the above compression/decompression
protocols by changing $\varepsilon\to-\varepsilon/2$ every time the energy
crosses the threshold value $V_{N}/N=10^{-16}$. We terminates the simulation
when $10^{-16}<V_{N}/N<2\times 10^{-16}$. In Fig. 4, we show configurations at
$\varphi_{J}$ generated by the above algorithm.

Figure 4:  Configurations at $\varphi_{J}$ for $N=32$, $\Delta=0.2$ and
$n=2,\dots,5$.

When calculate $z_{J}$ and $g(h)$, we remove the rattles, for which the
contact number is less than $d+1=3$ O’Hern  _et al._ (2003). To improve the
statistics, we take the average for $50$ independent samples.

##  VI Comparison of our model and ellipses

Figure 5:  Comparison of ellipses and our model for $n=2$. Data for ellipses
are taken from Ref. VanderWerf  _et al._ (2018).

In this section, we compare our results with a previous numerical simulation
of ellipses VanderWerf  _et al._ (2018), where $\varphi_{J}$ and $z_{J}$ were
calculated as functions of $A$. The shape of an ellipse is defined by the
following equation:

$\displaystyle 1=\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}.$ (37)

When the aspect ratio is close to one, say $a=R^{0}$ and
$b=R^{0}(1-\varepsilon)$ with $\varepsilon\ll 1$, we get

$\displaystyle R(\theta)$ $\displaystyle\equiv\sqrt{x^{2}+y^{2}}\approx
R^{0}\sqrt{1+2\varepsilon y^{2}}$ $\displaystyle\approx
R^{0}\left[1+\varepsilon+\varepsilon\sin\left(2\theta-\frac{\pi}{2}\right)\right].$
(38)

As the constant term and phase shift do not change the shapes of particles,
Eq. (38) suggests that ellipses can be identified with our model of $n=2$ in
the lowest order. In Fig. 13, we compare our result for $n=2$ and $N=512$ with
that of ellipses for $N=480$ VanderWerf  _et al._ (2018) (the finite $N$
effect is not visible in the semi-log plot.). We find a reasonable agreement
for $A\sim 1$ or equivalently $\Delta\ll 1$, as expected.

##  VII Packing fraction and fraction of rattles

Figure 6: $\Delta$ dependence of $\varphi_{J}$ for $N=1024$.

In Fig. 6, we show the $\Delta$ dependence of $\varphi_{J}$ for $n=2,\dots,5$.
On increasing of $\Delta$, $\varphi_{J}$ increases for all $n$. This is
consistent with the previous numerical results of non-spherical particles,
such as ellipsoids Donev  _et al._ (2004), dimers Schreck  _et al._ (2010),
spherocylinders Williams and Philipse (2003), and circulo-polygons VanderWerf
_et al._ (2018). The more gradual increase of $\varphi_{J}$ is observed for
the larger $n$. It is interesting future work to see what will happen for much
larger value of $n$. For $n\to\infty$, we expect that $\varphi_{J}$ decreases
on increasing of $\Delta$, because the surface of particles becomes very rough
and the system may be mapped into frictional particles of the effective
friction coefficient
$\mu\sim\left|\partial_{\theta}V_{N}/\partial_{\bm{r}}V_{N}\right|\sim\Delta$
Ikeda  _et al._ (2020b).

Figure 7: $\Delta$ dependence of the fraction of rattles $f_{r}$ for $N=1024$.

In Fig. 7, we show the $\Delta$ dependence of the fraction of rattles:

$\displaystyle f_{r}=\frac{N_{r}}{N},$ (39)

where $N_{r}$ denotes the number of rattles. We find that $f_{r}$ decreases on
increasing of $\Delta$ for all $n$. The similar results have been reported for
ellipses and ellipsoids Donev  _et al._ (2007).

##  VIII Contact number at jamming

In this section, we present our numerical results for the contact number at
the jamming transition point $z_{J}\equiv z(\varphi_{J})$.

###  VIII.1 $N$ dependence for $n=2$

We first perform the finite size scaling analysis for $n=2$. In Fig. 8, we
show $\delta z_{J}=z_{J}-z_{\rm iso}$ of our model for $n=2$ and
$N=64,\dots,1024$. We find a power-law region $\delta z_{J}\sim\Delta^{1/2}$
for intermediate values of $\Delta$. The power-law region becomes wider on
increasing $N$.

Figure 8: $\delta z_{J}$ for $n=2$ and several $N$. (a) Markers denote the
numerical results. (b) Scaling plot for the same data. Black solid line and
gray dashed line denote $\delta z_{J}\propto\Delta^{1/2}$ and $\delta
z_{J}=2/N_{\rm nr}$, respectively

Inspired by the finite $N$ scaling analysis for frictionless spherical
particles Goodrich  _et al._ (2012), we assume the following scaling form:

$\displaystyle\delta z_{J}\sim N_{\rm nr}^{-1}\tilde{z}(N_{\rm
nr}^{2}\Delta),$ (40)

where $N_{\rm nr}=N-N_{r}$ denotes the number of non-rattler particles, and

$\displaystyle\tilde{z}(x)\sim\begin{cases}x^{0}&x\ll 1\\\ x^{1/2}&x\gg 1\\\
\end{cases}.$ (41)

In Fig. 8 (b), we test the above scaling. A good scaling collapse confirms Eq.
(40). Also, we find that $\delta z_{J}\to 2/N_{\rm nr}$ for $N_{\rm
nr}^{2}\Delta\to 0$, see the dashed horizontal line in Fig. 8 (b). This means
that the system has just one extra contact than the number of degrees of
freedom, which is also consistent with the previous finite size analysis of
frictionless disks Goodrich  _et al._ (2012).

###  VIII.2 $n$ dependence for $N=1024$

Figure 9: $\delta z_{J}$ for $N=1024$ and $n=2,\dots,5$. Markers denote the
numerical results, and solid lines denote $\delta z_{J}\propto\Delta^{1/2}$.

Now we focus on the data of the largest size $N=1024$. We only show the
results for $\delta z_{J}\gg 1/N$ so that the finite $N$ effects do not
appear. In Fig. 9, we plot our numerical results of $\delta z_{J}$ for
$n=2,\dots,5$. We find that $\delta z_{J}\propto\Delta^{1/2}$ for all $n$,
which confirms the mean-field prediction Eq. (14).

##  IX Gap distribution

Figure 10:  (a–d) $Z(h)$ for $N=1024$, and $n=2,\dots,5$. Markers denote the
numerical results, solid line denotes $Z(h)\sim h^{1-\gamma}$, and dashed line
denotes $Z(h)\sim h$. (e–h) Scaling plot for the same data.

In this section, we discuss the gap distribution $g(h)$ at $\varphi_{J}$. To
improve the statistics, instead of $g(h)$ itself, we observe the cumulative
distribution function:

$\displaystyle Z(h)=\frac{\int_{0}^{h}g(h)dh}{\int_{0}^{h_{\rm cut}}g(h)dh}.$
(42)

By definition $Z(0)=0$ and $Z(h_{\rm cut})=1$. We set $h_{\rm cut}=2$, which
is large enough to observe the scaling behavior. In Fig. 10 (a–d), we show our
numerical results of $Z(h)$ for $n=2,\dots,5$. We find that for small $\Delta$
and $h$, $Z(h)$ exhibits the power-law $Z(h)\sim h^{1-\gamma}$, suggesting
$g(h)\sim h^{-\gamma}$. On the contrary, for large $\Delta$, $Z(h)$ exhibits
the liner behavior $Z(h)\sim h$ for $h\ll 1$, suggesting $g(h)\sim h^{0}$.
These results are consistent with the mean-field prediction Eq. (15).

By using Eq. (15), we can deduce the scaling form of $Z(h)$ as Ikeda  _et al._
(2020a)

$\displaystyle Z(h)\sim\Delta^{1/2}\tilde{Z}(\Delta^{-\mu}h),$ (43)

where $\tilde{Z}(x)$ satisfies

$\displaystyle\tilde{Z}(x)\sim\begin{cases}x&x\ll 1\\\ x^{1-\gamma}&x\gg
1\end{cases}.$ (44)

In Fig. 10 (e–h), we test the above equation. We find a reasonable data
collapse.

##  X Vibrational density of states

Finally, we investigate the vibrational density of states $D(\omega)$ at
$\varphi_{J}$. We define the Hessian of the interaction potential as

$\displaystyle\mathcal{H}_{X_{i}Y_{j}}=\frac{\partial^{2}V_{N}}{\partial
X_{i}\partial Y_{j}}=K_{X_{i}Y_{j}}+T_{X_{i}Y_{j}},$ $\displaystyle
K_{X_{i}Y_{j}}=\sum_{i<j}v^{\prime\prime}(h_{ij})\frac{\partial
h_{ij}}{\partial X_{i}}\frac{\partial h_{ij}}{\partial Y_{j}},$ $\displaystyle
T_{X_{i}Y_{j}}=\sum_{i<j}v^{\prime}(h_{ij})\frac{\partial^{2}h_{ij}}{\partial
X_{i}\partial Y_{j}},$ (45)

where $X_{i}\in\\{\bm{r}_{i},\theta_{i}\\}$ and
$Y_{j}\in\\{\bm{r}_{j},\theta_{j}\\}$. At the jamming transition point,
$v^{\prime}(h_{ij})=0$, and thus

$\displaystyle\mathcal{H}_{X_{i}Y_{j}}\to K_{X_{i}Y_{j}}$
$\displaystyle=\delta_{ij}\sum_{k\neq i}\Theta(-h_{ik})\frac{\partial
h_{ik}}{\partial X_{i}}\frac{\partial h_{ik}}{\partial
Y_{i}}+(1-\delta_{ij})\Theta(-h_{ij})\frac{\partial h_{ij}}{\partial
X_{i}}\frac{\partial h_{ij}}{\partial Y_{j}}.$ (46)

Using the eigenvalues of $H_{X_{i}Y_{j}}$, $\\{\lambda_{n}\\}_{n=1,\dots,3N}$,
$D(\omega)$ is calculated as

$\displaystyle
D(\omega)=\frac{1}{3N}\sum_{n=1}^{3N}\delta(\omega-\sqrt{\lambda_{n}}).$ (47)

As mentioned below Eq. (17), $D(\omega)$ has zero modes at $\varphi_{J}$, in
addition to the trivial zero modes related to the rattler particles Mailman
_et al._ (2009); VanderWerf  _et al._ (2018). In practice, however, the zero
modes have finite frequencies depending on the accuracy of the numerical
simulation. Hereafter, we focus on the range $\omega>10^{-5}$, which is large
enough to remove the zero modes.

Figure 11: $D(\omega)$ for $N=1024$ and $n=2,\dots,5$. Here we do not show the
zero modes.

In Fig. 11, we show our numerical results for $D(\omega)$. We find that the
behavior of the high $\omega$ region ($\omega>0.1$) does not much depend on
$\Delta$. On decreasing $\Delta$, $D(\omega)$ develops a plateau down to the
characteristic frequency $\omega_{*}$. $D(\omega)$ has the separated band at
$\omega_{1}\ll\omega_{*}$. These results are consistent with the mean-field
prediction shown in Fig. 1 (b). Note that the lowest band in Fig. 1 (b) does
not appear, since $\omega_{0}=0$ at $\varphi_{J}$, and we do not show the zero
modes.

Figure 12: $\Delta$ dependence of characteristic frequencies. (a) Markers
denote numerical results of $\omega_{1}$, while solid lines denote
$\omega_{1}\sim\Delta$. (b) (a) Markers denote numerical results of
$\omega_{*}$, while solid lines denote $\omega_{*}\sim\Delta^{1/2}$.

We want to calculate $\omega_{1}$ and $\omega_{*}$ from the numerical data of
$D(\omega)$. For this purpose, we define $\omega_{1}$ as the point that
maximizes $D(\omega)$, and $\omega_{*}$ as the point where $D(\omega_{*})=0.1$
in the range $\omega<0.1$. In Fig. 12, we show $\Delta$ dependence of
$\omega_{1}$ and $\omega_{*}$. We find $\omega_{1}\sim\Delta$ and
$\omega_{*}\sim\Delta^{1/2}$, which are consistent with the mean-field
prediction, Eq. (17). The similar results have been previously reported for
ellipses and ellipsoids Schreck  _et al._ (2012); Brito  _et al._ (2018).

##  XI Summary and discussions

In this work, we performed a systematic numerical investigation for the
jamming of nearly spherical particles in two dimensions. Starting from perfect
disks, we systematically deformed the shapes of particles by the $n$-th order
term of the Fourier series $\propto\sin(n\theta)$ and observed its effects on
the physical quantities at the jamming transition point. For an efficient
numerical simulation, we derived an analytic formula of the gap function by
using the perturbation expansion from the reference disks. By using the
approximated gap function, we numerically generated configurations at the
jamming transition point, and calculated the contact number, (cumulative) gap
distribution, and vibrational density of states for $n=2,\dots 5$. We found
the qualitatively the same scaling behaviors, which are fully consistent with
the mean-field predictions, for all $n$. This means that mean-field prediction
is applicable to general-shaped convex particles whose particle shape can be
represented by the Fourier series.

There are still several important points that deserve further investigation.
Here we give a tentative list:

  * •

As mentioned before, our approximation does not hold for non-convex particles,
such as dimers, where particles may have multiple contacts. It is important
future work to extend the approximation for the gap function so as to take
into account the effects of the multiple contacts.

  * •

In this work, we investigate the physical quantities only at $\varphi_{J}$. It
is of course important to investigate the behavior for $\varphi>\varphi_{J}$.
For instance, the mean-field theory predicts that the shear modulus $G$
behaves as

$\displaystyle G$
$\displaystyle\sim\Delta^{1/2}\tilde{G}\left(\Delta^{-1}\delta\varphi\right),$
(48)

where $\tilde{G}(x)\sim x$ for $x\ll 1$ and $\tilde{G}(x)\sim x^{1/2}$ for
$x\gg 1$ Ikeda  _et al._ (2020a). So far the above scaling is confirmed only
for ellipsoids Ikeda  _et al._ (2020a). It is important to test if the same
scaling holds for other shapes of particles.

  * •

The variational argument predicts that the correlation volume $v_{\rm corr}$
behaves as $v_{\rm corr}\sim\left|\delta z^{-1}\right|$ Yan  _et al._ (2016).
For frictionless spherical particles, $\delta z=0$, thus $v_{\rm corr}$
diverges at $\varphi_{J}$. On the contrary, for non-spherical particles,
$\delta z\neq 0$, therefore $v_{\rm corr}$ remains finite even at
$\varphi_{J}$. Recently, it has been reported that $v_{\rm corr}$ can be
extracted from the participation ratio of the lowest frequency mode of the
vibrational density of states Shimada  _et al._ (2018). It is interesting to
repeat the same analysis for non-spherical particles.

  * •

In this work, we focus on a system in two dimensions $d=2$. It is important
future work to extend the current approximation and analysis to higher $d$.

  * •

The mean-field theory of non-spherical particles predicts that the replica
symmetry breaking (RSB) occurs near the jamming transition point Ikeda  _et
al._ (2019), as in the case of spherical particles Charbonneau  _et al._
(2014). It is important future work to find out the signature of the RSB for
non-spherical particles by numerical simulations and experiments.

  * •

For frictionless spherical particles, the different critical exponents appear
in the quasi-one-dimensional system Ikeda (2020). It is interesting future
work to repeat the similar analysis in Ref. Ikeda (2020) for non-spherical
particles.

###### Acknowledgements.

This project has received funding from the JSPS KAKENHI Grant Number
JP20J00289.

##  Appendix A Isostaticity of particles consisting of spherical particles

To keep the generality, we consider $N$ particle system connected by $M$
bonds. For instance, $N/2$ dimers can be considered as $N$ spherical particles
with $M=N/2$ bonds. We consider the harmonic potential:

$\displaystyle
V_{N}=\sum_{i<j}^{1,N}\frac{h_{ij}^{2}}{2}\Theta(-h_{ij})+k\sum_{a=1}^{M}\frac{t_{i_{a}j_{a}}^{2}}{2},$
(49)

where

$\displaystyle h_{ij}=\left|\bm{r}_{i}-\bm{r}_{j}\right|-R_{i}-R_{j},$
$\displaystyle
t_{i_{a}j_{a}}=\left|\bm{r}_{i_{a}}-\bm{r}_{j_{a}}\right|-l_{i_{a}j_{a}}.$
(50)

$\bm{r}_{i}=\\{x_{1}^{i},\dots,x_{d}^{i}\\}$ denotes the position, and $R_{i}$
denotes the diameter of the $i$-th particle, and $l_{i_{a}j_{a}}$ denotes the
length of the $a$-th bond connecting particles $i_{a}$ and $j_{a}$.

Here we show that the system is isostatic at $\varphi_{J}$ by using the same
argument for frictionless spherical particles Wyart  _et al._ (2005); Wyart
(2005).

The number of degrees of freedom of the system is

$\displaystyle N_{f}=Nd.$ (51)

At the jamming transition point, we have

$\displaystyle\left|\bm{r}_{i_{\mu}}-\bm{r}_{j_{\mu}}\right|=R_{i_{\mu}}+R_{i_{\mu}},$
$\displaystyle\mu=1,\dots,N_{c},$
$\displaystyle\left|\bm{r}_{i_{a}}-\bm{r}_{j_{a}}\right|=l_{i_{a}j_{a}},$
$\displaystyle a=1,\dots,M.$ (52)

where $N_{c}$ denotes the number of contacts, $i_{\mu}$ and $j_{\mu}$ denote
particles of the $\mu$-th contact. One can find
$\\{\bm{r}_{i}\\}_{i=1,\dots,N}$ satisfying the above equation if

$\displaystyle Nd\geq N_{\rm const},$ (53)

where

$\displaystyle N_{\rm const}=N_{c}+M,$ (54)

denotes the number of constraints at the jamming transition point. On the
contrary, the force balance requires

$\displaystyle\frac{\partial V_{N}}{\partial\bm{r}_{i}}=\sum_{j\neq
i}\Theta(-h_{ij})h_{ij}\bm{n}_{ij}+k\sum_{a=1}^{M}t_{ij_{a}}\bm{n}_{ij_{a}}=0,$
(55)

where $\bm{n}_{ij}$ denotes the normal vector connecting particles $i$ and
$j$. This can be regarded as $Nd$ linear equations for
$\\{h_{i_{\mu}j_{\mu}}\\}_{\mu=1,\dots,N_{c}}$ and
$\\{t_{i_{a}j_{a}}\\}_{a=1,\dots,M}$. One can find a solution if

$\displaystyle N_{c}+M\geq Nd.$ (56)

From Eqs. (53) and (56), we get

$\displaystyle N_{\rm const}=N_{f}\leftrightarrow N_{c}+M=Nd,$ (57)

meaning that the system is isostaticity at the jamming transition point. For
$N/2$ dimers, the total number of contacts is written as $N_{c}=(N/2)z_{J}/2$,
leading to

$\displaystyle z_{J}=4d-2,$ (58)

which is consistent with the numerical results in $d=2$ Schreck  _et al._
(2010); Shiraishi  _et al._ (2019) and $d=3$ Shiraishi  _et al._ (2020).

##  Appendix B Derivation of Eq. (34)

Figure 13:  Schematic picture of two non-spherical particles. Red solid lines
denote particles shape, and blue dashed lines denote reference disks.

We write the gap function $h_{ij}$ as

$\displaystyle h_{ij}$ $\displaystyle=\left|\bm{u}_{i}-\bm{u}_{j}\right|,$
(59)

where $\bm{u}_{i}$ and $\bm{u}_{j}$ are points on the surfaces of particles
$i$ and $j$ that minimize $h_{ij}$, see Fig. 13. We expand $\bm{u}_{i}$ and
$\bm{u}_{j}$ from those of the reference disks as

$\displaystyle\bm{u}_{i}=\bm{r}_{i}+R_{i}(\bm{n}_{i})\bm{n}_{i}=\bm{u}_{i}^{0}+\delta\bm{u}_{i},$
$\displaystyle\bm{u}_{j}=\bm{r}_{j}+R_{j}(\bm{n}_{j})\bm{n}_{j}=\bm{u}_{j}^{0}+\delta\bm{u}_{j},$
(60)

where

$\displaystyle\bm{u}_{i}^{0}=\bm{r}_{i}+R_{i}^{0}\bm{n}_{i}^{0},\
\bm{u}_{j}^{0}=\bm{r}_{j}+R_{j}^{0}\bm{n}_{j}^{0},$
$\displaystyle\delta\bm{u}_{i}=R_{i}(\bm{n}_{i})\bm{n}_{i}-R_{i}^{0}\bm{n}_{i}^{0},\delta\bm{u}_{j}=R_{j}(\bm{n}_{j})\bm{n}_{j}-R_{j}^{0}\bm{n}_{j}^{0},$
$\displaystyle\bm{n}_{i}=\frac{\bm{r}_{i}-\bm{u}_{i}}{\left|\bm{r}_{i}-\bm{u}_{i}\right|},\bm{n}_{j}=\frac{\bm{r}_{j}-\bm{u}_{j}}{\left|\bm{r}_{j}-\bm{u}_{j}\right|},$
$\displaystyle\bm{n}_{i}^{0}=\frac{\bm{r}_{i}-\bm{r}_{j}}{\left|\bm{r}_{i}-\bm{r}_{j}\right|},\bm{n}_{j}^{0}=-\bm{n}_{i}^{0},$
(61)

see Fig. 13. $R_{i}(\bm{n}_{i})$ denotes the radius of particle $i$ along the
direction $\bm{n}_{i}$. In particular,

$\displaystyle R_{i}(\bm{n}_{i}^{0})=R_{i}^{0}\left[1+\Delta
F(\theta_{i}-\theta_{ij})\right],$ (62)

where $\theta_{i}$ denotes the direction of particle $i$, $\theta_{ij}$
denotes the relative angle between particles $i$ and $j$, see Fig. 3. For
$\Delta\ll 1$, we can expand $h_{ij}$ w.r.t $\Delta$ as

$\displaystyle h_{ij}$
$\displaystyle=\left|\bm{u}_{i}-\bm{u}_{j}\right|=h_{ij}^{0}+\bm{n}_{i}^{0}\cdot\left(\delta\bm{u}_{i}-\delta\bm{u}_{j}\right)+O(\delta\bm{u}_{i}^{2},\delta\bm{u}_{j}^{2})$
$\displaystyle=h_{ij}^{0}+R_{i}^{0}-R_{i}(\bm{n}_{i})+R_{j}^{0}-R_{j}(\bm{n}_{j})+O(\Delta^{2})$
$\displaystyle=h_{ij}^{0}-\Delta\left[R_{i}^{0}F(\theta_{i}-\theta_{ij})+R_{j}^{0}F(\theta_{j}-\theta_{ji})\right]+O(\Delta^{2}),$
(63)

where we used $\bm{n}_{i}^{0}\cdot\bm{n}_{i}=1+O(\Delta^{2})$, and
$R_{i}(\bm{n}_{i})=R_{i}(\bm{n}_{i}^{0})+O(\Delta^{2})$. $h_{ij}^{0}$ denotes
the gap function of the reference disks:

$\displaystyle
h_{ij}^{0}=\left|\bm{u}_{i}^{0}-\bm{u}_{j}^{0}\right|=\left|\bm{r}_{i}-\bm{r}_{j}\right|-R_{i}^{0}-R_{j}^{0}.$
(64)

## References

  * van Hecke (2009) M. van Hecke, Journal of Physics: Condensed Matter 22, 033101 (2009).
  * Liu and Nagel (2010) A. J. Liu and S. R. Nagel, Annu. Rev. Condens. Matter Phys. 1, 347 (2010).
  * O’Hern  _et al._ (2003) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
  * Nishimori and Ortiz (2010) H. Nishimori and G. Ortiz,  _Elements of phase transitions and critical phenomena_ (OUP Oxford, 2010).
  * Vågberg  _et al._ (2011) D. Vågberg, D. Valdez-Balderas, M. Moore, P. Olsson, and S. Teitel, Physical Review E 83, 030303 (2011).
  * Charbonneau  _et al._ (2014) P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, and F. Zamponi, Nature communications 5, 1 (2014).
  * Ikeda (2020) H. Ikeda, Phys. Rev. Lett. 125, 038001 (2020).
  * Parisi  _et al._ (2020) G. Parisi, P. Urbani, and F. Zamponi,  _Theory of simple glasses: exact solutions in infinite dimensions_ (Cambridge University Press, 2020).
  * Mézard  _et al._ (1987) M. Mézard, G. Parisi, and M. Virasoro,  _Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications_ , Vol. 9 (World Scientific Publishing Company, 1987).
  * Nishimori (2001) H. Nishimori,  _Statistical physics of spin glasses and information processing: an introduction_ , 111 (Clarendon Press, 2001).
  * Wyart  _et al._ (2005) M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten, Physical Review E 72, 051306 (2005).
  * Yan  _et al._ (2016) L. Yan, E. DeGiuli, and M. Wyart, EPL (Europhysics Letters) 114, 26003 (2016).
  * DeGiuli  _et al._ (2014a) E. DeGiuli, A. Laversanne-Finot, G. Düring, E. Lerner, and M. Wyart, Soft Matter 10, 5628 (2014a).
  * DeGiuli  _et al._ (2014b) E. DeGiuli, E. Lerner, C. Brito, and M. Wyart, Proceedings of the National Academy of Sciences 111, 17054 (2014b).
  * Beltukov (2015) Y. Beltukov, JETP Letters 101, 345 (2015).
  * Ikeda and Shimada (2020) H. Ikeda and M. Shimada, arXiv preprint arXiv:2009.12060 (2020).
  * Brito  _et al._ (2018) C. Brito, H. Ikeda, P. Urbani, M. Wyart, and F. Zamponi, Proceedings of the National Academy of Sciences 115, 11736 (2018).
  * Ikeda  _et al._ (2019) H. Ikeda, P. Urbani, and F. Zamponi, Journal of Physics A: Mathematical and Theoretical 52, 344001 (2019).
  * Lu  _et al._ (2015) G. Lu, J. Third, and C. Müller, Chemical Engineering Science 127, 425 (2015).
  * Ikeda  _et al._ (2020a) H. Ikeda, C. Brito, and M. Wyart, Journal of Statistical Mechanics: Theory and Experiment 2020, 033302 (2020a).
  * Ikeda  _et al._ (2020b) H. Ikeda, C. Brito, M. Wyart, and F. Zamponi, Phys. Rev. Lett. 124, 208001 (2020b).
  * Goodrich  _et al._ (2014) C. P. Goodrich, S. Dagois-Bohy, B. P. Tighe, M. van Hecke, A. J. Liu, and S. R. Nagel, Phys. Rev. E 90, 022138 (2014).
  * Goodrich  _et al._ (2012) C. P. Goodrich, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. 109, 095704 (2012).
  * Bernal and Mason (1960) J. Bernal and J. Mason, Nature 188, 910 (1960).
  * Donev  _et al._ (2005) A. Donev, S. Torquato, and F. H. Stillinger, Phys. Rev. E 71, 011105 (2005).
  * Charbonneau  _et al._ (2012) P. Charbonneau, E. I. Corwin, G. Parisi, and F. Zamponi, Phys. Rev. Lett. 109, 205501 (2012).
  * Charbonneau  _et al._ (2020) P. Charbonneau, E. Corwin, C. Dennis, R. D. H. Rojas, H. Ikeda, G. Parisi, and F. Ricci-Tersenghi, arXiv preprint arXiv:2011.10899 (2020).
  * Franz  _et al._ (2017) S. Franz, G. Parisi, M. Sevelev, P. Urbani, and F. Zamponi, (2017).
  * Schreck  _et al._ (2010) C. F. Schreck, N. Xu, and C. S. O’Hern, Soft Matter 6, 2960 (2010).
  * Shiraishi  _et al._ (2019) K. Shiraishi, H. Mizuno, and A. Ikeda, Phys. Rev. E 100, 012606 (2019).
  * Shiraishi  _et al._ (2020) K. Shiraishi, H. Mizuno, and A. Ikeda, arXiv preprint arXiv:2005.02598 (2020).
  * Donev  _et al._ (2004) A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Stillinger, R. Connelly, S. Torquato, and P. M. Chaikin, Science 303, 990 (2004).
  * Donev  _et al._ (2007) A. Donev, R. Connelly, F. H. Stillinger, and S. Torquato, Physical Review E 75, 051304 (2007).
  * Zeravcic  _et al._ (2009) Z. Zeravcic, N. Xu, A. Liu, S. Nagel, and W. van Saarloos, EPL (Europhysics Letters) 87, 26001 (2009).
  * Mailman  _et al._ (2009) M. Mailman, C. F. Schreck, C. S. O’Hern, and B. Chakraborty, Physical review letters 102, 255501 (2009).
  * Schreck  _et al._ (2012) C. F. Schreck, M. Mailman, B. Chakraborty, and C. S. O’Hern, Physical Review E 85, 061305 (2012).
  * Williams and Philipse (2003) S. Williams and A. Philipse, Physical Review E 67, 051301 (2003).
  * Blouwolff and Fraden (2006) J. Blouwolff and S. Fraden, EPL (Europhysics Letters) 76, 1095 (2006).
  * Azéma and Radjaï (2010) E. Azéma and F. Radjaï, Phys. Rev. E 81, 051304 (2010).
  * Marschall and Teitel (2018) T. Marschall and S. Teitel, Physical Review E 97, 012905 (2018).
  * Jiao  _et al._ (2010) Y. Jiao, F. H. Stillinger, and S. Torquato, Phys. Rev. E 81, 041304 (2010).
  * Delaney and Cleary (2010) G. W. Delaney and P. W. Cleary, EPL (Europhysics Letters) 89, 34002 (2010).
  * VanderWerf  _et al._ (2018) K. VanderWerf, W. Jin, M. D. Shattuck, and C. S. O’Hern, Phys. Rev. E 97, 012909 (2018).
  * Tarama  _et al._ (2013) M. Tarama, A. M. Menzel, B. ten Hagen, R. Wittkowski, T. Ohta, and H. Löwen, The Journal of Chemical Physics 139, 104906 (2013).
  * Bitzek  _et al._ (2006) E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Physical review letters 97, 170201 (2006).
  * Shimada  _et al._ (2018) M. Shimada, H. Mizuno, M. Wyart, and A. Ikeda, Physical Review E 98, 060901 (2018).
  * Wyart (2005) M. Wyart, arXiv preprint cond-mat/0512155 (2005).