Brendan Walsh*1, Mark Hamilton*1,2, Greg Newby3, Xi Wang1, Serena Ruan1, Sheng
Zhao1,
Lei He1, Shaofei Zhang1, Eric Dettinger1, William T. Freeman2,4, Markus
Weimer1

# Large-Scale Automatic Audiobook Creation

###### Abstract

An audiobook can dramatically improve a work of literature's accessibility and
improve reader engagement. However, audiobooks can take hundreds of hours of
human effort to create, edit, and publish. In this work, we present a system
that can automatically generate high-quality audiobooks from online e-books.
In particular, we leverage recent advances in neural text-to-speech to create
and release thousands of human-quality, open-license audiobooks from the
Project Gutenberg e-book collection. Our method can identify the proper subset
of e-book content to read for a wide collection of diversely structured books
and can operate on hundreds of books in parallel. Our system allows users to
customize an audiobook's speaking speed and style, emotional intonation, and
can even match a desired voice using a small amount of sample audio. This work
contributed over five thousand open-license audiobooks and an interactive demo
that allows users to quickly create their own customized audiobooks. To listen
to the audiobook collection visit https://aka.ms/audiobook.

##  1 Introduction

Audiobooks have become a popular way to consume literature, news, and other
publications. Audiobooks not only allow existing readers to be able to enjoy
content on the go, but can help make content accessible to communities such as
children, the visually impaired, and new language learners. Traditional
methods of audiobook production, such as professional human narration or
volunteer-driven projects like LibriVox, are time-consuming, expensive, and
can vary in recording quality. These factors make it difficult to keep up with
an ever-increasing rate of book publication. In contrast, automatic audiobook
creation is orders of magnitude faster, cheaper, and more consistent but has
historically suffered from the robotic nature of text-to-speech systems and
the challenge of deciding what text should not be read aloud (e.g. tables of
contents, page numbers, figures, and footnotes).

We present a system that overcomes both of the aforementioned challenges by
generating high-quality audiobooks from heterogeneous collections of online
e-books. In particular, our system combines recent advances in neural text-to-
speech, emotive reading, scalable computing, and automatic detection of
relevant text to create thousands of reasonable-sounding audiobooks. We
contribute over five thousand audiobooks totaling approximately thirty-five
thousand hours of speech to the open source. We also contribute a
demonstration app that allows conference attendees to create a custom
audiobook, read aloud in their own voice, from any book from the collection
using only a few seconds of example sound.

##  2 Related Work

LibriVox is a well-known project that creates open-license audiobooks using
human volunteers. Although it has made significant contributions to the
accessibility of audiobooks, the quality of the produced audiobooks can be
inconsistent due to the varying skills and recording environments of the
volunteers. Furthermore, the scalability of the project is limited by the
availability of volunteers and the time it takes to record and edit a single
audiobook. Private platforms such as Audible create high-quality audiobooks
but do not release their works openly and charge users for their audiobooks.
Project Gutenberg hosts a broad collection of free e-books and a few
audiobooks. Their existing audiobooks feature a robotic text-to-speech voice
which limits listen-ability.

Text-to-speech is a well-studied problem and recent deep learning methods such
as WaveNet [1], Tacotron [2], and Fastspeech [3] have shown considerable
progress towards generating speech that rivals human quality and naturalness.
In contrast, the problem of selecting which text to read from an e-book has
received considerably less attention. Nevertheless, recent work by [4] has
explored whether it's possible to predict the ``start reading location'' using
LSTM-based models but does not tackle the cleaning of other irrelevant text
throughout the body of an e-book.

##  3 Methods

This work introduces a scalable system capable of converting HTML-based
e-books to high-quality audiobooks. Our pipeline is built using SynapseML[5],
a scalable machine learning framework that enables distributed orchestration
of the entire audiobook creation process.

###  3.1 Parsing e-Book HTML

Figure 1: t-SNE Representation of Clustered Ebooks. Colored areas represent
uniformly formatted clusters of books.

Our pipeline begins with thousands of free e-books provided by Project
Gutenberg. These e-books are provided in several different formats, and our
work focuses on their HTML format which is most amenable to automated parsing.
Parsing this extremely heterogeneous and diverse collection of e-books was the
most significant challenge we encountered. Project Gutenberg does not
standardize the contents of its HTML files and its e-books contain a
significant amount of text that would not be relevant for audio readers
including pre-ambles, tables of contents, tables, illustrations, in-text page
numbers, footnotes, transcriber notes, and other strange artifacts.

To create a high-quality subset of e-books we first featurize each e-book's
HTML Document Object Model (DOM) tree using a combination of automated (the
TF-IDF statistic on HTML Components) and hand-crafted HTML features. This
allowed us to cluster and visualize the entire collection of Project Gutenberg
HTML files and allowed us to find several large groups of commonly structured
files. We used these clusters of HTML files to build a rule-based HTML
normalizer that converted the largest classes of e-books into a standard
representation that could be automatically parsed. This analysis allowed us to
create a system that could quickly and deterministically parse a large number
of books. Most importantly it also allowed us to restrict attention to a
subset of files that would generate high-quality recordings when read. Figure
1 shows the results of this clustering procedure, demonstrating that several
clusters of similarly structured e-books naturally emerge in the Project
Gutenberg collection. Once parsed we can extract a stream of plain text to
feed to text-to-speech algorithms.

###  3.2 Generating High Quality Speech

Different audiobooks require different reading styles. Nonfiction works
benefit from a clear and neutral voice while fictional works with dialogue can
benefit from an emotive reading and some ``acting''. For the majority of the
books, we use a clear and neutral neural text-to-speech voice, However, in our
live demonstration we will present users with the ability to customize the
voice, speed, pitch, and intonation of the text.

To clone a user's voice we utilize zero-shot text-to-speech methods [6] to
efficiently transfer the voice characteristics from limited enrolled
recordings. This allows a user to quickly create an audiobook in their own
voice using a small amount of recorded audio.

To create an emotive reading of the text, we use an automatic speaker and
emotion inference system to dynamically change the reading voice and tone
based on context. This makes passages with multiple characters and emotional
dialogue more life-like and engaging. To this end, we first segment the text
into narration and dialogue and identify the speaker for each dialogue
section. We then predict the emotion of each dialogue using [7] in a self-
supervised manner. Finally, we assign separate voices and emotions to the
narrator and the character dialogues using the multi-style and contextual-
based neural text-to-speech model proposed in [8].

##  4 The Project Gutenberg Open Audiobook Collection

We introduce the Project Gutenberg Open Audiobook Collection: over five
thousand high-quality audiobooks generated from the Project Gutenberg
collection and available for free download and open use. We host these files
as a single zip file for the research community as well as on the major
podcast and audio file hosting platforms for use by the broader community.
This collection offers over thirty-five thousand hours of content including
classic literature, non-fiction, plays, and biographical works narrated in a
clear and consistent voice. We hope this contribution can provide value to
both the research community, and the broader community of audiobook listeners.

##  5 Demonstration

We plan to host a live demonstration application that allows conference
attendees to create their own custom audiobooks using our system. Users will
first start by selecting a book from the 5,000 titles in our collection using
a simple search interface. They can then select what voice they would like to
use for the recording from a large collection of existing neutral and emotion-
aware voices or even their own voice. If a user wants to create a custom
audiobook using their own voice, they will be asked to speak a few sentences
to quickly train a custom voice profile. Users will be able to listen to a
preview of their audiobook in real time and add an optional custom dedication
before submitting a larger job that reads the entire book. Once the pipeline
finishes we will email the user a link to download their custom-made
audiobook.

##  6 Conclusions

