# Domain-Specific Pretraining of Language Models:
A Comparative Study in the Medical Field

Tobias Kerner
Technische Hochschule Ingolstadt
https://orcid.org/0009-0008-5485-1018

###### Abstract

There are many cases where LLMs are used for specific tasks in a single
domain. These usually require less general, but more domain-specific
knowledge. Highly capable, general-purpose state-of-the-art language models
like GPT-4 or Claude-3-opus can often be used for such tasks, but they are
very large and cannot be run locally, even if they were not proprietary. This
can be a problem when working with sensitive data. This paper focuses on
domain-specific and mixed-domain pretraining as potentially more efficient
methods than general pretraining for specialized language models. We will take
a look at work related to domain-specific pretraining, specifically in the
medical area, and compare benchmark results of specialized language models to
general-purpose language models.

###### Index Terms:

language model, domain-specific pretraining, medicine

##  I Introduction

Large Language Models (LLMs) have gained increasing popularity in recent
times. It is not unusual to see companies employ them as assistants on
websites or employ them directly in workflows like the development of internal
IT systems. They also have growing relevance in specialized fields like
medicine. LLMs could, for example, assist doctors with diagnoses or explain
complicated medical situations to patients in simple language. For models to
be used reliably, it is important for them to perform well in the domains of
their intended usage. High-performance general-purpose models currently
include GPT-4 and claude-3-opus, which can be applied to most tasks with
usually good results.

However, due to their large size and most state-of-the-art (SOTA) general-
purpose models being proprietary and only accessible through an API, it is
generally not possible to run these LLMs locally. This is a big problem when
reliability and privacy are the main concerns. When working with sensitive
patient data or with an unreliable internet connection, using an external API
is often not feasible.

A possible solution to this problem would be to use smaller, specialized LLMs
that perform well on tasks in a certain domain and whose performance in other
domains is irrelevant. Small models benefit from comparatively fast inference
times, lower latency, and less expensive training. Their reasonably small size
allows them to be hosted on local, consumer-grade hardware. Due to their lower
parameter count and limited learning capacity, small models inherently perform
worse than large models, given similar training conditions such as the amount
and quality of data used for training. Focusing the training on a specific
domain and ignoring capabilities in all other domains lets us circumvent this
limitation to some degree by not forcing the model to learn irrelevant or out-
of-domain information and allowing it to fully focus on memorizing domain-
related information.

But how can we create specialized models, and can they really be as good as
large general-purpose models in their domain? To answer these questions, we
will take a look at different methods of training LLMs, domain-specific
datasets and compare the benchmark results of specialized models to general-
purpose models in benchmarks related to their domain. This paper will focus on
the medical domain for all examples and comparisons.

##  II Pretraining

Pretraining is the most fundamental step in creating intelligent and capable
LLMs. During this training step, the model learns the structure of natural
language and tries to memorize as much of the training data as possible. For
causal language models, for example, this happens by inputting token sequences
into the model and letting it predict the next token in the sequence. The
closer the prediction is to the true next token in the sequence, the smaller
the loss. The loss is used to adjust the weights of the model, in order to
minimize the loss. When the model learns the contents of the data and improves
its predictions, the loss decreases. Pretraining teaches the model general
knowledge and language understanding, but not to answer in a chat-like way, as
many end-users are familiar with. For this usage, a further pretraining step
is often necessary, where the model is trained to reply in a certain format.
Usually, finetuning impacts the during pretraining learned knowledge only
slightly, since finetuning is a relatively short process compared to
pretraining and should not change the models weights a lot.

###  II-A General Pretraining

General pretraining is considered the default method for training large
language models and works well in most cases. There are many large datasets
that can be used for general pretraining without much modification. They
contain enough varied data to teach the model general language understanding
and general knowledge, even in multiple languages. During training, the model
will try to remember as much of the knowledge in the provided data as possible
in order to minimize its training loss. Large models may be able to learn most
of the data, but smaller models will struggle with this due to their inherent
parameter limit. This causes them to memorize all information in the data a
little bit, but nothing well. However, in the context of training a
specialized model, there are ways to circumvent this issue.

###  II-B Domain-Specific Pretraining

If the model will be used for tasks in a specific domain only, we can mostly
avoid this issue by simply removing everything from the data that is not
related to the target domain, leaving us with a domain-specific dataset. When
training a medical LLM for example, it is to be expected that a lack of
training data on poems or music, for example, should not affect its
performance on medical-related tasks. This is only possible if there is enough
data left in the domain-specific dataset for the model to get an understanding
of natural language from it while not overfitting [1].

The assumption is that the less irrelevant information the model has to learn,
the better it can learn the relevant information. This should make it possible
to decrease model size and accelerate training without a drastic decrease in
domain-related task performance. Even though a reduced parameter count will
decrease the model’s complexity and likely limit its ability to understand
very complex data structures, this should not be an issue for simpler tasks or
tasks that rely mostly on knowledge. An example of this in the medical field
would be the task of translating a detailed medical report into simple
language for the patient to understand, where the model has to know the
concept and meaning of medical terms in order to rephrase the report.

###  II-C Mixed-Domain Pretraining

Mixed-domain pretraining, also called continued pretraining, can be used if
there is not enough data in the domain-specific dataset for the model to gain
an understanding of natural language from it. We can first train the model on
a general-purpose dataset. This will enable the model to learn natural
language and, as mentioned before, a bit of all the different domains in the
data. This model with a general understanding of language can be used as a
base for continued pretraining on the domain-specific dataset. Since this
dataset will likely not include many of the things that the model has
memorized previously, training on it will result in the model gradually
forgetting most of them, which we can ignore when pretraining for domain-
specific tasks. If physical resources are limited and training an LLM from
scratch is not feasible, this method can also be applied using a third-party
general LLM as a base for continued pretraining.

continued pretraining on domain-specific data after training on general data
can lead to strong improvements of domain-specific models compared to
pretraining on general data only. It sets a stronger knowledge foundation for
in-domain finetuning. Their experiments show that small models struggle with
large amounts of knowledge and perform better on data specialized in a single
domain. [2]

However, if there is enough domain-specific data for domain-specific
pretraining without general pretraining before, performance on domain-specific
tasks may be better with domain-specific pretraining than mixed-domain
pretraining. While transfer-learning by pretraining on general data is useful
for models specialized on domains with low amounts of data, it could be
harmful for domain-specific models like medical models. This may be because
biomedical texts are substantially different from general texts or different
domains, so the transfer-learning might actually have a negative effect. [1]

This negative effect may arise from the embedding layer and the first few
layers in the model not benefiting from an early focus on domain-related data,
which could enable more efficient processing of complex domain-related data in
further layers.

##  III Datasets

There are many publicly available datasets that can be used for pretraining
LLMs, including general datasets as well as domain-specific ones. They are
created by collecting data from web scraping, books, articles, or other text-
based sources and then processing / cleaning that data. Their sizes can range
from a few kilobytes to terabytes.

###  III-A General Datasets

General-purpose datasets contain various types of content and tend to be a lot
larger than domain-specific ones. Table I shows some of the largest open-
source datasets currently available for pretraining. Their extensive size
allows the model to gain a solid understanding of language, as well as broad
general knowledge.

Dataset | Type | Size (GiB) | Language | Year |  Description
---|---|---|---|---|---
CommonCrawl[3] | Web | 391,168 | misc | 2024 |  2.7B webpages, June 2024 crawl
(first release in 2008)
The Stack[4] | Code | 3,136 | misc | 2022 |  Source codes in 30 different
programming languages
The Pile[5] | Mixed | 825 | en | 2020 |  Consisting of 22 sub-sets related
to various domains and tasks
RefinedWeb[6] | Web | 2,800 | en | 2023 |  Based on cleaned and refined data
from Common Crawl

TABLE I: Pretraining Datasets

It is to note that CommonCrawl-data is used partially in many other datasets,
such as RefinedWeb [6]. It contains an extensive collection of raw scraped
websites, which can be filtered, cleaned, and processed for the creation of
new, refined datasets. Figure 1 shows how the creators of RefinedWeb processed
CommonCrawl to create their dataset. A large part of the raw data gets
filtered out during language identification, during which they keep only
English texts, and removal of duplicates. [6]

Figure 1: Processing of CommonCrawl Dataset for RefinedWeb, taken from [6] and
edited for simplification.

###  III-B Medical Datasets Overview

There are publicly available domain-specific datasets as well. They can be
used to further refine a model’s knowledge and understanding in a specific
area. Table II shows select open-source medicine-related datasets, some of
which will also be used for benchmarking in the next section of this paper.
These are a lot smaller than the general datasets in table I due to the domain
limitation and the requirement for high-quality content. The quality
requirement for domain-specific datasets is higher since they tend to be used
to refine a base-model’s capabilities in a specific domain.

Some of the well-known medical datasets include the following: MeDAL consists
of PubMed-abstracts with around 3 abbreviations per abstract. The authors goal
was to create a dataset that teaches LLMs to handle medical-related
abbreviations better [7]. MedDialog contains English and Chinese dialogues.
The English dialogues have exactly 2 utterances per dialogue, with an average
length of 86.5 tokens per utterance. The Chinese dialogues are longer, with an
average of 3.3 utterances per dialogue up to a maximum of 198 utterances in
one dialogue and an average of 55.6 tokens per utterance [8]. MedQA consists
of English and Chinese question-answer pairs, collected from medical board
examinations in the USA, China, and Taiwan [9]. PubMedQA contains question-
answer pairs that are based on PubMed abstracts. It is split into 1k expert-
annotated, 61.2k unlabeled, and 211.3k artificially generated QA-pairs [10].
MedMCQA is designed to imitate real-world medical entrance exam questions. It
contains questions in natural language with four possible answers, one correct
answer, and an explanation for each [11].

Dataset | Type | Entries | Language | Release Year
---|---|---|---|---
MeDAL | Articles | 14.4m | en | 2020
MedDialog | Dialogues | 3.4m, 257.3k | en, ch | 2020
MedQA | QA | 12.7k, 48.4k | en, ch | 2021
PubMedQA | QA | 1.0k, 61.2k, 211.3k | en | 2019
MedMCQA | multiple-choice-QA | 193.1k | en | 2022

TABLE II: Overview and comparison of medical datasets

###  III-C Creating Domain-Specific Datasets

One of the biggest problems regarding training language models is the amount
of high-quality data available. For large-scale domain-specific training, we
will likely need more domain-related data than is available. For some domains,
there may not even be a publicly available dataset yet, or it does not fulfill
some requirements. In this case, we would have to create our own dataset.
There are several steps we have to take for creating our own domain-specific
dataset, including the acquisition of raw data, filtering and cleaning it, and
processing it into a format suitable for training.

For acquiring raw data, we can either collect it ourselves by means of
webscraping, using a dump of already scraped data such as CommonCrawl, or
collecting multiple datasets that are at least slightly related to our domain.
Large-scale webscraping is a time- and resource-intensive task, so it is
recommended to consider using data from CommonCrawl instead. The team behind
it releases regular updates to the dataset with newly scraped websites. Due to
its extensive size, it is more than likely to contain data related to our
domain. [3]

Assuming we chose to use CommonCrawl as the source for our raw data, we can
now filter it for data related to our domain, which would be medicine in our
case. One way to do this is by iterating through the scraped websites and
looking for medicine-related keywords, saving the websites containing them
separately. Other methods to filter for relevant data include Named Entity
Recognition, topic modeling (using techniques such as Latent Dirichlet
Analysis) or using a pretrained language model to tell us if a text is of
medical nature or not. These methods can yield better results than simple
keyword-matching, but are more complex. The best filtering can likely be
achieved by using a combination multiple methods.

After collecting enough in-domain data, it will have to be cleaned. This means
removing duplicate sites, as well as noise and irrelevant data embedded inside
the relevant text. Noise can, for example, be caused by websites inserting
advertisements in the middle of their articles. Keeping this would be harmful
to our training since it would likely distract or confuse the model and hence
has to be removed.

This cleaned data containing medical text could be used for training now, but
we can further refine it by creating question-answer pairs for it. This can be
especially useful for teaching causality and the relationship between concepts
and entities in the text. For this, there are tools like Augmentoolkit.
Provided with data in text format, it splits the data into chunks and employs
an LLM to extract suitable questions and the corresponding answers from these
chunks [12]. Due to the nature of language models, it is not guaranteed that
all extracted question-answer pairs are correct, even assuming that the
provided data was flawless. But this disadvantage would likely be negligible
in comparison to the potential benefit of training on data in this format. Of
course, we can first train our LLM on the normal cleaned data and then use
this refined question-answer data later for finetuning.

##  IV Performance of specialized LLMs

As mentioned previously, specialized LLMs are usually smaller than general-
purpose LLMs. They are trained on less out-of-domain data, speeding up
training time and reducing compute cost. In the following, we will take a look
at medical LLMs trained with domain-specific pretraining and mixed-domain
pretraining and compare their performance to general-purpose LLMs of equal or
larger size in medical-related tasks.

###  IV-A Domain-Specific Pretrained Models

A well-known medicine-focused model is PubMedBERT. It is based on the BERT
architecture [13] and trained on 14 million PubMed-abstracts with a total
length of 3.1 billion words. The model was trained from scratch, only using
the BERT architecture but not the pretrained weights of the original BERT
model. Its tokenizer uses vocabulary newly generated with the WordPiece
algorithm, based on the PubMed-data that was also used for training. [1]

BioMedLM is a 2.7B model trained exclusively on biomedical text, including
PubMed-abstracts and full articles. When finetuned to reply in the format
required by a specific task, BioMedLM is competitive with much larger models,
as shown by its score of 57.3% on MedMCQA (dev) and 69.0% on the MMLU Medical
Genetics exam. Even though BioMedLM only has 2.7B parameters, it beats GPT-3.5
with 175B parameters on MedMCQA and is on par with it in the MedQA and
PubMedQA benchmarks. However, it seems to struggle against larger models in
MMLU general medical knowledge benchmarks. [14]

Apollo is a family of small LLMs ranging from 0.5B to 7B parameters in size.
It is trained using a custom multilingual medical dataset called
ApolloCorpora. This dataset is based on books, papers, and encyclopedias
related to the medical domain, as well as doctor-patient dialogues with a
total of 2.5B training tokens. [15]

###  IV-B Mixed-Domain Pretrained

HEAL is based on LLaMA2-13B, applying continued pretraining to the weights of
the original model with 14.89B tokens of medical data. It surpasses GPT-4 on
PubMedQA, with a score closely behind Med-PaLM-2. It seems to struggle on the
MedQA benchmark, with only a small improvement over the base LLaMA2-13B and
quite behind another LLaMA-13B-based model. It is worth noting that its score
of 47.2 on MedQA is around the same level as GPT-3.5. [16, 17]

###  IV-C Benchmark Comparison

Figure 2: Comparison of models on medical benchmarks, sorted by model size.
The number in brackets behind the score represents the evaluation method
(x-shot). f stands for finetune, the model was finetuned for this task. if the
reference paper does not provide the evaluation method or specifies ’few-
shot’, there brackets will contain a ’?’.

Model | Size | MedMCQA | PubMedQA | MedQA | MMLU
---|---|---|---|---|---
Apollo | 0.5B | 37.8 (3) | - | - | 45.9 (3)
Apollo | 1.8B | 45.0 (3) | - | - | 49.1 (3)
BioMedLM | 2.7B | 57.3 (f) | 74.4 (f) | 50.3 (f) | 63.1 (f)
LLaMA 2 | 7B | 36.6 (3) | - | - | 40.1 (3)
Apollo | 7B | 58.2 (3) | - | - | 71.9 (3)
LLaMA 2 | 13B | 37.4 (f) | 76.4 (3) | 45.5 (3) | -
HEAL | 13B | - | 78.4 (3) | 47.2 (3) | -
LLaMA 2 | 70B | 48.3 (3) | - | - | 64.6 (3)
GPT-3.5 | 175B | 51.0 (0) | 71.6 (0) | 44.6 (0) | 70.6 (0)
GPT-4 |  1.7T1 | 73.4 (0) | 80.4 (0) | 81.4 (0) | 92.6 (0)

  * 1

rumoured parameter count [18]

As shown in table 2, the specialized LLMs perform well in in-domain tasks with
respect to their size when compared to not domain-specifically trained models.
Especially BioMedLM shows strong performance, beating GPT-3.5 in MedMCQA,
PubMedQA, and MedQA with only 2.7B parameters. The Apollo models perform good
as well, with Apollo-7B improving upon BioMedLM’s MMLU score, slightly beating
GPT-3.5. This is surprising, considering it was trained on 2.5B tokens of data
only, compared to BioMedLM’s 300B tokens. With a size of just 0.5B parameters,
Apollo-0.5B reaches solid scores in MedMCQA and MMLU, beating the about 14
times larger LLaMA-2-7B. Interestingly, both GPT-3.5 and GPT-4 perform worse
on PubMedQA with 5-shot than with 0-shot, having a score of 60.2 and 77.4 with
5-shot, respectively [17]. ††Data for table 2: The parts of MMLU used for
calculating the score for all Apollo models, LLaMA-2-70B, and LLaMA-2-7B were
specified as ’medical-related parts of MMLU’ [15]. For BioMedLM, GPT-3.5, and
GPT-4, the MMLU score is the average of the parts ’Clinical Knowledge’,
’Professional Medicine’, ’College Biology’, and ’Medical Genetics’. The scores
of GPT-4 are for the base-version of GPT-4, not to be confused with the
publicly available version of GPT-4 that was finetuned for more safe and
responsible answers and generally performs worse on medical benchmarks. The
MedQA score for GPT-3.5 and GPT-4 is based on ’US (5-option)’. [17] The
LLaMA-2-13B MedMCQA score was taken from [19].

$10^{0}$$10^{1}$$10^{2}$$10^{3}$$10^{4}$$10^{5}$$65$$70$$75$$80$$85$$90$$95$$100$GPT-4GPT-3.5Llama-2-13BBioMedLMHEALmodel
parameter count (B)PubMedQA score (inverted)General pretrainedDomain
specificMixed-domain

$10^{0}$$10^{1}$$10^{2}$$10^{3}$$10^{4}$$40$$60$$80$$100$GPT-4GPT-3.5Llama-2-70BLlama-2-13BLlama-2-7BBioMedLMApollo-7BApollo-2.7BApollo-0.5Bmodel
parameter count (B)MedMCQA score (inverted)General pretrainedDomain
specificMixed-domain

Figure 3: Comparison of model scores on PubMedQA and MedMCQA. Circle size is
logarithmically proportional to the number of tokens the model was trained
with. Number of training tokens for general models taken from [20, 18, 21].

###  IV-D Comparing Model-Size to Benchmark Score

$10^{0}$$10^{1}$$10^{2}$$10^{3}$$10^{4}$$40$$60$$80$$100$GPT-4GPT-3.5Llama-2-70BLlama-2-13BLlama-2-7BBioMedLMApollo-7BApollo-2.7BApollo-0.5Bmodel
parameter count (B)MedMCQA score (inverted)General pretrainedDomain
specificMixed-domain Figure 4: Comparison of model scores on MedMCQA with a
trend line based on general pretrained model data. Circle size is
logarithmically proportional to the number of tokens the model was trained
with. Trend line function: 9.226*ln(24.64*x+665.24)-25.59

Given a high-quality dataset, a model’s performance seems to be directly
related to its parameter count. According to our visualization in table 3,
especially the MedMCQA benchmark, we can see a strong relationship between
model size and benchmark score for the general pretrained LLMs. In 4, we can
see the MedMCQA benchmark plot with a line based on the general pretrained
model scores. The domain-specific and mixed-domain pretrained models seem to
be able to ’escape’ this line and achieve a better score compared to their
size. Based on this figure, BioMedLM has the best performance-to-parameter-
count relationship. This is likely possible due to the large quantity of
medical data it was trained on, enabling it to saturate its 2.7B weights with
medical knowledge. It is possible that this relationship is not as strong in
different benchmarks, but this example is an indication that it might be at
least similar.

###  IV-E Further Resource-Optimization

In the context of our initial question, we asked if LLMs small enough to be
run locally on consumer-grade hardware can match the performance of large,
general-purpose models that cannot be run locally in certain tasks. We should
be able to run a model with 7B parameters on local hardware with around 16GB
of VRAM, assuming the weights are saved as (b)float16. But a further reduction
in memory requirements and increase in inference speed are possible by
quantization of the model’s weights.

Current popular quantization formats include gguf, exl2, awq, and more. Gguf
allows for some of the model’s layers to be loaded into cpu-memory if the
model does not fully fit into gpu-memory, at the cost of inference speed due
to slower inference on the cpu. While the original Apollo-7B with (b)float16
weights would require slightly less than 16GB of gpu-memory due to each of the
slightly more than 7 billion weights each taking up 16 bits / 2 bytes, the
Q8_K quantized version only takes up 9.2GB. Most commonly used quantizations
go down to Q4, which stores each weight with an average of 4 bits. The fewer
bits are used to store weights, the stronger the degradation of the model’s
performance becomes. Below 4-bit quantization, especially for smaller models
(7B and below), the precision is often too low to be of reliable use. The
Q4_K_S and Q4_K_M quantized versions of Apollo-7B take up only 5.1 and 5.4 GB
of memory, respectively. The difference between the K_S and the K_M versions
is that the K_M version employs a less strict quantization, with some more
performance-impacting layers being kept at higher quanitization levels,
resulting in a slightly higher memory requirement. Some weights are not
quantized at all in the gguf model, such as those of normalization layers.
Because of this, we do not see a perfect 50% reduction in memory requirements
when switching from the 16-bit model to the 8-bit model. Although the gguf-
quantized models huggingface-page does not specifically state the precision
changes of the model in the quantized versions, the 8-bit gguf quantization is
generally considered to have a precision very similar to the original model,
making it a viable and, in some cases, even better-suited option than the
16-bit model. [22]

##  V Summary

Depending on the use case and available resources, domain-specific and mixed-
domain pretraining can be a viable and preferable alternative to general
pretraining. With a low amount of in-domain data available, mixed-domain
pretraining on a large general dataset and then switching to the small in-
domain dataset will likely result in a better-performing model than training
on a small domain-specific dataset only. These specialized models can have a
much lower parameter count but still perform impressively well on in-domain
tasks compared to larger, general pretrained models. This makes it possible to
run them on consumer-grade hardware, enabling independent local inference.
Therefore, domain-specific pretrained LLMs can be a good alternative to api-
based SOTA general-purpose models in some cases, such as supporting symptom
analysis and rephrasing documents in a medical setting.