In this work, we present a novel pipeline to automate the creation of high-
quality audiobooks from heterogeneous e-books. Our system uses new advances in
neural text-to-speech, emotion recognition, custom voice cloning, and
distributed computing to create engaging and lifelike audiobooks. We apply
this system to donate over five thousand audiobooks to the open-source
community and aim to demonstrate this system by allowing conference attendees
to create custom audiobooks. We believe that this work has the potential to
greatly improve the accessibility and availability of audiobooks.

## References

  * [1] A. van den Oord _et al._ , ``Wavenet: A generative model for raw audio,'' 2016\.
  * [2] Y. Wang _et al._ , ``Tacotron: Towards end-to-end speech synthesis,'' in _Proc. INTERSPEECH 2017 Annual Conference of the International Speech Communication Association_ , 2017, pp. 4006––4010.
  * [3] Y. Ren _et al._ , ``Fastspeech: Fast, robust and controllable text to speech,'' in _Advances in Neural Information Processing Systems_ , 2019, pp. 3165––3174.
  * [4] S. B. Bodapati _et al._ , ``A machine learning approach to detecting start reading location of ebooks,'' in _2018 IEEE International Conference on Data Mining Workshops (ICDMW)_. IEEE, 2018, pp. 1522–1529.
  * [5] M. Hamilton _et al._ , ``Flexible and scalable deep learning with mmlspark,'' in _Proceedings of The 4th International Conference on Predictive Applications and APIs_ , ser. Proceedings of Machine Learning Research, vol. 82. PMLR, 24–25 Oct 2018, pp. 11–22.
  * [6] Y. Wu _et al._ , ``Adaspeech 4: Adaptive text to speech in zero-shot scenarios,'' in _Proc. INTERSPEECH 2022 Annual Conference of the International Speech Communication Association_ , 2022, pp. 2568––2572.
  * [7] ——, ``Self-supervised context-aware style representation for expressive speech synthesis,'' in _Proc. INTERSPEECH 2022 Annual Conference of the International Speech Communication Association_ , 2022, pp. 5503––5507.
  * [8] H. Guo _et al._ , ``Conversational end-to-end tts for voice agent,'' in _Proc. SLT 2021 IEEE Spoken Language Technology Workshop_ , 2021, pp. 403––409.


# Rheology of periodically sheared suspensions undergoing
reversible-irreversible transition

Zhouyang Ge1,2 Gwynn J. Elfring1 1Department of Mechanical Engineering and
Institute of Applied Mathematics, University of British Columbia, Vancouver
V6T 1Z4, BC, Canada  2Department of Engineering Mechanics, KTH Royal Institute
of Technology, 100 44 Stockholm, Sweden

(September 3, 2024)

###### Abstract

The rheology of non-colloidal suspensions under cyclic shear is studied
numerically. The main findings are a strain amplitude ($\gamma_{0}$) dependent
response in the shear stress and second normal stress difference ($N_{2}$).
Specifically, we find a reduced viscosity, an enhanced intracycle shear
thinning, the onset of a finite $N_{2}$ and its frequency doubling, all near a
critical strain amplitude $\gamma_{c}$ that scales with the volume fraction
$\phi$ as $\gamma_{c}\sim\phi^{-2}$. These rheological changes also signify a
reversible-irreversible transition (RIT), dividing stroboscopic particle
dynamics into a reversible absorbing phase (for $\gamma_{0}<\gamma_{c}$) and a
persistently diffusing phase (for $\gamma_{0}>\gamma_{c}$). We explain the
results based on two flow-induced mechanisms and elucidate their connection in
the context of RIT through the underlying microstructure, which tends towards
hyperuniformity near $\gamma_{0}=\gamma_{c}$. Overall, we expect this
correspondence between rheology and emergent dynamics to hold in a wide range
of settings where structural organizations are dominated by volume exclusions.

##  I Introduction

When a concentrated suspension is driven by shear, the suspending particles
can undergo diffusive motion even in the absence of inertial or Brownian
effects. Although this has long been observed Eckstein  _et al._ (1977) and
was subsequently understood as shear-induced diffusion Davis (1996), recent
studies on periodically sheared suspensions revealed a remarkable dynamical
phase transition that can occur without complex hydrodynamic interactions Pine
_et al._ (2005); Corté  _et al._ (2008). This phenomenon, now called
_reversible-irreversible transition_ (RIT) Menon and Ramaswamy (2009) or
absorbing-phase transition, refers to the change in particle dynamics viewed
stroboscopically: depending on the driving strain amplitude, the particle
motions are either reversible (absorbing) or irreversible (diffusive), with a
transition characterized by a critical strain amplitude. Because of the
ubiquity of suspensions in nature and their increasing relevance in
technology, new understanding of suspension dynamics under cyclic shear has
since led to broad developments in the studies of jamming and yielding Das
_et al._ (2020); Ness and Cates (2020); Wilken  _et al._ (2021), of self-
organization Corté  _et al._ (2009); Franceschini  _et al._ (2011); Royer and
Chaikin (2015) and memory formation Keim and Nagel (2011); Paulsen  _et al._
(2014); Keim  _et al._ (2019) in disordered systems, as well as in the
fabrication of hyperuniform materials Hexner and Levine (2015); Tjhung and
Berthier (2015); Wilken  _et al._ (2020) with potential applications in
photonics, among others Torquato (2018).

Unlike the clear connection between emergent dynamical and structural
properties of periodically sheared suspensions, less is known about the
relation between these dynamics and the suspension _rheology_. Conventionally,
rheological properties of complex fluids are characterized in either simple
shear flows (to obtain the nonlinear viscosity at different shear rates) or
oscillatory shear flows (to obtain the linear viscoelasticity at different
oscillation frequencies). For suspensions of non-colloidal spheres, in
particular, viscous effects usually dominate over elasticity Mewis and Wagner
(2012), and previous experiments on such systems have found that the
suspension viscosity in oscillatory shear can depend on either the strain
amplitude Breedveld  _et al._ (2001); Bricker and Butler (2006); Lin  _et al._
(2013); Ness  _et al._ (2017); Martone  _et al._ (2020), or both the strain
amplitude and the oscillation frequency Breedveld  _et al._ (2001); Narumi
_et al._ (2005); Martone  _et al._ (2020). Despite the complex dependence,
rheological properties of a suspension are theoretically fully determined from
its microstructure and the interactions among constituent particles. The
general understanding is that shearing induces anisotropies in the
microstructure Gadala-Maria and Acrivos (1980); Brady and Morris (1997), which
further leads to non-Newtonian behaviors (e.g. normal stress differences)
under simple shear Sierou and Brady (2002); Morris (2009), as well as strain-
dependent asymmetric stress responses under shear reversal Gadala-Maria and
Acrivos (1980). In addition, irreversible particle dynamics and rate-dependent
suspension rheology may arise from weak non-hydrodynamic interactions (e.g.
van der Waals forces), where the effect is typically stronger at smaller
strain amplitudes Ge  _et al._ (2021). Therefore, an indirect connection must
exist between rheology and emergent dynamics for suspensions under cyclic
shear, as both are related to the underlying microstructure. Although this
might have been realized to some extent in the past Breedveld  _et al._
(2001); Corté  _et al._ (2008); Lin  _et al._ (2013), a clear physical picture
of their relationship in the context of RIT has not been established so far.