## References

  * [1] Y. Gu, R. Tinn, H. Cheng, M. Lucas, N. Usuyama, X. Liu, T. Naumann, J. Gao, and H. Poon, “Domain-specific language model pretraining for biomedical natural language processing,” _ACM Trans. Comput. Healthcare_ , vol. 3, no. 1, oct 2021. [Online]. Available: https://doi.org/10.1145/3458754
  * [2] S. Gururangan, A. Marasović, S. Swayamdipta, K. Lo, I. Beltagy, D. Downey, and N. A. Smith, “Don’t stop pretraining: Adapt language models to domains and tasks,” in _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , D. Jurafsky, J. Chai, N. Schluter, and J. Tetreault, Eds. Online: Association for Computational Linguistics, Jul. 2020, pp. 8342–8360. [Online]. Available: https://aclanthology.org/2020.acl-main.740
  * [3] “Common Crawl - commoncrawl.org,” 2024, [Accessed 14-07-2024]. [Online]. Available: https://commoncrawl.org
  * [4] L. Gao, S. Biderman, S. Black, L. Golding, T. Hoppe, C. Foster, J. Phang, H. He, A. Thite, N. Nabeshima, S. Presser, and C. Leahy, “The pile: An 800gb dataset of diverse text for language modeling,” 2020. [Online]. Available: https://arxiv.org/abs/2101.00027
  * [5] D. Kocetkov, R. Li, L. B. Allal, J. Li, C. Mou, C. M. Ferrandis, Y. Jernite, M. Mitchell, S. Hughes, T. Wolf, D. Bahdanau, L. von Werra, and H. de Vries, “The stack: 3 tb of permissively licensed source code,” 2022. [Online]. Available: https://arxiv.org/abs/2211.15533
  * [6] G. Penedo, Q. Malartic, D. Hesslow, R. Cojocaru, H. Alobeidli, A. Cappelli, B. Pannier, E. Almazrouei, and J. Launay, “The refinedweb dataset for falcon llm: Outperforming curated corpora with web data only,” in _Advances in Neural Information Processing Systems_ , A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, Eds., vol. 36. Curran Associates, Inc., 2023, pp. 79 155–79 172.
  * [7] Z. Wen, X. H. Lu, and S. Reddy, “MeDAL: Medical abbreviation disambiguation dataset for natural language understanding pretraining,” in _Proceedings of the 3rd Clinical Natural Language Processing Workshop_ , A. Rumshisky, K. Roberts, S. Bethard, and T. Naumann, Eds. Online: Association for Computational Linguistics, Nov. 2020, pp. 130–135. [Online]. Available: https://aclanthology.org/2020.clinicalnlp-1.15
  * [8] G. Zeng, W. Yang, Z. Ju, Y. Yang, S. Wang, R. Zhang, M. Zhou, J. Zeng, X. Dong, R. Zhang, H. Fang, P. Zhu, S. Chen, and P. Xie, “MedDialog: Large-scale medical dialogue datasets,” in _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , B. Webber, T. Cohn, Y. He, and Y. Liu, Eds. Online: Association for Computational Linguistics, Nov. 2020, pp. 9241–9250. [Online]. Available: https://aclanthology.org/2020.emnlp-main.743
  * [9] D. Jin, E. Pan, N. Oufattole, W.-H. Weng, H. Fang, and P. Szolovits, “What disease does this patient have? a large-scale open domain question answering dataset from medical exams,” _Applied Sciences_ , vol. 11, no. 14, 2021. [Online]. Available: https://www.mdpi.com/2076-3417/11/14/6421
  * [10] Q. Jin, B. Dhingra, Z. Liu, W. Cohen, and X. Lu, “PubMedQA: A dataset for biomedical research question answering,” in _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , K. Inui, J. Jiang, V. Ng, and X. Wan, Eds. Hong Kong, China: Association for Computational Linguistics, Nov. 2019, pp. 2567–2577. [Online]. Available: https://aclanthology.org/D19-1259
  * [11] A. Pal, L. K. Umapathi, and M. Sankarasubbu, “Medmcqa: A large-scale multi-subject multi-choice dataset for medical domain question answering,” in _Proceedings of the Conference on Health, Inference, and Learning_ , ser. Proceedings of Machine Learning Research, G. Flores, G. H. Chen, T. Pollard, J. C. Ho, and T. Naumann, Eds., vol. 174. PMLR, 07–08 Apr 2022, pp. 248–260. [Online]. Available: https://proceedings.mlr.press/v174/pal22a.html
  * [12] E. Armstrong, cocktailpeanut, darkacorn, A. T. (afterSt0rm), abasgames, juanjopc, and RyanGreenup, “e-p-armstrong/augmentoolkit: Classifier Creator & General Overhaul,” Jul. 2024. [Online]. Available: https://doi.org/10.5281/zenodo.12702005
  * [13] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” 2019. [Online]. Available: https://arxiv.org/abs/1810.04805
  * [14] E. Bolton, A. Venigalla, M. Yasunaga, D. Hall, B. Xiong, T. Lee, R. Daneshjou, J. Frankle, P. Liang, M. Carbin, and C. D. Manning, “Biomedlm: A 2.7b parameter language model trained on biomedical text,” 2024.
  * [15] X. Wang, N. Chen, J. Chen, Y. Hu, Y. Wang, X. Wu, A. Gao, X. Wan, H. Li, and B. Wang, “Apollo: An lightweight multilingual medical llm towards democratizing medical ai to 6b people,” 2024.
  * [16] D. Yuan, E. Rastogi, G. Naik, S. P. Rajagopal, S. Goyal, F. Zhao, B. Chintagunta, and J. Ward, “A continued pretrained LLM approach for automatic medical note generation,” in _Proceedings of the 2024 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies (Volume 2: Short Papers)_ , K. Duh, H. Gomez, and S. Bethard, Eds. Mexico City, Mexico: Association for Computational Linguistics, Jun. 2024, pp. 565–571. [Online]. Available: https://aclanthology.org/2024.naacl-short.47
  * [17] H. Nori, N. King, S. M. McKinney, D. Carignan, and E. Horvitz, “Capabilities of gpt-4 on medical challenge problems,” 2023.
  * [18] P. McGuinness, “GPT-4 Details Revealed — patmcguinness.substack.com,” 2023, [Accessed 15-07-2024]. [Online]. Available: https://patmcguinness.substack.com/p/gpt-4-details-revealed
  * [19] C. Wu, W. Lin, X. Zhang, Y. Zhang, W. Xie, and Y. Wang, “PMC-LLaMA: toward building open-source language models for medicine,” _Journal of the American Medical Informatics Association_ , p. ocae045, 04 2024. [Online]. Available: https://doi.org/10.1093/jamia/ocae045
  * [20] T. B. Brown, B. Mann, N. Ryder, M. Subbiah, J. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, S. Agarwal, A. Herbert-Voss, G. Krueger, T. Henighan, R. Child, A. Ramesh, D. M. Ziegler, J. Wu, C. Winter, C. Hesse, M. Chen, E. Sigler, M. Litwin, S. Gray, B. Chess, J. Clark, C. Berner, S. McCandlish, A. Radford, I. Sutskever, and D. Amodei, “Language models are few-shot learners,” 2020. [Online]. Available: https://arxiv.org/abs/2005.14165
  * [21] H. Touvron, L. Martin, K. Stone, P. Albert, A. Almahairi, Y. Babaei, N. Bashlykov, S. Batra, P. Bhargava, S. Bhosale, D. Bikel, L. Blecher, C. C. Ferrer, M. Chen, G. Cucurull, D. Esiobu, J. Fernandes, J. Fu, W. Fu, B. Fuller, C. Gao, V. Goswami, N. Goyal, A. Hartshorn, S. Hosseini, R. Hou, H. Inan, M. Kardas, V. Kerkez, M. Khabsa, I. Kloumann, A. Korenev, P. S. Koura, M.-A. Lachaux, T. Lavril, J. Lee, D. Liskovich, Y. Lu, Y. Mao, X. Martinet, T. Mihaylov, P. Mishra, I. Molybog, Y. Nie, A. Poulton, J. Reizenstein, R. Rungta, K. Saladi, A. Schelten, R. Silva, E. M. Smith, R. Subramanian, X. E. Tan, B. Tang, R. Taylor, A. Williams, J. X. Kuan, P. Xu, Z. Yan, I. Zarov, Y. Zhang, A. Fan, M. Kambadur, S. Narang, A. Rodriguez, R. Stojnic, S. Edunov, and T. Scialom, “Llama 2: Open foundation and fine-tuned chat models,” 2023. [Online]. Available: https://arxiv.org/abs/2307.09288
  * [22] “FreedomIntelligence/Apollo-7B-GGUF - huggingface.co,” 2024, [Accessed 16-07-2024]. [Online]. Available: https://huggingface.co/FreedomIntelligence/Apollo-7B-GGUF


#  Softening and residual loss modulus of jammed grains under oscillatory
shear in an absorbing state

Michio Otsuki<EMAIL_ADDRESS>Graduate School of Engineering
Science, Osaka University, Toyonaka, Osaka 560-8531, Japan  Hisao Hayakawa
Yukawa Institute for Theoretical Physics, Kyoto University,
Kitashirakawaoiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan

###### Abstract

From a theoretical study of the mechanical response of jammed materials
comprising frictionless and overdamped particles under oscillatory shear, we
find that the material becomes soft, and the loss modulus remains non-zero
even in an absorbing state where any irreversible plastic deformation does not
exist. The trajectories of the particles in this region exhibit hysteresis
loops. We succeed in clarifying the origin of the softening of the material
and the residual loss modulus with the aid of Fourier analysis. We also
clarify the roles of the yielding point in the softening to distinguish the
plastic deformation from reversible deformation in the absorbing state.

Introduction— The mechanical response of jammed disordered materials, such as
granular materials, foams, emulsions, and colloidal suspensions, garners much
attention Hecke ; Behringer . For vanishingly small strain, the shear stress
$\sigma$ is proportional to the shear strain $\gamma$, which is characterized
by the shear modulus satisfying a critical scaling law near the jamming point
$\phi_{J}$ OHern02 ; Tighe11 ; Otsuki17 . However, the region of the linear
response is quite narrow near $\phi_{J}$ Coulais ; Otsuki14 . Hence, revealing
the nonlinear response is essential for understanding the dynamics of
disordered materials.

In crystalline materials, the nonlinear response originates from yielding
associated with irreversible plastic deformation. Yielding also takes place in
disordered materials when the strain is sufficiently large Nagamanasa ;
Knowlton ; Kawasaki16 ; Leishangthem ; Clark ; Boschan19 . The yielding
transition attracts much attention among researchers as an example of the
reversible-irreversible transition Hinrichsen ; Henkel ; Pine ; Corte . When
plastic deformation causes rearrangements of contact networks, the mechanical
response becomes nonlinear. It had been believed that plastic deformation is
always necessary for the nonlinear response. Unlike this expectation, recent
studies have revealed that plastic deformation is not always necessary for the
nonlinear response Boschan ; Nakayama ; Kawasaki20 ; Bohy ; Ishima . Under
steady shear, $\sigma$ becomes hypoelastic before the yielding Boschan ;
Kawasaki20 , and the storage modulus in the steady state after applying a
sufficient number of cyclic shears decreases as the strain amplitude increases
without any irreversible plastic deformation Bohy . The decrease of the
storage modulus is called softening.

It is known that plastic deformation causes dissipation characterized by the
loss modulus Bohy ; Ishima . It is natural that the loss modulus disappears in
quasi-static strains without any plastic deformation. However, we need careful
check of this naive picture, because the loss modulus might be related to the
softening observed without any plastic deformation.

The mechanical response should be related to the motion of particles
constituting the disordered materials. This suggests that the trajectories of
particles provide information on the softening of the materials. Several
studies have reported that the trajectories of dense particles form closed
loops under oscillatory shear below the yielding point associated with
reversible contact changes where there are some cyclic open and closed
contacts between particles Lundberg ; Schreck ; Keim13 ; Keim14 ; Regev13 ;
Regev15 ; Priezjev ; Lavrentovich ; Nagasawa ; Das ; Deen ; Khirallah . The
formation of closed loops means that the system is reduced to an absorbing
state after some time has passed. A previous study numerically showed that the
softening in the absorbing state becomes significant when there are closed
loops associated with many contact changes. However, the quantitative
relationship remains unclear Bohy .

In this study, we numerically investigate jammed materials comprising $N$
frictionless and overdamped particles under oscillatory shear to clarify the
origin of the softening. For this purpose, we focus on the roles of the
trajectories to clarify the relationship between the softening in the
absorbing state and the softening in the plastic regime. We find that the
shear modulus exhibits softening, and the loss modulus remains non-zero even
in the absorbing state below the yielding point. The trajectory of a test
particle forms a nontrivial loop in this region. With the aid of Fourier
analysis, we investigate the geometric structure of the trajectories and
reveal the role of Fourier components for the storage and loss moduli. We also
present the theoretical expressions for the storage and loss moduli, whose
quantitative validities are numerically confirmed.

Setup— Let us consider a jammed two-dimensional system consisting of
frictionless particles under oscillatory shear. The particles are driven by
the overdamped equation with Stokes’ drag under Lees–Edwards boundary
conditions Evans , where the equation of motion is given by

$\zeta\left\\{\frac{d}{dt}{\bm{r}}_{i}-\dot{\gamma}(t)y_{i}\bm{e}_{x}\right\\}=-\sum_{j\neq
i}\frac{\partial}{\partial\bm{r}_{i}}U(r_{ij})$ (1)

with the position ${\bm{r}}_{i}=(x_{i},y_{i})$ of particle $i$. Here, $\zeta$
and $\dot{\gamma}(t)$ are the drag coefficient and strain rate, respectively.
The interaction potential $U(r_{ij})$ is assumed to be

$U(r_{ij})=\frac{k}{2}(d_{ij}-r_{ij})^{2}\Theta(d_{ij}-r_{ij}),$ (2)

where $\Theta(x)$, $k$, $d_{ij}$, and
$r_{ij}=|\bm{r}_{ij}|=|\bm{r}_{i}-\bm{r}_{i}|$ are the Heaviside step function
satisfying $\Theta(x)=1$ for $x\geq 0$ and $\Theta(x)=0$ otherwise, the spring
constant, the average diameter of particles $i$ and $j$, and the distance
between particles $i$ and $j$, respectively. The system is bidisperse and
consists of an equal number of particles with diameters $d_{0}$ and
$d_{0}/1.4$. We have verified that particles with inertia and damping at
contact, which corresponds to the model in Ref. Bohy , exhibit almost
identical behavior in our system Supple .

We prepare the initial state with a given packing fraction $\phi$ by slowly
compressing the system from a state below the jamming point $\phi_{\rm
J}\simeq 0.841$ Otsuki17 . The oscillatory shear strain is applied for $n_{c}$
cycles as

$\gamma(\theta)=\gamma_{0}\sin\theta$ (3)

with the phase $\theta=\omega t$, where $\gamma_{0}$ and $\omega$ are the
strain amplitude and angular frequency, respectively. Note that the shear rate
satisfies $\dot{\gamma}(t)=(d\theta/dt)(d/d\theta)\gamma(\theta)$. In the last
cycle, we measure the storage and loss moduli $G^{\prime}$ and
$G^{\prime\prime}$, respectively, given by Doi

$\displaystyle G^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int_{0}^{2\pi}\ d\theta\
\frac{\left\langle\sigma(\theta)\right\rangle\sin\theta}{\gamma_{0}},$ (4)
$\displaystyle G^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int_{0}^{2\pi}\ d\theta\
\frac{\left\langle\sigma(\theta)\right\rangle\cos\theta}{\gamma_{0}},$ (5)

with shear stress

$\sigma=\frac{1}{L^{2}}\sum_{i}\sum_{j>i}\frac{x_{ij}y_{ij}}{r_{ij}}U^{\prime}(r_{ij}),$
(6)

where $x_{ij}=x_{i}-x_{j}$, $y_{ij}=y_{i}-y_{j}$, $\langle\cdot\rangle$
represents the ensemble average, and $L$ is the linear system size. See Ref.
Supple  for the stress-strain curves in our system. We have verified that
$G^{\prime}$ and $G^{\prime\prime}$ are independent of $N$ and $n_{c}$ for
$N\geq 1000$ and $n_{c}\geq 20$ Supple . We use $N=1000$ and $n_{c}=20$ in our
numerical analysis. We adopt the Euler method using the time step $\Delta
t=0.05\tau_{0}$ with $\tau_{0}=\zeta/k$.

Closed Trajectories— As the number of cycles increases, the system reaches a
statistically steady state through a transient regime as shown in Ref. Supple
. Figure 1 displays typical non-affine trajectories of a particle

$\tilde{\bm{r}}_{i}(\theta)=\bm{r}_{i}(\theta)-\gamma(\theta)y_{i}(\theta)\bm{e}_{x}$
(7)

in the last two cycles with $\phi=0.87$ and $\omega=10^{-4}\tau_{0}^{-1}$ in
the steady state. In Fig.1(a) ($\gamma_{0}=0.02$), the trajectories are
closed, and the particle returns to its original position after every cycle.
This indicates that irreversible plastic deformation does not occur, at least
in the last two cycles. The closed trajectories form nontrivial loops, which
differ from ellipses or lines observed for small $\gamma_{0}$ as shown in Ref.
Supple . In Fig. 1(b) ($\gamma_{0}=0.1$), the particle moves away from its
original positions after a cycle, as a characteristic behavior of plastic
deformation. Here, we define the absorbing state where the displacement of
each particle after several cycles is smaller than $d_{c}=10^{-4}d_{0}$ in the
statistically steady state. We also define the plastic state where the
displacement after several cycles exceeds $d_{c}$. It should be noted that
some rare samples exhibit trajectories where particles return to their
original positions after more than one cycle Regev13 ; Regev15 ; Lavrentovich
; Nagasawa ; Khirallah . However, our theoretical results shown below are
unchanged even if such samples exist Supple .

---
Figure 1:  Non-affine particle trajectories in the last two cycles for
$\gamma_{0}=0.02$ (a) and $0.1$ (b) with $\omega=10^{-4}\tau_{0}^{-1}$ and
$\phi=0.87$, which corresponds to $\phi-\phi_{J}=0.029$. The circles represent
the trajectory in the last cycle. The line represents the trajectory in the
second to the last cycle.

Shear Modulus— We plot the storage modulus $G^{\prime}$ against the strain
amplitude $\gamma_{0}$ for $\omega=10^{-4}\tau_{0}^{-1}$ with
$\phi=0.870,0.860,0.850,$ and $0.845$ in Fig. 2. The yielding points to
distinguish the absorbing state from the plastic state for various $\phi$ are
shown by open pentagons Supple . The storage modulus $G^{\prime}$ decreases as
$\gamma_{0}$ increases, but the yielding point is not identical to the point
where $G^{\prime}$ starts to decrease. We call the decrease for
$\gamma_{0}<\gamma_{c}$, the yielding strain amplitude, the softening in the
absorbing state (SAS). We also call the decrease for $\gamma_{0}>\gamma_{c}$
the softening in the plastic state (SPS). It is remarkable that SAS is
continuously connected to SPS, while a shoulder in $G^{\prime}$ appears in SPS
for $0.04\leq\gamma_{0}\leq 0.1$ with $\phi=0.845$. In the inset of Fig. 2, we
demonstrate that $G^{\prime}$ and $\gamma_{0}$ can be scaled by
$\sqrt{\phi-\phi_{J}}$ and $\phi-\phi_{J}$, respectively, as indicated in
Refs. OHern02 ; Bohy . We have confirmed that $G^{\prime}$ is independent of
$\omega$ for $\omega\leq 10^{-3}\tau_{0}^{-1}$.

Figure 2:  Storage modulus $G^{\prime}$ obtained in our simulation (filled
symbols) against $\gamma_{0}$ for $\omega=10^{-4}\tau_{0}^{-1}$ with
$\phi=0.870,0.860,0.850,$ and $0.845$, which corresponds to
$\phi-\phi_{J}=0.029,0.019,0.009$, and $0.004$, respectively. The legends
represent the packing fraction $\phi$. The data in the absorbing (plastic)
state obtained in our simulation are shown in larger (smaller) filled symbols.
The open pentagons represent the yielding strain amplitude $\gamma_{c}$, while
other open symbols represent the theoretical expression using $G^{\prime}_{\rm
T}$ in Eq. (14). (Inset) Scaled storage modulus
$\tilde{G}^{\prime}=G^{\prime}/\sqrt{\phi-\phi_{J}}$ obtained in our
simulation (filled symbols) and its theoretical expression using
$G^{\prime}_{\rm T}$ (open symbols) in Eq. (14) against scaled strain
amplitude $\tilde{\gamma}_{0}=\gamma_{0}/(\phi-\phi_{J})$ in the absorbing
state.

Figure 3(a) displays the loss modulus $G^{\prime\prime}$ in the absorbing
state against $\gamma_{0}$ for $\omega=10^{-4}\tau_{0}^{-1}$ with
$\phi=0.870,0.860,0.850,$ and $0.845$, in which $G^{\prime\prime}$ does not
strongly depend on $\phi$ and $\gamma_{0}$. See Ref. Supple  for
$G^{\prime\prime}$ in the plastic state. In Fig. 3(b), we plot the loss
modulus $G^{\prime\prime}$ in the absorbing state against $\omega$ for
$\phi=0.87$ with $\gamma=0.01$. Remarkably, $G^{\prime\prime}$ in Fig. 3(b)
seems to converge to a non-zero value in the limit $\omega\to 0$, which
contrasts with the behavior of the Kelvin–Voigt model (i.e.,
$G^{\prime\prime}\propto\omega$ Meyers ). This behavior indicates that
dissipation remains even in the quasi-static limit in the absorbing state.
Note that $G^{\prime\prime}\propto\omega$ is recovered when we adopt a
sufficiently small $\gamma_{0}$ Supple .

Figure 3: (a) Loss modulus $G^{\prime\prime}$ in the absorbing state obtained
in our simulation (filled symbols) and its theoretical expression
$G^{\prime\prime}_{\rm T}$ (open symbols) in Eq. (15) against $\gamma_{0}$ for
$\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870,0.860,0.850,$ and $0.845$,
which corresponds to $\phi-\phi_{J}=0.029,0.019,0.009$, and $0.004$,
respectively. (b) Loss modulus $G^{\prime\prime}$ against $\omega\tau_{0}$ for
$\phi=0.87$ with $\gamma_{0}=0.01$.

Fourier Analysis— In the absorbing state, the non-affine trajectory
$\tilde{\bm{r}}_{i}(\theta)$ of particle $i$ can be expressed in a Fourier
series as

$\tilde{\bm{r}}_{i}(\theta)=\bm{R}_{i}+\sum_{n=1}^{\infty}\left(\bm{a}_{i}^{(n)}\sin
n\theta+\bm{b}_{i}^{(n)}\cos n\theta\right)$ (8)

with the center of the trajectory

$\displaystyle\bm{R}_{i}$ $\displaystyle=$
$\displaystyle(X_{i},Y_{i})=\frac{1}{2\pi}\int_{0}^{2\pi}\ d\theta\
\tilde{\bm{r}}_{i}(\theta),$ (9)

and the Fourier coefficients

$\displaystyle\bm{a}_{i}^{(n)}$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int_{0}^{2\pi}\ d\theta\ \sin n\theta\
\tilde{\bm{r}}_{i}(\theta),$ (10) $\displaystyle\bm{b}_{i}^{(n)}$
$\displaystyle=$ $\displaystyle\frac{1}{\pi}\int_{0}^{2\pi}\ d\theta\ \cos
n\theta\ \tilde{\bm{r}}_{i}(\theta).$ (11)

If $\bm{a}_{i}^{(n)}=\bm{b}_{i}^{(n)}=\bm{0}$ for all $n$, the particle motion
is affine. When only $\bm{a}_{i}^{(1)}$ is non-zero, the non-affine trajectory
is a straight line, as shown in Fig. 4(a). In contrast, the trajectory
exhibits an ellipse when $\bm{b}_{i}^{(1)}$ is also non-zero, as shown in Fig.
4(b). A nontrivial trajectory, as shown in Fig. 1(a), contains modes with
$n\geq 2$. See Ref. Supple  for the relationship between the trajectories and
the Fourier coefficients.

Figure 4:  Schematics of the non-affine trajectory when only
$\bm{a}_{i}^{(1)}$ is non-zero (a) and only $\bm{a}_{i}^{(1)}$ and
$\bm{b}_{i}^{(1)}$ are non-zero (b).

In. Fig. 5(a), we plot the magnitudes of the Fourier components

$a^{(n)}=\sum_{i}\left\langle\left|\bm{a}_{i}^{(n)}\right|\right\rangle/N,\ \
b^{(n)}=\sum_{i}\left\langle\left|\bm{b}_{i}^{(n)}\right|\right\rangle/N$ (12)

obtained from our numerical data using Eqs. (10) and (11) against $n$ for
$\phi=0.87$ and $\gamma_{0}=0.01$ with $\omega\tau_{0}=10^{-4}$ and $10^{-5}$.
The Fourier components do not strongly depend on $\omega$, which indicates
that the nontrivial loops do not disappear in the limit $\omega\to 0$. For
different $\phi>\phi_{J}$ and $\gamma_{0}\geq 10^{-3}$, we have confirmed that
$a^{(1)}$ is always the largest Comment , the other modes are non-zero to make
loops with non-zero areas, and the Fourier components are independent of
$\omega$. In Fig. 5(b), we plot $a^{(n)}/\gamma_{0}$ and $b^{(n)}/\gamma_{0}$
against $\gamma_{0}$ for $\phi=0.87$ and $\omega\tau_{0}=10^{-4}$ with $n=1$,
where $a^{(n)}/\gamma_{0}$ and $b^{(n)}/\gamma_{0}$ are almost independent of
$\gamma_{0}$. This behavior is consistent with that for the number of contact
changes Supple .