In this paper, we identify generic rheological responses of periodically
driven suspensions that are also undergoing RIT. Specifically, we numerically
show that non-colloidal suspensions under cyclic shear have at least four
rheological signatures near their critical transition points: (I) reduced
relative viscosities, (II) enhanced intracycle shear thinning, (III) the onset
of a finite second normal stress difference and (IV) its peculiar frequency
doubling within a cycle. In addition, the critical strain amplitude
$\gamma_{c}$ over which the particle dynamics change most significantly can be
inferred from the above rheological measures and scales with the volume
fraction $\phi$ as $\gamma_{c}\sim\phi^{-2}$. The scaling suggests the
importance of many-body couplings which, counterintuitively, do not require
far-field hydrodynamic interactions. We explain these results on the basis of
the stress asymmetry caused by flow reversal Gadala-Maria and Acrivos (1980)
and the phenomenon of shear-induced diffusion Leighton and Acrivos (1987) that
were known experimentally but not fully recognized in the context of RIT.
Furthermore, microstructural analysis of the suspension near its critical
point links these observations to a common configuration characterized by a
local anisotropy and suppression of long-range density fluctuations; the
latter further implies hyperuniformity Torquato and Stillinger (2003);
Torquato (2018). Overall, our results highlight a nontrivial connection
between rheology and collective dynamics in periodically sheared suspensions.
More generally, we believe rheology may serve as an indicator of both the
underlying structure and the emergent dynamics.

##  II Model

We simulate particle suspensions in shear flow using a method called Fast
Stokesian Dynamics (FSD) Fiore and Swan (2019). FSD is a recent version of the
Stokesian Dynamics (SD) that is commonly used for modeling of (non-)colloidal
suspensions Brady and Bossis (1988). The core of SD is a multipole expansion
of the boundary integral equation for particles in Stokes flow, linearly
relating the hydrodynamic forces and torques ${\bm{F}}^{\text{H}}$ on the
particles to their translational and angular velocities ${\bm{U}}$ through a
resistance tensor ${\bm{R}}_{\text{FU}}$, viz.
${\bm{F}}^{\text{H}}=-{\bm{R}}_{\text{FU}}\cdot{\bm{U}}$. Assuming negligible
particle inertia, ${\bm{F}}^{\text{H}}$ is always in balance with the
prescribed external force ${\bm{F}}^{\text{Ext}}$ on each particle, viz.
${\bm{F}}^{\text{H}}+{\bm{F}}^{\text{Ext}}=0$; thus, the above equations are
closed provided ${\bm{R}}_{FU}$ and its inverse are known. All SD methods
decompose the hydrodynamic resistance into a far-field and a near-field part,
viz.
${\bm{R}}_{\text{FU}}={\bm{R}}_{\text{FU}}^{\text{ff}}+{\bm{R}}_{\text{FU}}^{\text{nf}}$,
allowing for independent control of the two contributions. Specifically, the
far-field interactions can be turned off by setting the off-diagonal terms of
${\bm{R}}_{\text{FU}}^{\text{ff}}$ to zero. The same decomposition also
applies to the stress calculation, which depends on both far-field and near-
field hydrodynamic interactions but does not affect the particle dynamics
Brady and Bossis (1988). FSD is a fast algorithm to numerically solve the
particle dynamics utilizing graphical processing units (GPU).

In our simulations we consider monodisperse suspensions made of 500 spheres
(except in one case where 4096 particles are considered) of radius $a$ at
three volume fractions $\phi=30\%$, 40%, and 50%; see Fig. 1. The
computational domain is a cubic box with Lees-Edwards boundary conditions,
where the particles are driven by oscillatory strain,
$\gamma(t)=\gamma_{0}\sin(\omega t)$, under various driving amplitudes
$\gamma_{0}$ and frequencies $\omega$. This naturally defines a force scale,
$F_{\text{stk}}\equiv 6\pi\eta_{0}\gamma_{0}\omega a^{2}$, being the Stokes
drag of an isolated particle in a fluid with dynamic viscosity $\eta_{0}$.
Since particles in liquids are usually charged (due to surface acids or bases,
adsorption of free ions, etc.), we impose a pairwise repulsive force between
nearby particles, viz. ${\bm{F}}^{\text{El}}={\bm{F}}_{0}\exp(-\kappa h)$,
where ${\bm{F}}_{0}$ is the maximal repulsion, $\kappa$ the inverse Debye
length, and $h$ the surface gap between two particles Mewis and Wagner (2012);
Israelachvili (2011). $F_{0}$ can be equated with $F_{\text{stk}}$ to obtain a
reference time scale to compare with that of the flow; cf. Ref. Ge  _et al._
(2021). Throughout this work we fix $F_{0}=10F_{\text{stk}}$ and
$\kappa^{-1}=0.01a$ to model a strong, short-range repulsion. Effectively,
this also keeps the non-dimensional shear rate (proportional to
$F_{\text{stk}}/F_{0}$) constant under different $\gamma_{0}$ and $\omega$.
The model produces quasi-Newtonian rheology with weak shear thickening under
steady shear (the suspension viscosity increases by less than a factor of 1.2
when the shear rate increases by 100 times). The source code of the FSD
algorithm was originally published in the Supplementary Materials of Ref.
Fiore and Swan (2019). The updated version we use is also available at
https://github.com/GeZhouyang/FSD.

Figure 1: Periodically sheared suspensions. The schematic shows 500 particles
in a cubic box of length $L$ at volume fraction $\phi=50\%$. The particles are
projected onto the shear plane, where the color indicates the instataneous
shear stress (reddish colors for higher stresses). Here, the stresses are
close to zero at maximum strains because the shear rate therein is zero.

##  III Results

In what follows, we discuss the rheology, collective dynamics, and underlying
microstructural evolution of periodically sheared suspensions. All simulations
at the same volume fraction are started from the same random initial
condition, and the results are averaged over all particles (since they are
identical) and the last few cycles (if the steady or quasi-steady state
statistics are reported). We have verified that the steady state results do
not depend on the initial condition appreciably, except when the suspension is
driven by very small strain amplitudes in which case the microstructure does
not evolve. In general, we expect the insights gained from our analysis to be
valid regardless of the model details.

###  III.1 Rheology

(a)(b)(c)

Figure 2: Viscosity of periodically sheared suspensions. (a) Evolution of the
relative viscosity $\eta_{r}^{\prime}$ and relative elastic response
$\eta_{r}^{\prime\prime}$ at $\phi=40\%$. Inset shows the corresponding total
viscous responses during the last half cycle; see Eq. (2). (b) Steady-state
$\eta_{r}^{\prime}$ under various driving amplitudes and volume fractions.
Gray circles are the experimental results at $\phi=40\%$ from Ref. Bricker and
Butler (2006). (c) Magnitude of the intracycle shear thinning quantified by
$-\eta^{\prime}_{3}/\eta^{\prime}_{1}$ (results at $\phi=40\%$ and 50% are
shifted up by 0.05 and 0.1, respectively, for clarity). In (b,c), filled
symbols are with far-field hydrodynamic interactions, empty symbols without.

####  III.1.1 Viscosity

One may decompose the shear stress $\sigma_{xy}(t)$ within an oscillation
cycle in Fourier series,

$\displaystyle\sigma_{xy}(t)=\gamma_{0}\sum_{n=1}^{\infty}\big{(}G_{n}^{\prime}\sin(n\omega
t)+\omega\eta_{n}^{\prime}\cos(n\omega t)\big{)},$ (1)

where $G_{n}^{\prime}$ and $\eta_{n}^{\prime}$ denote the $n$th elastic and
viscous coefficients, respectively. These coefficients vary from cycle to
cycle; but at steady state (when $\sigma_{xy}(t)$ becomes fully periodic), the
leading-order even modes are always much smaller than their preceding odd
modes, $G_{n-1}^{\prime}$ and $\eta_{n-1}^{\prime}$, reflecting the fore-aft
symmetry of the flow-induced microstructure Gadala-Maria and Acrivos (1980).
Further inspection of the power spectrum shows that $\eta_{n}^{\prime}\gg
G_{n}^{\prime}/\omega$ for all $n$, as well as rapid decays of both
$G_{n}^{\prime}$ and $\eta_{n}^{\prime}$ as $n$ increases (e.g.
$\eta_{5}^{\prime}/\eta_{1}^{\prime}\sim\mathcal{O}(10^{-2})$ typically).
These observations allow us to dramatically simplify Eq. (1) and approximate
the instantaneous shear viscosity as