---
Figure 5: (a) Magnitudes of Fourier coefficients $a^{(n)}$ and $b^{(n)}$
against $n$ for $\phi=0.87$ and $\gamma_{0}=0.02$ with
$\omega\tau_{0}=10^{-4}$ (filled symbols) and $10^{-5}$ (open symbols). (b)
Magnitudes of the Fourier coefficients $a^{(n)}$ and $b^{(n)}$ normalized by
$\gamma_{0}$ against $\gamma_{0}$ for $\phi=0.87$ and $\omega\tau_{0}=10^{-4}$
with $n=1$. $\phi=0.87$ corresponds to $\phi-\phi_{J}=0.029$.

Theoretical Analysis— Now, let us reproduce the numerical results by a simple
analytic calculation. Substituting Eq. (8) into Eq. (7),
${\bm{r}}_{ij}(\theta)$ is given by

$\displaystyle{\bm{r}}_{ij}(\theta)$ $\displaystyle=$
$\displaystyle\bm{R}_{ij}+\gamma_{0}Y_{ij}\sin\theta\bm{e}_{x}$ (13)
$\displaystyle+\sum_{n=1}^{\infty}\left(\bm{a}_{ij}^{(n)}\sin
n\theta+\bm{b}_{ij}^{(n)}\cos n\theta\right)$

Here, we define $\bm{a}^{(n)}_{ij}=\bm{a}^{(n)}_{i}-\bm{a}^{(n)}_{j}$,
$\bm{b}^{(n)}_{ij}=\bm{b}^{(n)}_{i}-\bm{b}^{(n)}_{j}$, and
$\bm{R}_{ij}=(X_{ij},Y_{ij})=\bm{R}_{i}-\bm{R}_{j}$. Substituting Eq. (13)
into Eq. (4) with Eq. (6) and neglecting the terms of $O(\gamma_{0})$, we
obtain the expression $G^{\prime}_{\rm T}$ of the storage modulus in SAS as
Supple

$\displaystyle G^{\prime}_{\rm T}$ $\displaystyle=$
$\displaystyle-\frac{1}{L^{2}}\sum_{i,j}\left\langle\frac{X_{ij}^{2}Y_{ij}^{2}}{R_{ij}}\Psi^{\prime}(R_{ij})\right\rangle$
(14) $\displaystyle-\frac{1}{L^{2}}\sum_{i,j}\left\langle
Y_{ij}^{2}\Psi(R_{ij})\right\rangle$
$\displaystyle-\frac{1}{L^{2}}\sum_{i,j}\left\langle\left(\frac{a_{ij,x}^{(1)}}{\gamma_{0}}Y_{ij}+X_{ij}\frac{a_{ij,y}^{(1)}}{\gamma_{0}}\right)\Psi(R_{ij})\right\rangle$
$\displaystyle-\frac{1}{L^{2}}\sum_{i,j}\left\langle
X_{ij}Y_{ij}\Psi^{\prime}(R_{ij})\frac{\bm{R}_{ij}\cdot\bm{a}_{ij}^{(1)}}{\gamma_{0}R_{ij}}\right\rangle,$

where $\Psi(r)=-U^{\prime}(r)/r$. Here, we have assumed
$|a_{i}^{(n)}|\sim|b_{i}^{(n)}|\sim\gamma_{0}$ and $\gamma_{0}\ll 1$. In the
expression of Eq. (14), only the first harmonic contribution from
$\bm{a}_{i}^{(1)}$ can survive because of Eq. (4). Note that $\bm{R}_{i}$ and
$\bm{a}^{(1)}_{i}$ cannot be determined within the theory but are determined
by our simulation data. In Fig. 2, we plot the theoretical prediction
$G^{\prime}_{\rm T}$ as open symbols. The theoretical prediction
quantitatively reproduces the numerical results except for large $\gamma_{0}$,
which is out of the scope of our theory. The first and second terms on the
right-hand side (RHS) of Eq. (14) represent the contributions from the affine
transformation depending only on $\bm{R}_{i}$, while the third and fourth
terms including $\bm{a}^{(1)}_{ij}/\gamma_{0}$ indicate the contributions from
the non-affine trajectories. As shown in Ref. Supple , the contributions from
the non-affine trajectories are almost independent of $\gamma_{0}$, which is
consistent with the behavior of $a^{(1)}/\gamma_{0}$ shown in Fig. 5(b).
Numerical evaluation in Ref. Supple  reveals that SAS is dominated by the
first term on RHS of Eq. (14) through the change of $\bm{R}_{i}$. The center
of the non-affine trajectories $\bm{R}_{i}$ is changed by the rearrangement of
the configuration during the transient to the absorbing state, which is
consistent with the memory formation of dense particles during oscillatory
shear Fiocco ; Paulsen ; Adhikari .

The theoretical expression $G^{\prime\prime}_{\rm T}$ of the loss modulus in
SAS is given by Supple

$\displaystyle G_{\rm T}^{\prime\prime}$ $\displaystyle=$
$\displaystyle-\frac{1}{L^{2}}\sum_{i,j}\left\langle\left(\frac{b_{ij,x}^{(1)}}{\gamma_{0}}Y_{ij}+X_{ij}\frac{b_{ij,y}^{(1)}}{\gamma_{0}}\right)\Psi(R_{ij})\right\rangle$
(15) $\displaystyle-\frac{1}{L^{2}}\sum_{i,j}\left\langle
X_{ij}Y_{ij}\Psi^{\prime}(R_{ij})R_{ij}\frac{\bm{R}_{ij}\cdot\bm{b}_{ij}^{(1)}}{\gamma_{0}R_{ij}^{2}}\right\rangle,$

where we have used the same assumption to obtain Eq. (14). Similar to the case
of $G_{T}^{\prime}$, only the contribution of the first harmonics
$\bm{b}_{i}^{(1)}$ in the expression of Eq. (8) can survive because of Eq.
(5). Note that $\bm{b}_{i}^{(1)}$ cannot be determined within the theory but
is evaluated by the simulation data. The loss modulus depends only on the non-
affine contribution including ${\bm{b}}_{i}^{(1)}$. The amplitude $b^{(1)}$
remains non-zero in the limit $\omega\to 0$, which leads to the residual loss
modulus as in Fig. 3(b). We plot the theoretical expression
$G^{\prime\prime}_{\rm T}$ using the open symbols in Fig. 3(a).
$G^{\prime\prime}_{\rm T}$ also reproduces the numerical results except for
large $\gamma_{0}$. Thus, our theory reveals the quantitative relationship
between the loss modulus and closed trajectories, which was suggested in Ref.
Keim14 .

Conclusion— We numerically studied the mechanical response of jammed materials
consisting of frictionless and overdamped particles under oscillatory shear.
The shear modulus exhibits SAS and the residual loss modulus exists in the
quasi-static limit in the absorbing state. Through Fourier analysis of the
closed trajectories, the theoretical expressions for the storage and loss
moduli quantitatively agree with the numerical results.

Reference Tighe11  reported that the loss modulus vanishes in the absorbing
jammed states in the limit $\omega\to 0$, which is inconsistent with our
result. It is noteworthy that Ref. Tighe11  did not consider any transient
state associated with contact changes before the system reaches the absorbing
state. Since the loss modulus is expected to be given by the generalized
Green-Kubo formula Chong ; Hayakawa , the origin of the residual loss modulus
might be plastic events in the transient dynamics.

Recent studies of large amplitude oscillatory shear (LAOS) reveal that there
are contributions from higher harmonics in the mechanical response of
nonlinear viscoelastic materials Wagner ; Hyun . We calculate nonlinear
viscoelastic moduli $G^{\prime}_{n}$ and $G^{\prime\prime}_{n}$ with $n\geq 2$
and confirm that such higher order moduli are negligible in our system as
shown in Ref. Supple .

In this Letter, we focus only on the nonlinear response of disordered
frictionless particles. However, even frictional grains and exhibit SAS
depending on the friction coefficient Otsuki21 . Therefore, an extension of
our theory to these systems will be our future work.

###### Acknowledgements.

The authors thank K. Saitoh, D. Ishima, T. Kawasaki, K. Miyazaki, and K.
Takeuchi for fruitful discussions. This work was supported by JSPS KAKENHI
Grants No. JP16H04025 and No. JP19K03670 and ISHIZUE 2020 of the Kyoto
University Research Development Program.

## References

  * (1) M. van Hecke, J. Phys. Condens. Matter 22, 033101 (2009)
  * (2) R. P. Behringer and B. Chakraborty, Rep. Prog. Phys. 82 012601 (2019)
  * (3) C. S. O’Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. 88, 075507 (2002).
  * (4) B. P. Tighe, Phys. Rev. Lett. 107, 158303 (2011).
  * (5) M. Otsuki and H. Hayakawa, Phys. Rev. E 95, 062902 (2017).
  * (6) C. Coulais, A. Seguin, and O. Dauchot, Phys. Rev. Lett. 113, 198001 (2014).
  * (7) M. Otsuki and H. Hayakawa, Phys. Rev. E 90, 042202 (2014).
  * (8) K. Hima Nagamanasa, S. Gokhale, A. K. Sood, and R. Ganapathy, Phys. Rev. E 89, 062308 (2014).
  * (9) E. D. Knowlton, D. J. Pine, and L. Cipelletti, Soft Matter 10, 6931 (2014).
  * (10) T. Kawasaki and L. Berthier, Phys. Rev. E 94, 022615 (2016).
  * (11) P. Leishangthem, A. D. S. Parmar, and S. Sastry, Nat. Commun. 8, 14653 (2017).
  * (12) A. H. Clark, J. D. Thompson, M. D. Shattuck, N. T. Ouellette, and C. S. O’Hern, Phys. Rev. E 97, 062901 (2018).
  * (13) J. Boschan, S. Luding, and B. P. Tighe, Granul. Matter 21, 58 (2019).
  * (14) H. Hinrichsen, Adv. Phys. 49, 815 (2000).
  * (15) M. Henkel, H. Hinrichsen, and S. Lubeck, Non-equilibrium Phase Transition I: Absorbing Phase Transitions (Springer, Heidelberg, 2008).
  * (16) D. J. Pine, J. P. Collub, J. F. Brady, and A. M. Leshansky, Nature (London) 438, 997 (2005).
  * (17) L. Corté, P. M. Chaikin, J. P. Gollub, and D. J. Pine, Nature Phys. 4, 420 (2008).
  * (18) J. Boschan, D. Vågberg, E. Somfai, and B. P. Tighe, Soft Matter 12, 5450 (2016).
  * (19) D. Nakayama, H. Yoshino, and F. Zamponi,J. Stat. Mech. 2016 104001 (2016).
  * (20) T. Kawasaki and K. Miyazaki, arXiv:2003.10716.
  * (21) S. Dagois-Bohy, E. Somfai, B. P. Tighe, and M. van Hecke, Soft Matter 13, 9036 (2017).
  * (22) D. Ishima and H. Hayakawa, Phys. Rev. E 101, 042902 (2020).
  * (23) M. Lundberg, K. Krishan, N. Xu, C. S. O’Hern, and M. Dennin, Phys. Rev. E 77, 041505 (2008).
  * (24) C. F. Schreck, R. S. Hoy, M. D. Shattuck, and C. S. O’Hern, Phys. Rev. E 88, 052205 (2013).
  * (25) N. C. Keim and P. E. Arratia, Soft Matter 9, 6222 (2013).
  * (26) N. C. Keim and P. E. Arratia, Phys. Rev. Lett. 112, 028302 (2014).
  * (27) I. Regev, T. Lookman, and C. Reichhardt Phys. Rev. E 88, 062401 (2013).
  * (28) I. Regev, J. Weber, C. Reichhardt, K. A. Dahmen, and T. Lookman, Nat. Commun. 6, 8805 (2015).
  * (29) N. V. Priezjev, Phys. Rev. E 93, 013001 (2016).
  * (30) M. O. Lavrentovich, A. J. Liu, and S. R. Nagel, Phys. Rev. E 96, 020101(R) (2017).
  * (31) K. Nagasawa, K. Miyazaki and T. Kawasaki, Soft Matter 15, 7557 (2019).
  * (32) P. Das, H. A. Vinutha, and S. Sastry, Proc. Natl. Acad. Sci. USA 117, 10203 (2020).
  * (33) K. Khirallah, B. Tyukodi, D. Vandembroucq, and C. E. Maloney, Phys. Rev. Lett. 126, 218005 (2021).
  * (34) M. S. van Deen, J. Simon, Z. Zeravcic, S. Dagois-Bohy, B. P. Tighe, and M. van Hecke, Phys. Rev. E 90, 020202(R) (2014).
  * (35) D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids 2nd ed. (Cambridge University Press, Cambridge, 2008).
  * (36) See Supplemental Material.
  * (37) M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Oxford University Press, Oxford, 1986).
  * (38) M. Meyers and K. Chawla, Mechanical Behavior of Materials (Cambridge University Press, Cambridge, 2008).
  * (39) We suppose $a^{(1)}$ is the largest because the mode proportional to $a^{(1)}$ is synchronized with the external oscillation $\sin\theta$.
  * (40) D. Fiocco, G. Foffi, and S. Sastry, Phys. Rev. Lett. 112, 025702 (2014).
  * (41) J. D. Paulsen, N. C. Keim, and S. R. Nagel, Phys. Rev. Lett. 113, 068301 (2014).
  * (42) M. Adhikari and S. Sastry, Eur. Phys. J. E 41, 105 (2018).
  * (43) S-H. Chong, M. Otsuki, and H. Hayakawa, Phys. Rev. E 81, 041130 (2010).
  * (44) H. Hayakawa and M. Otsuki, Phys. Rev. E 88, 032117 (2013).
  * (45) M. H. Wagner , V. H. Rolón-Garrido, K. Hyun, and M. Wilhelm, J. Rheol. 55, 495 (2011).
  * (46) K. Hyun, M. Wilhelm, C. O. Klein, K. S. Cho, J. G. Nam, K. H. Ahn, S. J. Lee, R. H. Ewoldt, and G. H. McKinley, Prog. Polym. Sci. 36, 1697 (2011).
  * (47) M. Otsuki and H. Hayakawa, Eur. Phys. J. E 44, 70 (2021).

Supplemental Material:

This Supplemental Material provides some details that are not written in the
main text. The results for underdamped frictionless granular particles without
background friction are presented in Sec. I. In Sec. II, we show the
dependence of $G^{\prime}$ and $G^{\prime\prime}$ on the number of particles
$N$ and the number of cycles $n_{c}$. In Sec. III, we present the time
evolution of the displacements of particles before reaching the absorbing
state and the evaluation of the yielding strain amplitude $\gamma_{c}$. In
Sec. IV, we illustrates the time evolutions of the stress-strain curves in the
absorbing and plastic states. In Sec. V, we show how particle trajectories
depend on $\gamma_{0}$ and $\omega$. In Sec. VI, we show that trajectories in
the absorbing state with longer periods do not affect our theoretical results
based on the absorbing trajectories whose periods are identical to the period
of the external oscillation. In Sec. VII, we present the loss modulus in the
absorbing and plastic states. In Sec. VIII, we demonstrate how the naive
result of the Kelvin–Voigt model can be recovered for sufficiently small
strain amplitude. In Sec. IX, we illustrate the relation between the Fourier
coefficients and the shape of particle trajectories. In Sec. X, we show the
number of contact changes during the last cycle in the absorbing state. In
Sec. XI, we derive Eqs. (14) and (15) in the main text. In Sec. XII, we
decompose the storage and loss moduli into several components, and clarify
what components are dominant contributions for the storage and loss moduli. In
Sec. XIII, we show the nonlinear viscoelastic moduli in our system to clarify
the roles of higher harmonics.

##  I Underdamped granular particles

In this section, we show that our results are qualitatively unchanged in
underdamped frictionless granular particles without background friction. Here,
we use the SLLOD equation given by [34]

$\displaystyle\frac{d}{dt}{\bm{r}}_{i}$ $\displaystyle=$
$\displaystyle\dot{\gamma}(t)y_{i}\bm{e}_{x}+\bm{p}_{i},$ (S1)
$\displaystyle\frac{d}{dt}{\bm{p}}_{i}$ $\displaystyle=$
$\displaystyle-\dot{\gamma}(t)p_{i,y}\bm{e}_{x}+\bm{F}_{i}$ (S2)

under the Lees–Edwards boundary condition, where
$\bm{p}_{i}=m(\dot{\bm{r}}_{i}-\dot{\gamma}(t)y_{i})\bm{e}_{x}$ and

$\bm{F}_{i}=-\sum_{j\neq
i}\frac{\partial}{\partial\bm{r}_{i}}U(r_{ij})-\sum_{j\neq i}\eta
v^{\rm(n)}_{ij}\Theta(d_{ij}-r_{ij})\frac{\bm{r}_{ij}}{r_{ij}}$ (S3)

with mass $m$, the interaction potential $U(r_{ij})$ given by Eq. (2), the
viscous constant $\eta$, and the normal velocity

$\displaystyle
v^{\rm(n)}_{ij}=\left\\{\frac{d}{dt}\bm{r}_{i}-\frac{d}{dt}\bm{r}_{j}\right\\}\cdot\frac{\bm{r}_{ij}}{r_{ij}}.$
(S4)

We adopt $\eta=\sqrt{mk}$. This model corresponds to frictionless granular
particles with the restitution coefficient $e=0.043$. We adopt the leapfrog
algorithm using the time step $\Delta t=0.05t_{0}$ with the characteristic
time with $t_{0}=\sqrt{m/k}$.

In Fig. S1, we plot non-affine trajectories in the last two cycles for a
particle with $\gamma_{0}=0.02$, $\phi=0.87$, and
$\omega=10^{-4}\tau_{0}^{-1}$ in the absorbing state. The trajectories are
closed, but there is a region where the position of the particle depends on
the number of cycles $n_{c}$. The result of Fig. S1 is a typical one from the
inertia effect in the underdamped system.

Figure S1:  Non-affine trajectories in the last two cycles for underdamped
particles with $\gamma_{0}=0.02$, $\phi=0.87$, and
$\omega=10^{-4}\tau_{0}^{-1}$ in the absorbing state. The circles represent
the trajectory in the last cycle. The line represents the trajectory in the
second to the last cycle.

Figure S2 plots the magnitude of the Fourier coefficients $a^{(n)}$ and
$b^{(n)}$ against $n$ for $\phi=0.87$ and $\gamma_{0}=0.02$ with
$\omega\tau_{0}=10^{-4}$. As in the overdamped system, $a^{(1)}$ takes the
largest value, and the other components are non-zero. In the inset of Fig. S2,
we show $a^{(1)}$ and $b^{(1)}$ against $\gamma_{0}$ for $\phi=0.87$ and
$\gamma_{0}=0.02$ with $\omega\tau_{0}=10^{-4}$, which are proportional to
$\gamma_{0}$.

Figure S2: Magnitudes of Fourier coefficients $a^{(n)}$ and $b^{(n)}$ of the
underdamped particles against $n$ for $\phi=0.87$ and $\gamma_{0}=0.02$ with
$\omega\tau_{0}=10^{-4}$. (Inset) Magnitudes of the Fourier coefficients
$a^{(n)}$ and $b^{(n)}$ against $\gamma_{0}$ for $\phi=0.87$ and
$\omega\tau_{0}=10^{-4}$ with $n=1$. The dashed line represents
$a^{(n)},b^{(n)}\propto\gamma_{0}$.

In Fig. S3, we plot the scaled storage modulus $G^{\prime}$ of the underdamped
particles in the absorbing state against the scaled amplitude $\gamma_{0}$ for
$\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870$ and $0.860$. The storage
modulus exhibits SAS. The corresponding theoretical expression
$G^{\prime}_{\rm T}$ in Eq. (14) as open symbols is also presented in Fig. S3,
which quantitatively reproduces the numerical results except for the region of
quite large $\gamma_{0}$

Figure S3: Scaled storage modulus $G^{\prime}/\sqrt{\phi-\phi_{J}}$ of
underdamped particles (filled symbols) and its theoretical expression using
$G^{\prime}_{\rm T}$ (open symbols) in Eq. (14) against scaled
$\gamma_{0}/(\phi-\phi_{J})$ in the absorbing state scaled by the distance
$\phi-\phi_{J}$ from the jamming point for $\omega=10^{-4}\tau_{0}^{-1}$ with
$\phi=0.870$ and $0.860$.

In Fig. S4(a), we present the loss modulus $G^{\prime\prime}$ obtained in our
simulation and its theoretical expression $G^{\prime\prime}_{\rm T}$ in Eq.
(15) against $\gamma_{0}$ for $\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870$
and $0.860$ in the underdamped system. The loss modulus $G^{\prime\prime}$
does not strongly depend on $\phi$ and $\gamma_{0}$, and the theoretical
expression $G^{\prime\prime}_{\rm T}$ reproduces the numerical results. In
Fig. S4(b), we plot the loss modulus $G^{\prime\prime}$ against $\omega$ for
$\phi=0.87$ with $\gamma=0.01$. The loss modulus seems to converge to a non-
zero value in the limit $\omega\to 0$.

Figure S4: (a) Loss modulus $G^{\prime\prime}$ of underdamped particles
obtained in our simulation (filled symbols) and its theoretical expression
$G^{\prime\prime}_{\rm T}$ (open symbols) in Eq. (14) against $\gamma_{0}$ for
$\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870$ and $0.860$. (b) Loss modulus
$G^{\prime\prime}$ against $\omega$ for $\phi=0.87$ with $\gamma_{0}=0.01$.

The results in this section are consistent with those in the main text for the
overdamped system. This indicates that the results presented in the main text
are universal for jammed disordered materials.

##  II Dependence of $G^{\prime}$ and $G^{\prime\prime}$ on $N$ and $n_{c}$

In this section, we show the dependence of $G^{\prime}$ and $G^{\prime\prime}$
on the numbers of particles $N$ and cycles $n_{c}$ for the overdamped dynamics
discussed in the main text. In Figs. S5(a) and (b), we plot $G^{\prime}$ and
$G^{\prime\prime}$ against $\gamma_{0}$ for $\omega=10^{-4}\tau_{0}^{-1}$,
$\phi=0.870$, and $n_{c}=20$ with $N=1000$ and $4000$, respectively. The shear
moduli $G^{\prime}$ and $G^{\prime\prime}$ for $N=1000$ and $4000$ are
consistent within error bars.

---
Figure S5: (a) Storage modulus $G^{\prime}$ against $\gamma_{0}$ for
$\omega=10^{-4}\tau_{0}^{-1}$, $\phi=0.870$, and $n_{c}=20$ with $N=1000$ and
$4000$. (b) Loss modulus $G^{\prime\prime}$ against $\gamma_{0}$ for
$\omega=10^{-4}\tau_{0}^{-1}$, $\phi=0.870$, and $n_{c}=20$ with $N=1000$ and
$4000$.

Figures S6(a) and (b) show $G^{\prime}$ and $G^{\prime\prime}$ against $n_{c}$
for $\omega=10^{-4}\tau_{0}^{-1}$, $\phi=0.870$, and $N=1000$ with
$\gamma_{0}=0.02,0.04$, and $0.08$, respectively. The shear moduli
$G^{\prime}$ and $G^{\prime\prime}$ reach statistical steady states for
$n_{c}\geq 20$ within error bars.

---
Figure S6: (a) Storage modulus $G^{\prime}$ against $n_{c}$ for
$\omega=10^{-4}\tau_{0}^{-1}$, $\phi=0.870$, and $N=1000$ with
$\gamma_{0}=0.02,0.04$, and $0.08$. (b) Loss modulus $G^{\prime\prime}$
against $n_{c}$ for $\omega=10^{-4}\tau_{0}^{-1}$, $\phi=0.870$, and $N=1000$
with $\gamma_{0}=0.02,0.04$, and $0.08$.

##  III Particle displacement and yielding strain amplitude

In this section, we show the time evolution of the displacements of particles
before reaching the absorbing state and the evaluation of the yielding strain
amplitude $\gamma_{c}$. Here, we introduce the particle displacement between
$n_{c}$-th and $(n_{c}-m)$-th cycles as