$\displaystyle\eta(t)\approx\eta_{1}^{\prime}\cos(\omega
t)+\eta_{3}^{\prime}\cos(3\omega t),$ (2)

where $\eta\equiv\sigma_{xy}/\gamma_{0}\omega$. In Eq. (2),
$\eta_{1}^{\prime}$ and $\eta_{3}^{\prime}$ correspond to the _average_
viscosity and the leading order _intracycle_ viscosity variation,
respectively. While $\eta_{1}^{\prime}$ is strictly positive for suspensions
of passive particles, $\eta_{3}^{\prime}$ can be either positive or negative.
Under the oscillatory strain $\gamma(t)=\gamma_{0}\sin(\omega t)$, thus shear
rate $\dot{\gamma}(t)=\gamma_{0}\omega\cos(\omega t)$, it can be shown that a
negative (or positive) $\eta_{3}^{\prime}$ indicates shear thinning (or
thickening) within a cycle Ewoldt  _et al._ (2008). Note that this is not
because $\eta_{3}^{\prime}$ is the next Fourier coefficient in Eq. (2), but
rather it is deduced from an orthogonal stress decomposition Cho  _et al._
(2005); Ewoldt  _et al._ (2008).

In the following, we rescale $\eta_{1}^{\prime}$ by the underlying fluid
viscosity $\eta_{0}$ to obtain the conventional relative viscosity
$\eta_{r}^{\prime}\equiv\eta_{1}^{\prime}/\eta_{0}$; $\eta_{3}^{\prime}$ is
rescaled by $-\eta_{1}^{\prime}$ to highlight the effect of shear thinning;
lastly, we define the relative elastic response (same unit as
$\eta_{r}^{\prime}$) as $\eta_{r}^{\prime\prime}\equiv
G_{1}^{\prime}/\eta_{0}\omega$.

Fig. 2(a) shows the evolution of $\eta_{r}^{\prime}$ and
$\eta_{r}^{\prime\prime}$, as well as the instantaneous viscous response to
the imposed shear rate (obtaining the so-called Lissajous curves), under three
representative strain amplitudes $\gamma_{0}=0.2$, 0.8, and 3.2 at
$\phi=40\%$. Here, the duration is measured in total strain units
$\gamma_{\text{tot}}$ at the increment of full cycles; i.e.
$\gamma_{\text{tot}}=4\gamma_{0}n_{\text{cyc}}$, where $4\gamma_{0}$ is the
accumulated strain over one cycle, and $n_{\text{cyc}}$ is the number of
cycles. The Lissajous curves are plotted only during a half cycle, since both
$\eta(t)$ and $\dot{\gamma}(t)$ are anti-symmetric about $t|_{\dot{\gamma}=0}$
at steady state. These plots reveal a few interesting features of the rheology
under cyclic shear. First, it confirms that the stress response is
predominantly viscous, as $\eta_{r}^{\prime}\gg\eta_{r}^{\prime\prime}$ for
all $\gamma_{0}$. Second, it shows that the times to reach steady state are
long; there are continuous, albeit small, changes in both $\eta_{r}^{\prime}$
and $\eta_{r}^{\prime\prime}$ even after $\gamma_{\text{tot}}>100$ for
$\gamma_{0}<3.2$. Third, the instantaneous viscous response within the last
cycle shows qualitative differences under the three strain amplitudes: the
response is linear at $\gamma_{0}=0.2$, exhibits shear thinning at
$\gamma_{0}=0.8$, and becomes quasi-linear and hysteretic at $\gamma_{0}=3.2$.
The slow convergence to steady state and the three stress response patterns
are consistent with previous experiments and simulations Gadala-Maria and
Acrivos (1980); Bricker and Butler (2006); Martone  _et al._ (2020); Ness  _et
al._ (2017), suggesting that there might be a common rheological signature of
periodically sheared suspensions. Furthermore, because stress is uniquely
determined from microstructure in non-colloidal suspensions, this implies
glassy dynamics of suspensions under periodic strain, which are indeed
generally seen in disordered soft materials Sollich  _et al._ (1997).

To better examine the rheology, we show the steady-state $\eta_{r}^{\prime}$
and $-\eta_{3}^{\prime}/\eta_{1}^{\prime}$ under various strain amplitudes and
volume fractions in Fig. 2(b,c). Here, apart from the different behavior of
$\eta_{r}^{\prime}$ in the densest packing where the suspension tends to
crystallize under shear 111Shear-induced crystallization has been observed in
monodisperse systems at high packing fractions, cf. Refs. Sierou and Brady
(2002); Duff and Lacks (2007); Ghosh  _et al._ (2022). Numerically, this makes
the steady state difficult to obtain as the structural changes are slow. The
results for $\gamma_{0}>0.6$ at $\phi=50\%$ are thus taken at
$\gamma_{\text{tot}}=256$ rather than the true steady states., both viscous
coefficients are non-monotonic in $\gamma_{0}$: at each $\phi$, the relative
viscosity tends to a minimal value near an intermediate amplitude, where the
effect of intracycle shear thinning is roughly at a maximum (cf. Fig. 2(a),
inset). $\gamma_{0}$-dependent viscosity has been previously observed in
oscillatory shear experiments Breedveld  _et al._ (2001); Bricker and Butler
(2006); Lin  _et al._ (2013); Ness  _et al._ (2017); Martone  _et al._ (2020).
For example, Bricker and Butler performed extensive experiments of particle
suspensions at $\phi=40\%$, using four different types of materials in two
rheometer geometries, and found that the complex viscosity is smallest at
$\gamma_{0}\approx 1$ Bricker and Butler (2006, 2007). Our results are in
qualitative agreement with theirs (the precise value of $\eta_{r}^{\prime}$
depends not only on $\gamma_{0}$, but also on the physicochmeical interaction
between the particles); see Fig. 2(b).

The physical origin of the observed shear rheology lies in the structural
changes caused by _repeated_ flow reversals. Microscopically, shearing a
suspension creates a non-uniform structure characterized by anisotropic pair
distributions, which are correlated with higher particle stresses Brady and
Morris (1997); Sierou and Brady (2002); Morris (2009). The deviation of this
shear-induced microstructure from a given initial configuration, as well as
the extent to which it can be reversed within a cycle, depends on the applied
strain and the straining _history_. Crucially, history matters only
transiently and when $\gamma_{0}$ is “small”. The microstructure formed under
large $\gamma_{0}$ is statistically invariant within each half cycle; that is,
it always reaches the same steady state twice every period. The threshold
$\gamma_{0}$ above which history plays diminishing role in the rheology is
related to the RIT and will be quantified in Sec. III.2. As we shall see, it
is the combined effects of shearing with a finite amplitude and the resulting
microstructural evolution that give rise to the non-monotonic viscosity
variations on average and within a cycle.

####  III.1.2 Normal stress differences

Apart from a varying viscosity, suspensions can also develop anisotropic
normal stresses under shear, leading to finite normal stress differences
(NSDs). Rheologically, NSDs account for a number of non-Newtonian phenomena
(most notably the Weissenberg effect) and are a consequence of the flow-
induced microstructure Morris (2009). Below, we examine NSD responses of
periodically sheared suspensions.

Similar to the shear stress, any component of the normal stresses,
$\sigma_{\alpha\alpha}(t)$, may be decomposed within a cycle in Fourier
series,

$\displaystyle\sigma_{\alpha\alpha}(t)=\sum_{n=0}^{\infty}\hat{\sigma}_{n}^{(\alpha)}\exp(in\omega
t),\quad\alpha=x\textrm{ or }y\textrm{ or }z,$ (3)

where $\hat{\sigma}_{n}^{(\alpha)}$ denotes the $n$th complex Fourier
coefficient of $\sigma_{\alpha\alpha}$. The first and second NSD, defined
respectively as $N_{1}\equiv\sigma_{xx}-\sigma_{yy}$ and
$N_{2}\equiv\sigma_{yy}-\sigma_{zz}$, are thus also decomposed