$dr_{m}(n_{c})=\sum_{i=1}^{N}\left|\bm{r}_{i}(n_{c}T)-\bm{r}_{i}((n_{c}-m)T)\right|/N$
(S5)

with the period $T=2\pi/\omega$. We define $dr(n_{c})$ as the minimum value of
$dr_{m}(n_{c})$ for $m$. In Fig. S7, we plot $dr(n_{c})$ against $n_{c}$ for
$\omega=10^{-4}\tau_{0}^{-1}$ and $\phi=0.860$ with $\gamma_{0}=0.08,0.04$,
and $0.02$. For $\gamma_{0}=0.08$, $dr(n_{c})$ remains non-zero, while it
approaches $0$ after a transient for $\gamma_{0}=0.04$ and $0.02$. It is
noteworthy that $dr(n_{c})$ reaches a steady state for $n_{c}>50$ in the case
of $\gamma_{0}=0.04$, which is much larger than the steady $n_{c}$ in Fig. S6.

Figure S7:  Displacements of particles $dr(n_{c})$ against $n_{c}$ for
$\omega=10^{-4}\tau_{0}^{-1}$ and $\phi=0.860$ with $\gamma_{0}=0.08,0.04$,
and $0.02$.

In Fig. S8, we plot $dr(n_{c})$ against $\gamma_{0}$ at $n_{c}=100$ for
$\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870,0.860,0.850$ and $0.845$. For
all $\phi$, $dr(n_{c})$ changes from $0$ to non-zero values as $\gamma_{0}$
increases. We call the absorbing state for $dr(n_{c})<d_{c}$ with smaller
$\gamma_{0}$ and the plastic state for $dr(n_{c})>d_{c}$ with larger
$\gamma_{0}$. The yielding strain amplitude $\gamma_{c}$ is defined as the
boundary between these states. From Fig. S8, we estimate
$0.04<\gamma_{c}<0.08$ for $\phi=0.870$, $0.04<\gamma_{c}<0.05$ for
$\phi=0.860$, $0.02<\gamma_{c}<0.03$ for $\phi=0.850$, and
$0.01<\gamma_{c}<0.02$ for $\phi=0.845$.

Figure S8:  Displacement of particles $dr(n_{c})$ against $\gamma_{0}$ at
$n_{c}=100$ for $\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870,0.860,0.850$
and $0.845$.

##  IV Stress-strain curve

In this section, we present typical stress strain curves in the absorbing and
plastic states including their time evolution. Figure S9 displays the shear
stress $\sigma$ against the strain $\gamma$ with $\gamma_{0}=0.02$ for
different $n_{c}$. For $n_{c}\leq 5$, the stress-strain curves are not
convergent, which indicate the system is in a transient state. For $n_{c}=6$
and $7$, the stress-strain curves become identical in the absorbing state. We
plot the shear stress $\sigma$ against the strain $\gamma$ for
$\gamma_{0}=0.1$ in Fig. S10. All the stress-strain curves are different for
all $n_{c}$ because the system is in the plastic state.

Figure S9:  Plots of shear stress $\sigma$ against $\gamma$ for
$\gamma_{0}=0.02$, $\omega=10^{-4}\tau_{0}^{-1}$, and $\phi=0.87$
corresponding to $\phi-\phi_{J}=0.029$ with various $n_{c}$.  Figure S10:
Plots of shear stress $\sigma$ against $\gamma$ for $\gamma_{0}=0.1$,
$\omega=10^{-4}\tau_{0}^{-1}$ and $\phi=0.87$ corresponding to
$\phi-\phi_{J}=0.029$ with various $n_{c}$.

##  V Dependence of trajectories on $\gamma_{0}$ and $\omega$

In this section, we present how particle trajectories depend on $\gamma_{0}$
and $\omega$. In Fig. S11, we plot the non-affine particle trajectories in the
last cycle for $\omega=10^{-3}\tau_{0}^{-1}$ and $10^{-5}\tau_{0}^{-1}$ with
$\phi=0.87$ and $\gamma_{0}=0.01$. Let us introduce

${\bm{r}}^{\prime}_{i}=(x^{\prime}_{i},y^{\prime}_{i})=\tilde{\bm{r}}_{i}-{\bm{R}}_{i}.$
(S6)

The trajectories form nontrivial loops, which remain for smaller $\omega$.

---
Figure S11:  Non-affine particle trajectories in the last cycle for
$\omega=10^{-3}\tau_{0}^{-1}$ (a) and $10^{-5}\tau_{0}^{-1}$ (b) with
$\phi=0.87$ and $\gamma_{0}=0.01$.

Figure S12 represents the non-affine particle trajectories in the last cycle
for $\omega=10^{-3}\tau_{0}^{-1}$ and $10^{-5}\tau_{0}^{-1}$ with $\phi=0.87$
and $\gamma_{0}=1.0\times 10^{-7}$. In Fig. S12 (a), the trajectory with
$\omega=10^{-3}\tau_{0}^{-1}$ form an ellipse, but the trajectory becomes a
straight line for $\omega=10^{-5}\tau_{0}^{-1}$ in Fig S12 (b).

---
Figure S12:  Non-affine particle trajectories in the last cycle for
$\omega=10^{-3}\tau_{0}^{-1}$ (a) and $10^{-5}\tau_{0}^{-1}$ (b) with
$\phi=0.87$ and $\gamma_{0}=1.0\times 10^{-7}$.

##  VI Effect of trajectories with longer periods

In this section, we discuss the effect of closed trajectories with periods
longer than $2\pi$. As indicated by Refs. [27, 28, 30, 31, 33], some samples
exhibit non-trivial absorbing trajectories where particles return to their
original positions after more than one cycle of oscillatory shear. In these
samples, the non-affine trajectories of a particle $\bm{r}_{i}(\theta)$
satisfy

$\bm{r}_{i}(\theta)=\bm{r}_{i}(\theta+2M\pi)$ (S7)

with $M=2,3,4,\cdots$. In this case, $\bm{r}_{i}(\theta)$ for $0<\theta<2M\pi$
is expressed in the Fourier series as

$\tilde{\bm{r}}_{i}(\theta)=\bm{R}^{\prime}_{i}+\sum_{m=1}^{\infty}\left(\bm{A}_{i}^{(m)}\sin\frac{m\theta}{M}+\bm{B}_{i}^{(m)}\cos\frac{m\theta}{M}\right)$
(S8)

with

$\displaystyle\bm{R}^{\prime}_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{2M\pi}\int_{0}^{2M\pi}\ d\theta\
\tilde{\bm{r}}_{i}(\theta),$ (S9)

and the Fourier coefficients

$\displaystyle\bm{A}_{i}^{(m)}$ $\displaystyle=$
$\displaystyle\frac{1}{M\pi}\int_{0}^{2M\pi}\ d\theta\ \sin\frac{m\theta}{M}\
\tilde{\bm{r}}_{i}(\theta),$ (S10) $\displaystyle\bm{B}_{i}^{(m)}$
$\displaystyle=$ $\displaystyle\frac{1}{M\pi}\int_{0}^{2M\pi}\ d\theta\
\cos\frac{m\theta}{M}\ \tilde{\bm{r}}_{i}(\theta).$ (S11)

However, in Eqs. (4) and (5), we need $\bm{r}_{i}(\theta)$ for
$0\leq\theta<2\pi$ to calculate $G^{\prime}$ and $G^{\prime\prime}$. When
$\bm{r}_{i}(\theta)$ is restricted to $0\leq\theta<2\pi$, we can use Eq. (8)
with the Fourier coefficient given by Eqs. (10) and (11) as an expression of
the trajectory, and we obtain the theoretical expressions Eqs. (14) and (15)
even in this case. It should be noted that samples in the absorbing state with
longer periods are rare, and the probability of emerging such a trajectory is
smaller than $0.01$ for sufficiently packed systems above the jamming point,
as shown in Ref. [30]. Therefore, we can ignore the effect of rare samples.

##  VII Loss modulus

In Fig. S13, we plot the loss modulus $G^{\prime\prime}$ against $\gamma_{0}$
for $\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870$ and $0.860$ including the
data in the absorbing and plastic states. This figure corresponds to Fig. 3(a)
in the main text, but Fig. S13 contains the data for a wide range of
$\gamma_{0}$. The previous studies [21, 22] reported that the loss modulus has
a peak around the yield strain for an underdamped system, but the peak of
$G^{\prime\prime}$ is not clearly visible in our overdamped system.

Figure S13:  Loss modulus $G^{\prime\prime}$ against $\gamma_{0}$ for
$\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870$ and $0.860$. The larger
(smaller) filled symbols represent the data in the absorbing (plastic) state.
The open pentagons represent the yield strain amplitude $\gamma_{c}$.

##  VIII Shear modulus for small $\gamma_{0}$

In this section, we demonstrate that $G^{\prime}$ and $G^{\prime\prime}$ obey
the Kelvin–Voigt model for a sufficiently small $\gamma_{0}$. Figure S14 is a
set of plots of $G^{\prime}$ and $G^{\prime\prime}$ against $\omega\tau_{0}$
for $\phi=0.870$ and $\gamma_{0}=1.0\times 10^{-7}$, where $G^{\prime}$ is
almost independent of $\omega$ and $G^{\prime\prime}$ is proportional to
$\omega$. This behavior is consistent with that of the Kelvin–Voigt model.

Figure S14: Plots of $G^{\prime}$ and $G^{\prime\prime}$ against $\omega$ for
$\phi=0.870$ and $\gamma_{0}=1.0\times 10^{-7}$. The dashed line represents
$G^{\prime\prime}\propto\omega$.

##  IX Relationship between closed trajectories and the Fourier coefficients

In this section, we present how the trajectory of a particle depends on the
Fourier coefficients. Figure S15 compares the trajectory of a particle
corresponding to Fig. 1(a) with its approximate trajectory using Eq. (8) with
some restricted modes, where we estimate the coefficients using the true
trajectory. In Fig. S15 (a), we plot the approximate trajectory (blue filled
circles) using only $\bm{a}_{i}^{\rm(1)}$, where we set the other coefficients
to $0$. The approximate trajectory (blue filled circles) is a straight line.
Figure S15 (b) shows the approximate trajectory using $\bm{a}_{i}^{\rm(1)}$
and $\bm{b}_{i}^{\rm(1)}$, where the trajectory becomes an ellipse. As we
increase the number of modes, the approximate trajectory approaches the true
trajectory, as shown in Figs. S15 (c) and (d).

Figure S15: Trajectory shown in Fig. 1(a) of the main text and its approximate
trajectories with some restricted modes. The red solid lines represent the
original data, and the blue filled circles represent the approximate
trajectory using (a) $\bm{a}_{i}^{\rm(1)}$, (b) $\bm{a}_{i}^{\rm(1)}$ and
$\bm{b}_{i}^{\rm(1)}$, (c) $\bm{a}_{i}^{\rm(n)}$ and $\bm{b}_{i}^{\rm(n)}$
with $n=1$ and $2$, and (d) $\bm{a}_{i}^{\rm(n)}$ and $\bm{b}_{i}^{\rm(n)}$
with $n=1,2,$ and $3$.

##  X Cyclic contact changes

In this section, we show the number of contact changes during the last cycle
in the absorbing state. References [21, 23, 25, 26] demonstrate that the
nontrivial loops originate from cyclic open and close contacts. Here, we
define $N_{cc}$ as the number of events where the same contact opens and
closes again during the last cycle. In Fig. S16, we present $N_{\rm cc}$
during the last cycle for $\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870$
against $\gamma_{0}$ in the absorbing state. The number of cyclic contact
changes $N_{\rm cc}$ is nearly proportional to $\gamma_{0}$. This dependence
is consistent with the behaviors of $a^{(n)}$ and $b^{(n)}$ of the Fourier
components, which are almost proportional to $\gamma_{0}$.

Figure S16:  The number of cyclic contact changes $N_{\rm cc}$ during the last
cycle for $\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870$ against
$\gamma_{0}$. The solid line represents $N_{\rm cc}\sim\gamma_{0}$.

##  XI Details of the theoretical analysis