$\displaystyle
N_{1}(t)=\sum_{n=0}^{\infty}\big{(}\hat{\sigma}_{n}^{(x)}-\hat{\sigma}_{n}^{(y)}\big{)}\exp(in\omega
t),$ (4a) $\displaystyle
N_{2}(t)=\sum_{n=0}^{\infty}\big{(}\hat{\sigma}_{n}^{(y)}-\hat{\sigma}_{n}^{(z)}\big{)}\exp(in\omega
t).$ (4b)

Note that, NSDs can have nonzero averages over a cycle and are even functions
of the applied strain at steady state. This is in contrast to the shear
stress, though the different behaviors are all ultimately due to the
underlying time-reversal symmetry.

Figure 3: Average responses of the first normal stress difference of
periodically sheared suspensions. Inset shows the pure hydrodynamic
contributions.

Fig. 3 shows the average $N_{1}$ normalized by $\eta_{0}\gamma_{0}\omega$
under various $\gamma_{0}$ and $\phi$ at steady state. The most salient
feature is that $N_{1}$ remains negligible for small $\gamma_{0}$ until it
becomes finite as $\gamma_{0}$ further increases. A negligible $N_{1}$
suggests that the pair distribution function characterizing suspension
microstructure is nearly fore-aft symmetric, as would be the case in the
purely hydrodynamic limit Morris (2009). This ceases to be the case when
$\gamma_{0}$ is beyond a $\phi$-dependent threshold, where $N_{1}$ first
becomes positive (ever so slightly), then negative. A positive $N_{1}$ can be
due to either a higher tensile normal stress in the flow direction, or more
collisions (which contribute to a negative normal stress) in the velocity
gradient direction. Our simulations suggest the latter, as the hydrodynamic
contributions to $N_{1}$ are always negative (see Fig. 3, inset). Lastly, we
note that a transition from positive to negative $N_{1}$ has been observed
before and may have certain rheological implications Bricker and Butler
(2007). Whether this is a generic behavior or can be neglected remains to be
explored.

(a)(b)(c)

Figure 4: The second normal stress difference of periodically sheared
suspensions. (a) Average $N_{2}$ per cycle, where the dashed lines are linear
fits of the last few data points (filled symbols are the entire stress, empty
symbols are contributions from collisions and electrostatic repulsions). (b)
Ratio of the magnitudes of the fourth ($|\nu_{4}|$) and second ($|\nu_{2}|$)
harmonic responses (inset shows the ratio due to hydrodynamic interactions).
(c) Lissajous curves at $\phi=40\%$ for three representative $\gamma_{0}$
(arrows show the direction of the curves).

We now turn to $N_{2}$, which is consistently reported negative and dominates
over $N_{1}$ (i.e. $|N_{2}|>|N_{1}|$) for suspensions Guazzelli and Pouliquen
(2018), but has received somewhat less attention Maklad and Poole (2021).
Dividing $N_{2}$ by $\gamma_{0}\omega$ and keeping three leading-order even
modes, its instantaneous response is approximated as

$\displaystyle\nu(t)\approx\nu_{0}+\nu_{2}\exp(2i\omega
t)+\nu_{4}\exp(4i\omega t),$ (5)

where $\nu\equiv N_{2}/\gamma_{0}\omega$. Note that, the coefficients
$\nu_{0}$, $\nu_{2}$, and $\nu_{4}$ are not necessarily in decaying magnitude,
but the higher order terms are always decaying (i.e.
$|\nu|_{n+2}\ll|\nu|_{n}$, $n=4,6,8,\dots$) upon inspection of the power
spectrum. In the following, we rescale $\nu_{0}$ by $-\eta_{0}$ to examine the
average $N_{2}$ per cycle, and compare the magnitudes of the second
($|\nu_{2}|$) and fourth ($|\nu_{4}|$) harmonic responses to analyze the
nonlinear response.

Fig. 4(a) shows the average $N_{2}$ for all $\phi$ and $\gamma_{0}$ at steady
state. There are clearly two regimes: under small $\gamma_{0}$, $N_{2}\approx
0$ as in the case of $N_{1}$, indicating an overall isotropic microstructure
in all principal directions; under large $\gamma_{0}$, $N_{2}<0$ as particles
collide more in the shear plane than in the flow-vorticity plane, consistent
with typical observations in simple shear flow Guazzelli and Pouliquen (2018)
(notice the hollow markers in Fig. 4a). The onset of a finite $N_{2}$ occurs
at a well-marked $\gamma_{0}$, which reduces with $\phi$ and, importantly, is
within the range of the corresponding $\gamma_{0}$ where the minimum
$\eta_{r}^{\prime}$ (prior to crystallization, if any) and the maximum
$-\eta_{3}^{\prime}/\eta_{1}^{\prime}$ are also obtained, cf. Fig. 2(b,c). As
will be shown in Sec. III.2, this is not a coincidence but rather an
indication of RIT.

The most unexpected behavior of $N_{2}$ is its nonlinear response near the
threshold strain amplitude. As shown in Fig. 4(b), this nonlinear effect is so
strong that $|\nu_{4}|/|\nu_{2}|$ can increase by three orders-of-magnitude.
Therefore, the responding frequency of $N_{2}$ under an imposed frequency
$\omega$ will double for a range of $\gamma_{0}$ (i.e. $2\omega\to 4\omega$).

To understand this behavior, we examine the Lissajous curves of the normalized
instantaneous second normal stress difference ($\nu(t)/\eta_{0}$) vs. strain
($\gamma(t)/\gamma_{0}$) for three representative $\gamma_{0}$ at $\phi=40\%$
in Fig. 4(c). At $\gamma_{0}=0.2$, $\nu(t)/\eta_{0}$ exhibits the expected
periodic response with a $2\omega$ frequency: it is anti-symmetric about both
the zero strain (due to reversibility) and the zero shear rate (due to even
symmetry). At $\gamma_{0}=3.2$, $\nu(t)/\eta_{0}$ is no longer zero at the
zero strain (due to persistent particle collisions, consistent with its non-
zero average), but still responds to the driving strain at $2\omega$.
Interestingly, at $\gamma_{0}=0.8$ where the response frequency doubles,
$\nu(t)/\eta_{0}$ is initially positive as $\gamma(t)/\gamma_{0}$ increases
from 0, but only up to a fraction of the strain amplitude, then it becomes
negative as $\gamma(t)/\gamma_{0}$ further increases; see arrows in Fig. 4(c).
A positive $N_{2}$ can only be due to higher tension in the direction aligning
with the velocity gradient than the vorticity, since collisions occur more
frequently in the shear plane than the others, as mentioned earlier. Indeed,
the stress contributions from hydrodynamic interactions alone exhibit the same
nonlinear response (see Fig. 4(b) inset). This is in stark contrast to the
average response, which is dominated by the interparticle repulsion and
collisions.

Finally, we note that all aforementioned rheological responses, namely, the
reduced relative viscosity, the enhanced intracycle shear thinning, the onset
of a finite second normal stress difference and its frequency doubling within
a cycle, _all_ occur within a narrow range of $\gamma_{0}$ that depends on
$\phi$ only. We will show that the $\phi$-dependent $\gamma_{0}$ coincides
with the critical strain amplitude $\gamma_{c}$ of RIT next.

###  III.2 Dynamics

(a)(b)(c)

Figure 5: Dynamics of periodically sheared suspensions. (a) Evolution of the
per-cycle diffusivity $D_{\text{cyc}}$ in total strain $\gamma_{\text{tot}}$
under various strain amplitudes $\gamma_{0}$ (labelled aside) at $\phi=40\%$.
Inset shows the particle dispersion under $\gamma_{0}=0.8$ after half a cycle.
(b) Final $D_{\text{cyc}}$ evaluated at $\gamma_{\text{tot}}=256$ under
various $\gamma_{0}$ and $\phi$ (filled symbols are with far-field
hydrodynamic interactions, empty symbols without). The shaded regions
correspond to the range of $\gamma_{0}$ estimated from the rheological
signatures according to Figs. 2(b,c) and 4(a,b) (see text). (c) Phase diagram
of the reversible-irreversible transition. The empty circles are estimated
bounds (shades in b) of the critical amplitude $\gamma_{c}$, whereas the
dashed line corresponds to $\gamma_{c}=0.14\phi^{-1.93}$ Pine  _et al._
(2005).

The $\gamma_{0}$-dependent transient and steady-state rheology implies
irreversible particle dynamics. If the particle trajectories are completely
reversible from the beginning, we should not expect to observe any stress
evolution from cycle to cycle, and both $\eta_{r}^{\prime}$ and
$-\eta_{3}^{\prime}/\eta_{1}^{\prime}$ are likely to change monotonically with
$\gamma_{0}$ (assuming the starting configuration is random). Irreversibility,
however, may last only _transiently_ as particles can self-organize into
reversible absorbing states according to a so-called random organization (RO)
model Corté  _et al._ (2008); Wilken  _et al._ (2021). In that model, the only
physical assumptions are (I) particles do not overlap and (II) there is some
source of microscopic irreversibility (or noise). Despite the simplicity, the
model applies to a wide range of systems subject to cyclic shear, and is
believed to be the underlying reason for the observed reversible-irreversible
transition.

To check irreversibility, we sample the particle statistics directly from
their trajectories in the simulations. Fig. 5(a) shows the transient self-
diffusivity under various strain amplitudes $\gamma_{0}$ at volume fraction
$\phi=40\%$, where the diffusivity $D_{\text{cyc}}$ is computed _per cycle_
according to

$\langle(\Delta_{i}/a)^{2}\rangle=2D_{\text{cyc}}\gamma_{\text{cyc}},$ (6)

where $\Delta_{i}$ is the displacement of particle $i$ over a cycle,
$\langle\cdot\rangle$ averages over all particles, and
$\gamma_{\text{cyc}}=4\gamma_{0}$. Clearly, the initial dynamics are
irreversible under all amplitudes. As a further illustration, the inset in
Fig. 5(a) shows the distribution of particle dispersions after the first
_half_ cycle under $\gamma_{0}=0.8$. Here, the particles not only fail to
return to their original positions, but are also displaced anisotropically due
to shear. After another half cycle, the preferential direction of the
dispersion will reverse (i.e. $\Delta\to-\Delta$). Repeated oscillations
eventually leads to either decaying or persistent diffusion, depending on
$\gamma_{0}$, as predicted by the RO model.

The final per-cycle diffusivity is plotted as a function of strain amplitudes
for all volume fractions in Fig. 5(b), where $D_{\text{cyc}}$ is evaluated at
$\gamma_{\text{tot}}=256$. The common pattern of these curves is a transition
from reversible to irreversible dynamics as $\gamma_{0}$ increases. For each
$\phi$, we highlight the range of $\gamma_{0}$ within which the suspension
exhibits the aforementioned rheological signatures illustrated in Figs. 2 and
4. This is done by intersecting the range of $\gamma_{0}$ estimated
individually. For example, at $\phi=40\%$, we estimate that the minimum
$\eta_{r}^{\prime}$ is obtained under $\gamma_{0}\in(0.8,1.6)$, whereas the
maximum $-\eta_{3}^{\prime}/\eta_{1}^{\prime}$ is obtained under
$\gamma_{0}\in(0.4,1.0)$; thus, the overlap $\gamma_{0}\in(0.8,1.0)$ provides
a bound for the critical amplitude $\gamma_{c}$. The onset of
$-\nu_{0}/\eta_{0}$ and the peak of $|\nu_{4}|/|\nu_{2}|$ also fall in the
same range. Despite the sparsity of data, these bounds roughly correspond to
the region where the dynamics are the most susceptible to changes in
$\gamma_{0}$, i.e. where $\partial D_{\text{cyc}}/\partial\gamma_{0}$ is the
steepest. Therefore, our results suggest that rheological measures can also be
used to infer the dynamical phase transition of suspensions under cyclic
shear, as plotted in Fig. 5(c).

The critical amplitude $\gamma_{c}$ for periodically sheared suspensions
considered herein depends only on the volume fraction. Previous experiments
report $\gamma_{c}=0.14\phi^{-1.93}$ by extrapolating the effective particle
diffusivities (based on particle mean square displacements over all cycles,
rather than per cycle) from $\gamma_{0}>\gamma_{c}$ Pine  _et al._ (2005). In
our simulations it is difficult to track the particle displacements under
large $\gamma_{0}$ over long time (e.g. $\gamma_{\text{tot}}>256$) due to the
finite size of the computational domain. However, our rheologically estimated
$\gamma_{c}$ agrees rather well with the experimental fit, as seen in Fig.
5(c); hence, it validates the alternative approach.

The scaling of $\gamma_{c}$ with respect to $\phi$ can be explained
qualitatively by the phenomenon of shear-induced diffusion. Earlier
experiments on non-colloidal suspensions in simple shear flow found that the
self-diffusivity $D_{s}$ scales with the particle size $a$, the imposed shear
rate $\dot{\gamma}$, and the volume fraction $\phi$ as
$D_{s}\sim\phi^{2}\dot{\gamma}a^{2}$ Leighton and Acrivos (1987). Here, the
$\phi^{2}$ dependence is due to many-body interactions, since the probability
of one particle interacting with two others is proportional to $\phi^{2}$ in a
statistically uniform suspension (pair interactions would lead to a scaling of
$D_{s}\sim\phi\dot{\gamma}a^{2}$, but they are insufficient to cause
persistent diffusion when the non-hydrodynamic force is weak Drazer  _et al._
(2002)). Replacing $\dot{\gamma}$ with $\gamma_{0}\omega$ and rescaling
$D_{s}$ by $a^{2}\omega$, we have $\tilde{D}_{s}\equiv
D_{s}/(a^{2}\omega)\sim\phi^{2}\gamma_{0}$; thus the scaling of
$\gamma_{c}\sim\phi^{-2}$ at a measurable onset diffusivity.
Counterintuitively, our simulations suggest that far-field hydrodynamic
interactions are _not_ responsible for the many-body effect, as the results
are qualitatively the same with and without such interactions; see Figs.
2(b,c) and 5(b). Effectively, local interactions such as the lubrication and
electrostatic repulsion are sufficient to couple the particle motions to yield
the observed rheology and collective dynamics in concentrated suspensions. By
contrast, far-field hydrodynamic interactions promote reversibility (evidently
in the extreme case of two particles in shear flow Batchelor and Green (1972);
da Cunha and Hinch (1996)), thus acting to smear the dynamical phase
transition Metzger and Butler (2010).

Finally, we note that the consequence of RIT on the rheology is non-trivial.
When $\gamma_{0}<\gamma_{c}$, the per-cycle dynamics are not immediately
quiescent, but rather decay in time, thus changing the reference condition at
the beginning of each cycle; this is manifested in the intracycle rheological
responses of both the viscosity and second normal stress difference. When
$\gamma_{0}>\gamma_{c}$, the particle dynamics depend less and less on the
strain history as collisions become persistent, leading to anisotropic
microstructures that are evident in the increased viscosity and finite second
normal stress difference. These two regimes have distinct rheological
behaviors which, at the onset of transition, is highly nonlinear. Below, we
examine the suspension microstructure at $\gamma_{0}\approx\gamma_{c}$ to
elucidate the connection between the rheology and emergent dynamics.

###  III.3 Microstructure