In this section, we derive Eqs. (14) and (15) in the main text by assuming
$|a_{i}^{(n)}|\sim|b_{i}^{(n)}|\sim\gamma_{0}$ and $\gamma_{0}\ll 1$. From Eq.
(13) in the main text, $x_{ij}(\theta)$ and $y_{ij}(\theta)$ are given by

$\displaystyle x_{ij}=X_{ij}+\gamma_{0}\sin\theta
Y_{ij}+\sum_{n=1}^{\infty}\left(a_{ij,x}^{(n)}\sin n\theta+b_{ij,x}^{(n)}\cos
n\theta\right),$ (S12) $\displaystyle
y_{ij}=Y_{ij}+\sum_{n=1}^{\infty}\left(a_{ij,y}^{(n)}\sin
n\theta+b_{ij,y}^{(n)}\cos n\theta\right),$ (S13)

where
$\bm{a}^{(n)}_{ij}=(a_{ij,x}^{(n)},a_{ij,y}^{(n)})=\bm{a}^{(n)}_{i}-\bm{a}^{(n)}_{j}$,
$\bm{b}^{(n)}_{ij}=(b_{ij,x}^{(n)},b_{ij,y}^{(n)})=\bm{b}^{(n)}_{i}-\bm{b}^{(n)}_{j}$.
Using this equation and neglecting the terms of $O(\gamma_{0}^{2})$,
$|\bm{r}_{ij}(\theta)|^{2}=x_{ij}^{2}+y_{ij}^{2}$ is given by

$\displaystyle|\bm{r}_{ij}(\theta)|^{2}\simeq
R_{ij}^{2}\left\\{1+2E_{ij}(\theta)\right\\}$ (S14)

with

$\displaystyle E_{ij}(\theta)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{\bm{R}_{ij}\cdot\bm{a}_{ij}^{(n)}}{R_{ij}^{2}}\sin
n\theta+\sum_{n=1}^{\infty}\frac{\bm{R}_{ij}\cdot\bm{b}_{ij}^{(n)}}{R_{ij}^{2}}\cos
n\theta$ (S15)
$\displaystyle+\gamma_{0}\frac{X_{ij}Y_{ij}}{R_{ij}^{2}}\sin\theta.$

From Eq. (S14), $r_{ij}(\theta)$ is approximately obtained as

$\displaystyle r_{ij}(\theta)\simeq R_{ij}\left\\{1+E_{ij}(\theta)\right\\}$
(S16)

up to $O(\gamma_{0})$. Using this equation, we obtain
$\Psi(r)=-U^{\prime}(r)/r$ up to $O(\gamma_{0})$ as

$\displaystyle\Psi(r_{ij}(\theta))$ $\displaystyle\simeq$
$\displaystyle\Psi(R_{ij})+\Psi^{\prime}(R_{ij})R_{ij}E_{ij}(\theta).$ (S17)

Substituting Eqs. (S12)–(S17) into Eq. (6), we obtain

$\displaystyle\sigma(\theta)$ $\displaystyle=$
$\displaystyle-\frac{1}{L^{2}}\sum_{(i,j)}\left\\{\Psi(R_{ij})+\Psi^{\prime}(R_{ij})R_{ij}E_{ij}(\theta)\right\\}$
(S18) $\displaystyle\times\left\\{X_{ij}+\gamma_{0}\sin\theta
Y_{ij}+\sum_{n=1}^{\infty}\left(a_{ij,x}^{(n)}\sin n\theta+b_{ij,x}^{(n)}\cos
n\theta\right)\right\\}$
$\displaystyle\times\left\\{Y_{ij}+\sum_{n=1}^{\infty}\left(a_{ij,y}^{(n)}\sin
n\theta+b_{ij,y}^{(n)}\cos n\theta\right)\right\\}.$

Here, we abbreviate $\displaystyle\sum_{i}\sum_{j>i}$ as
$\displaystyle\sum_{(i,j)}$. Neglecting the terms of $O(\gamma_{0}^{2})$,
$\sigma(\theta)$ is approximated as

$\displaystyle\sigma(\theta)$ $\displaystyle\simeq$
$\displaystyle-\frac{1}{L^{2}}\sum_{(i,j)}X_{ij}Y_{ij}\Psi(R_{ij})-\frac{1}{L^{2}}\sum_{(i,j)}\gamma_{0}\sin\theta
Y_{ij}^{2}\Psi(R_{ij})$ (S19)
$\displaystyle-\frac{1}{L^{2}}\sum_{(i,j)}\sum_{n=1}^{\infty}\left(a_{ij,x}^{(n)}\sin
n\theta+b_{ij,x}^{(n)}\cos n\theta\right)Y_{ij}\Psi(R_{ij})$
$\displaystyle-\frac{1}{L^{2}}\sum_{(i,j)}\sum_{n=1}^{\infty}X_{ij}\left(a_{ij,y}^{(n)}\sin
n\theta+b_{ij,y}^{(n)}\cos n\theta\right)\Psi(R_{ij})$
$\displaystyle-\frac{1}{L^{2}}\sum_{(i,j)}X_{ij}Y_{ij}\Psi^{\prime}(R_{ij})R_{ij}E_{ij}(\theta).$

By substituting this equation into Eqs. (4) and (5) and using

$\displaystyle\frac{1}{\pi}\int_{0}^{2\pi}\ d\theta\ \sin m\theta\sin
n\theta=\delta_{mn},$ (S20) $\displaystyle\frac{1}{\pi}\int_{0}^{2\pi}\
d\theta\ \sin m\theta\cos n\theta=0,$ (S21)

we obtain Eqs. (14) and (15) in the main text.

##  XII Components of shear moduli

In this section, we clarify what terms of the theoretical expressions
$G^{\prime}_{\rm T}$ and $G^{\prime\prime}_{\rm T}$ in the absorbing state in
Eqs. (14) and (15) are dominant. Here, $G^{\prime}_{\rm T}$ consists of four
terms as

$\displaystyle G^{\prime}_{\rm T}$ $\displaystyle=$ $\displaystyle
G^{\prime}_{{\rm T},1}+G^{\prime}_{{\rm T},2}+G^{\prime}_{{\rm
T},3}+G^{\prime}_{{\rm T},4}$ (S22)

with

$\displaystyle G^{\prime}_{{\rm T},1}$
$\displaystyle=-\frac{1}{L^{2}}\sum_{i,j}\left\langle\frac{X_{ij}^{2}Y_{ij}^{2}}{R_{ij}}\Psi^{\prime}(R_{ij})\right\rangle,$
(S23) $\displaystyle G^{\prime}_{{\rm T},2}$
$\displaystyle=-\frac{1}{L^{2}}\sum_{i,j}\left\langle
Y_{ij}^{2}\Psi(R_{ij})\right\rangle,$ (S24) $\displaystyle G^{\prime}_{{\rm
T},3}$
$\displaystyle=-\frac{1}{L^{2}}\sum_{i,j}\left\langle\left(\frac{a_{ij,x}^{(1)}}{\gamma_{0}}Y_{ij}+X_{ij}\frac{a_{ij,y}^{(1)}}{\gamma_{0}}\right)\Psi(R_{ij})\right\rangle,$
(S25) $\displaystyle G^{\prime}_{{\rm T},4}$
$\displaystyle=-\left\langle\frac{1}{L^{2}}\sum_{i,j}X_{ij}Y_{ij}\Psi^{\prime}(R_{ij})\frac{\bm{R}_{ij}\cdot\bm{a}_{ij}^{(1)}}{\gamma_{0}R_{ij}}\right\rangle,$
(S26)

where $G^{\prime}_{{\rm T},1}$ and $G^{\prime}_{{\rm T},2}$ represent the
contributions from the affine motion, respectively, while $G^{\prime}_{{\rm
T},1}$ and $G^{\prime}_{{\rm T},2}$ are the contributions from the non-affine
motion, respectively. In Fig. S17, we show $G^{\prime}_{{\rm T},n}$ in the
absorbing state for $\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870$. We find
that $G^{\prime}_{{\rm T},1}$ and $G^{\prime}_{{\rm T},4}$ are dominant.
$G^{\prime}_{{\rm T},1}$ decreases with $\gamma_{0}$, while the other $G_{{\rm
T},n}^{\prime}$ with $n=2,3,4$ are almost independent of $\gamma_{0}$. This
indicates that SAS results from the behavior of $G^{\prime}_{{\rm T},1}$.

Figure S17:  $G^{\prime}_{{\rm T},n}$ with $n=1,2,3$ and $4$ against
$\gamma_{0}$ in the absorbing state for $\omega=10^{-4}\tau_{0}^{-1}$ with
$\phi=0.870$. The horizontal lines represent $G^{\prime}_{{\rm T},n}$ in the
limit $\gamma_{0}\to 0$, which is estimated at $\gamma_{0}=0.001$.

On the other hand, the loss modulus $G^{\prime\prime}_{\rm T}$ consists of two
terms as

$\displaystyle G_{\rm T}^{\prime\prime}$ $\displaystyle=$ $\displaystyle
G_{{\rm T},1}^{\prime\prime}+G_{{\rm T},2}^{\prime\prime}$ (S27)

with

$\displaystyle G_{{\rm T},1}^{\prime\prime}$
$\displaystyle=-\frac{1}{L^{2}}\sum_{i,j}\left\langle\left(\frac{b_{ij,x}^{(1)}}{\gamma_{0}}Y_{ij}+X_{ij}\frac{b_{ij,y}^{(1)}}{\gamma_{0}}\right)\Psi(R_{ij})\right\rangle,$
(S28) $\displaystyle G_{{\rm T},2}^{\prime\prime}$
$\displaystyle=-\frac{1}{L^{2}}\sum_{i,j}\left\langle
X_{ij}Y_{ij}\Psi^{\prime}(R_{ij})R_{ij}\frac{\bm{R}_{ij}\cdot\bm{b}_{ij}^{(1)}}{\gamma_{0}R_{ij}^{2}}\right\rangle.$
(S29)

In Fig. S17, we show $G^{\prime\prime}_{{\rm T},n}$ with $n=1$ and $2$ in the
absorbing state for $\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870$. The
result shows that $G^{\prime\prime}_{{\rm T},1}$ is dominant and almost
independent of $\gamma_{0}$. $G^{\prime\prime}_{{\rm T},2}$ depends on
$\gamma_{0}$, but it is much smaller than $G^{\prime\prime}_{{\rm T},1}$ for
$\gamma_{0}<0.1$.

Figure S18:  $G^{\prime}_{{\rm T},n}$ with $n=1$ and $2$ against $\gamma_{0}$
in the absorbing state for $\omega=10^{-4}\tau_{0}^{-1}$ with $\phi=0.870$.

##  XIII Non-linear viscoelastic moduli

In this section, we examine the non-liner viscoelastic moduli in our system.
The nonlinear elastic response is generally characterized by nonlinear
viscoelastic moduli $G^{\prime}_{n}$ and $G^{\prime\prime}_{n}$ satisfying
[41, 42]

$\sigma(t)=\gamma_{0}\sum_{n=1}\\{G^{\prime}_{n}\sin(n\omega
t)+G^{\prime\prime}_{n}\cos(n\omega t)\\}.$ (S30)

The storage and loss moduli are, respectively, given by
$G^{\prime}=G^{\prime}_{1}$ and $G^{\prime\prime}=G^{\prime\prime}_{1}$.
$G^{\prime}_{n}$ and $G^{\prime\prime}_{n}$ for $n\geq 2$ represent higher
harmonics. In Figs. S19 and S20, we plot $G^{\prime}_{n}$ and
$G^{\prime\prime}_{n}$ in the absorbing state for
$\omega=10^{-4}\tau_{0}^{-1}$ and $\phi=0.870$ with $n=1,2$, and $3$,
respectively. These figures indicate that the higher harmonics are negligible
in our system.

Figure S19:  $G_{n}^{\prime}$ against $\gamma_{0}$ in the absorbing state for
$\omega=10^{-4}\tau_{0}^{-1}$ and $\phi=0.870$ with $n=1,2$, and $3$.  Figure
S20:  $G_{n}^{\prime\prime}$ against $\gamma_{0}$ in the absorbing state for
$\omega=10^{-4}\tau_{0}^{-1}$ and $\phi=0.870$ with $n=1,2$, and $3$.