The correspondence between rheology and dynamics can only be due to the
underlying microstructure, which may be characterized by the pair correlation
function. As an example, we show in Fig. 6 the particle angular distributions
under $\gamma_{0}=0.8$ and $\phi=40$%, at three representative distances
corresponding to the first three peaks of the radial distribution function at
steady state. The initial random configuration shows the expected isotropic
distribution at all distances. However, the final distribution within a
boundary layer of thickness $0.1a$ is strongly anisotropic, with virtually no
particles at latitudes above $15\degree$ or below $-15\degree$, except for
those near longitude $\pm 90\degree$. Here, the equatorial ring corresponds to
regions of the smallest velocity difference, whereas the longitudinal ring is
perpendicular to the flow; particles in both of these regions experience less
relative motion under shear than elsewhere da Cunha and Hinch (1996). Since
particles tend to be displaced in the velocity gradient direction under shear
(see the inset in Fig. 5a), repeated oscillations lead to a local depletion of
particles from regions of greater encounter probability. The final effect is
the formation of a highly correlated short-range structure, with seemingly
little change in the bulk, as depicted in Fig. 6. We observe this phenomenon
in general for all cases; however, it is always the most pronounced at
$\gamma_{0}\approx\gamma_{c}$.

Figure 6: Particle angular distributions (Aitoff projection) under $\phi=40$%
and $\gamma_{0}=0.8$. (left) Initial distributions. (right) Final
distributions. In all cases, the angular distribution is averaged at distance
$r\pm 0.1a$, where darker colour indicates higher probability. The maps are
oriented such that the shear flow has the same velocity along latitudes.

To gain further insight into the bulk structure near the critical transition
point, we analyse the nearest-neighbour graph (NNG) obtained from snapshots of
the particle configuration. By definition, a NNG is formed by connecting each
particle to its nearest neighbour, resulting in a partition of the particles
into clusters (or minimum spanning trees as in graph theory) of various size
and shape. Specifically, we define the size as the number of particles in the
cluster, and characterize the shape from the moment-of-inertia tensor
${\bm{I}}_{c}$ of the (as if rigid) cluster; that is, we define the shape
factor as $\sqrt{I_{max}/I_{min}}$, where $I_{max}$ and $I_{min}$ are the
maximum and minimum eigenvalues of ${\bm{I}}_{c}$, respectively. Under these
definitions, the minimum size is 2 (i.e. a doublet), while the minimum shape
factor is 1 (i.e. an isotropic cluster). NNG is suitable for clustering
analysis due to its simplicity and analytical properties Eppstein  _et al._
(1997); here, we use it also because it is a unique and parameter-free
characterization of the suspension network.

(a)(b)

Figure 7: Size and shape distributions of particle clusters in its nearest-
neighbour graph under $\gamma_{0}=0.8$ at $\phi=40$%. (Top) Initial
statistics. (Bottom) Final statistics. (Inset, a) Average distance to the
$k$th nearest neighbour; blue squares for the initial distances, orange
circles final. (Inset, b) Evolution of the number variance scaling exponent
$\beta$ obtained under the same $\phi$ and $\gamma_{0}$, but in a much larger
computational box ($N=4096$). The dashed line is a fit of the data according
to Eq. (7).

Fig. 7 shows the size and shape distributions of particle clusters for a
suspension under $\gamma_{0}=0.8$ at $\phi=40\%$. The initial graph is
composed of clusters from a wide range of sizes, with the biggest cluster
approaching 10% of the sample size ($N=500$ in this case). The NNG at steady
state, however, is dominated by small clusters of less than 4 particles. The
dissolution of large clusters is a sign of collective organization, implying
that the distance between adjacent particles are enlarging throughout the
suspension. This can be seen directly in the inset (a) in Fig. 7, where we
plot the average distance to the $k$th nearest neighbour in the initial and
final configurations. The distance increases up to $k=7$, with the final
$r_{1}/a=2.1$. Correspondingly, a doublet made of two nearest spheres at
$r/a=2.1$ has a shape factor of 1.93, which explains the peak in the final
distribution of the shape factor.

Lastly, the uniform increase of local separations between particles at a fixed
volume implies suppression of density fluctuations at large length scales. The
latter is a hallmark of _hyperuniformity_ , which has been recently proposed
as a distinguishable state of matter Torquato and Stillinger (2003); Torquato
(2018). Here, density fluctuations can be quantified by the variance of the
number of particles within a randomly located sphere of radius $r$, viz.
$\sigma_{n}^{2}(r)\equiv\langle n(r)^{2}\rangle-\langle n(r)\rangle^{2}$, for
$r\in[2a,L/4]$, where $L$ is the length of the periodic box. Then, the
exponent of the scaling law $\sigma_{n}^{2}(r)\sim r^{\beta}$ measures the
degree of hyperuniformity: in $d$-dimensional space, $\beta=d$ corresponds to
non-hyperuniform structures (e.g. ordinary fluids), whereas materials with
$\beta\in[d-1,d)$ are considered hyperuniform (e.g. crystals have $\beta=d-1$;
other forms of scaling are also possible Torquato (2018)). We compute $\beta$
for the same parameters $\gamma_{0}=0.8$ and $\phi=40\%$, but without far-
field hydrodynamics in a much larger computational box ($N=4096$), as shown in
the inset (b) in Fig. 7. The data plotted in total strain
$\gamma_{\text{tot}}$ is well-fitted by

$\beta=(\beta_{0}-\beta_{\infty})(1+\gamma_{\text{tot}})^{\nu}+\beta_{\infty},$
(7)

where $\beta_{0}$ and $\beta_{\infty}$ are the initial and steady-state values
of $\beta$, respectively. The fitting yields $\nu=-0.29$ and
$\beta_{\infty}=2.0$, strongly suggesting that system will eventually become
hyperuniform. We note that the size of our largest simulation ($N=4096$,
$L/4=8.8a$) may still be too small to conclude on the large-scale structural
correlations. Nevertheless, the different characteristics of the suspension
microstructure are consistent and in agreement with previous predictions and
experiments Torquato (2018).

##  IV Concluding remarks

In summary, we have presented a numerical study of periodically sheared non-
colloidal suspensions under various oscillation amplitudes and at different
volume fractions. The main finding is a generic rheology characterized by a
reduced viscosity, an enhanced intracycle shear thinning, the onset of a
finite second normal stress difference as well as its frequency doubling
within a cycle, all near a volume fraction dependent oscillation amplitude.
This rheological behavior is closely related to the emergent dynamics of
suspensions under cyclic shear. In fact, we show that the amplitude inferred
from rheological measures also predicts the reversible-irreversible
transition. These findings are explained based on, firstly, the stress
asymmetry caused by flow reversal and, secondly, the shear-induced diffusion.
Although the basic phenomena have been studied separately in the past, the
connection between them in the context of RIT was unclear. This paper was
motivated to bridge the gap.

Finally, despite our focus on non-colloidal suspensions of frictionless hard
spheres of equal size, the general physical picture depicted here should hold
in a broader context. As long as the structural organization is dominated by
volume exclusions, there could always be a configuration of least density
fluctuation, thus the specific correspondence between rheology and emergent
dynamics, regardless of the particle shape, size distribution, surface
roughness, and so on. Furthermore, evolution towards a state of least density
fluctuation (or crystallization at high packing fractions Ghosh  _et al._
(2022)) seems favorable, since the microstructure associated with it reduces
the energy dissipation under shear. At present, there is a growing interest in
understanding the collective dynamics of _active suspensions_ that are made of
self-propelled particles such as bacteria or algae, where long-range
hydrodynamic interactions are of greater importance Elfring and Brady (2021).
Extending the current study in those settings may lead to fruitful future
research.

###### Acknowledgements.

The research is sponsored by Canada’s Department of National Defence through
Contract No. CFPMN1-026. Z.G. acknowledges support from the Swedish Research
Council (grant No. 2021-06669VR) and thanks S. Bagheri for hospitality.
Finally, the authors wish to thank S. Torquato, P.M. Chaikin, V. Shaik, M.
Boudina, and the anonymous referees for useful comments.

## References

  * Eckstein  _et al._ (1977) E. C. Eckstein, D. G. Bailey, and A. H. Shapiro, J. Fluid Mech. 79, 191 (1977).
  * Davis (1996) R. H. Davis, J. Fluid Mech. 310, 325 (1996).
  * Pine  _et al._ (2005) D. J. Pine, J. P. Gollub, J. F. Brady, and A. M. Leshansky, Nature 438, 1476 (2005).
  * Corté  _et al._ (2008) L. Corté, P. M. Chaikin, J. P. Gollub, and D. J. Pine, Nat. Phys. 4, 420 (2008).
  * Menon and Ramaswamy (2009) G. I. Menon and S. Ramaswamy, Phys. Rev. E 79, 061108 (2009).
  * Das  _et al._ (2020) P. Das, H. Vinutha, and S. Sastry, Proc. Natl. Acad. Sci. U. S. A. 117, 10203 (2020).
  * Ness and Cates (2020) C. Ness and M. E. Cates, Phys. Rev. Lett. 124, 088004 (2020).
  * Wilken  _et al._ (2021) S. Wilken, R. E. Guerra, D. Levine, and P. M. Chaikin, Phys. Rev. Lett. 127, 038002 (2021).
  * Corté  _et al._ (2009) L. Corté, S. J. Gerbode, W. Man, and D. J. Pine, Phys. Rev. Lett. 103, 248301 (2009).
  * Franceschini  _et al._ (2011) A. Franceschini, E. Filippidi, E. Guazzelli, and D. J. Pine, Phys. Rev. Lett. 107, 250603 (2011).
  * Royer and Chaikin (2015) J. R. Royer and P. M. Chaikin, Proc. Natl. Acad. Sci. U. S. A. 112, 49 (2015).
  * Keim and Nagel (2011) N. C. Keim and S. R. Nagel, Phys. Rev. Lett. 107, 010603 (2011).
  * Paulsen  _et al._ (2014) J. D. Paulsen, N. C. Keim, and S. R. Nagel, Phys. Rev. Lett. 113, 068301 (2014).
  * Keim  _et al._ (2019) N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry, and S. R. Nagel, Rev. Mod. Phys. 91, 035002 (2019).
  * Hexner and Levine (2015) D. Hexner and D. Levine, Phys. Rev. Lett. 114, 110602 (2015).
  * Tjhung and Berthier (2015) E. Tjhung and L. Berthier, Phys. Rev. Lett. 114, 148301 (2015).
  * Wilken  _et al._ (2020) S. Wilken, R. E. Guerra, D. J. Pine, and P. M. Chaikin, Phys. Rev. Lett. 125, 148001 (2020).
  * Torquato (2018) S. Torquato, Physics Reports 745, 1 (2018).
  * Mewis and Wagner (2012) J. Mewis and N. J. Wagner,  _Colloidal suspension rheology_ (Cambridge University Press, 2012).
  * Breedveld  _et al._ (2001) V. Breedveld, D. van den Ende, R. Jongschaap, and J. Mellema, J. Chem. Phys. 114, 5923 (2001).
  * Bricker and Butler (2006) J. M. Bricker and J. E. Butler, J. Rheol. 50, 711 (2006).
  * Lin  _et al._ (2013) Y. Lin, N. Phan-Thien, and B. C. Khoo, J. Rheol. 57, 1325 (2013).
  * Ness  _et al._ (2017) C. Ness, Z. Xing, and E. Eiser, Soft Matter 13, 3664 (2017).
  * Martone  _et al._ (2020) R. Martone, C. Carotenuto, and M. Minale, J. Rheol. 64, 1075 (2020).
  * Narumi  _et al._ (2005) T. Narumi, H. See, A. Suzuki, and T. Hasegawa, J. Rheol. 49, 71 (2005).
  * Gadala-Maria and Acrivos (1980) F. Gadala-Maria and A. Acrivos, J. Rheol. 24, 799 (1980).
  * Brady and Morris (1997) J. F. Brady and J. F. Morris, J. Fluid Mech. 348, 103 (1997).
  * Sierou and Brady (2002) A. Sierou and J. F. Brady, J. Rheol. 46, 1031 (2002).
  * Morris (2009) J. F. Morris, Rheol. acta 48, 909 (2009).
  * Ge  _et al._ (2021) Z. Ge, R. Martone, L. Brandt, and M. Minale, Phys. Rev. Fluids 6, L101301 (2021).
  * Leighton and Acrivos (1987) D. Leighton and A. Acrivos, J. Fluid Mech. 177, 109–131 (1987).
  * Torquato and Stillinger (2003) S. Torquato and F. H. Stillinger, Phys. Rev. E 68, 041113 (2003).
  * Fiore and Swan (2019) A. M. Fiore and J. W. Swan, J. Fluid Mech. 878, 544 (2019).
  * Brady and Bossis (1988) J. F. Brady and G. Bossis, Annu. Rev. Fluid Mech. 20, 111 (1988).
  * Israelachvili (2011) J. Israelachvili,  _Intermolecular and Surface Forces_ (Academic Press, 2011).
  * Ewoldt  _et al._ (2008) R. H. Ewoldt, A. Hosoi, and G. H. McKinley, J. Rheol. 52, 1427 (2008).
  * Cho  _et al._ (2005) K. S. Cho, K. Hyun, K. H. Ahn, and S. J. Lee, J. Rheol. 49, 747 (2005).
  * Sollich  _et al._ (1997) P. Sollich, F. Lequeux, P. Hébraud, and M. E. Cates, Phys. Rev. Lett. 78, 2020 (1997).
  * Note (1) Shear-induced crystallization has been observed in monodisperse systems at high packing fractions, cf. Refs. Sierou and Brady (2002); Duff and Lacks (2007); Ghosh  _et al._ (2022). Numerically, this makes the steady state difficult to obtain as the structural changes are slow. The results for $\gamma_{0}>0.6$ at $\phi=50\%$ are thus taken at $\gamma_{\text{tot}}=256$ rather than the true steady states.
  * Bricker and Butler (2007) J. M. Bricker and J. E. Butler, J. Rheol. 51, 735 (2007).
  * Guazzelli and Pouliquen (2018) E. Guazzelli and O. Pouliquen, J. Fluid Mech. 852, P1 (2018).
  * Maklad and Poole (2021) O. Maklad and R. Poole, J. Non-Newtonian Fluid Mech. 292, 104522 (2021).
  * Drazer  _et al._ (2002) G. Drazer, J. Koplik, B. Khusid, and A. Acrivos, J. Fluid Mech. 460, 307–335 (2002).
  * Batchelor and Green (1972) G. Batchelor and J.-T. Green, J. Fluid Mech. 56, 375 (1972).
  * da Cunha and Hinch (1996) F. da Cunha and E. Hinch, J. Fluid Mech. 309, 211 (1996).
  * Metzger and Butler (2010) B. Metzger and J. E. Butler, Phys. Rev. E 82, 051406 (2010).
  * Eppstein  _et al._ (1997) D. Eppstein, M. S. Paterson, and F. F. Yao, Discrete Comput. Geom. 17, 263 (1997).
  * Ghosh  _et al._ (2022) A. Ghosh, J. Radhakrishnan, P. M. Chaikin, D. Levine, and S. Ghosh, Phys. Rev. Lett. 129, 188002 (2022).
  * Elfring and Brady (2021) G. J. Elfring and J. F. Brady, arXiv:2107.03531 (2021).
  * Duff and Lacks (2007) N. Duff and D. J. Lacks, Phys. Rev. E 75, 031501 (2007).