# Relationships among structure, memory, and flow in sheared disordered
materials

K.L. Galloway1, E.G. Teich2, X-g Ma 3, Ch. Kammer1, I.R. Graham3, N.C. Keim4,
C. Reina1, D.J. Jerolmack5,1, A. G. Yodh3, P.E. Arratia∗,1
1Department of Mechanical Engineering and Applied Mechanics, University of
Pennsylvania
2Department of Bioengineering, University of Pennsylvania
3Department of Physics and Astronomy, University of Pennsylvania
4Department of Physics, Pennsylvania State University
5Department of Earth and Environmental Science, University of Pennsylvania

###### Abstract

A fundamental challenge for disordered solids is predicting macroscopic yield
from the microscopic arrangements of constituent particles. Yield is
accompanied by a sudden and large increase in energy dissipation due to the
onset of plastic rearrangements. This suggests that one path to understanding
bulk rheology is to map particle configurations to their mode of deformation.
Here, we perform laboratory experiments and numerical simulations that are
designed to do just that: 2D dense colloidal systems are subjected to
oscillatory shear, and particle trajectories and bulk rheology are measured.
We quantify particle microstructure using excess entropy. Results reveal a
direct relation between excess entropy and energy dissipation, that is
insensitive to the nature of interactions among particles. We use this
relation to build a physically-informed model that connects rheology to
microstructure. Our findings suggest a framework for tailoring the rheological
response of disordered materials by tuning microstructural properties.

##  I Introduction

Disordered solids are ubiquitous. They are found, for example, in our foods as
pastes and gels Nagel (2017), and amidst our homes in the form of concrete
Ioannidou _et al._ (2016) and mud Jerolmack and Daniels (2019); Nie _et al._
(2020). Frustratingly, these materials can experience sudden mechanical
failure, such as the collapse of soil during rapid mudslides. Indeed, when
sufficiently stressed, all disordered materials exhibit a swift decrease in
ability to support load. In the vicinity of this “yield” transition, the solid
material shifts from a state wherein energy is stored via internal elastic
forces, to a state in which energy is dissipated via irreversible plastic
rearrangements Larson (2010); Guazzelli and Pouliquen (2018); Chen _et al._
(2010); Buttinoni _et al._ (2017); Pham _et al._ (2008). Microscopic
spatiotemporal features are associated with this yield transition and affect
macroscopic material responses such as ductile versus brittle behavior.
Unfortunately, in contrast to the case for crystalline materials, our ability
to predict and control yield in disordered solids based on their constituents
and their interactions is still limited Guazzelli and Pouliquen (2018);
Cipelletti _et al._ (2020). To build such microstructural models, we need to
identify key microscopic metrics relevant to plasticity in disordered
materials Richard _et al._ (2020). Recently, excess entropy has been explored
for this purpose Galloway _et al._ (2020); Ingebrigtsen and Tanaka (2018);
Bonnecaze _et al._ (2020). In equilibrium systems, excess entropy has been
utilized to connect viscosity with interparticle structure Dzugutov (1996);
Rosenfeld (1999); Dyre (2018); Ma _et al._ (2019). Recently in far-from-
equilibrium systems, excess entropy scaling has been shown to facilitate a
relationship between microscopic structure and dynamics Ingebrigtsen and
Tanaka (2018); Bonnecaze _et al._ (2020); Galloway _et al._ (2020). Thus,
excess entropy offers an untapped signature for plasticity and a potential
tool for modeling the mechanical response of disordered solids.

The study of rheology and particle dynamics in disordered systems has a
venerable history Argon (1979). As a result of this research, theories Richard
_et al._ (2020) have proliferated in recent decades. Two of the most
successful are Mode Coupling Theory, wherein the interplay of dynamical modes
causes the emergence of rearrangements Siebenbürger _et al._ (2009); Fuchs
(2010), and Shear Transformation Zone theory, which posits that local
configurations determine where rearrangements occur Argon (1979); Falk and
Langer (2011, 1998); Slotterback _et al._ (2012). More recently, structural
signatures for rearrangement have been revealed by machine learning approaches
Cubuk _et al._ (2017, 2015); Bapst _et al._ (2020); Ma _et al._ (2019), by
study of low-frequency excitations Chen _et al._ (2010); Yunker _et al._
(2011); Chen _et al._ (2011); Seth _et al._ (2011); Khabaz _et al._ (2020); Xu
_et al._ (2007), and via local yield stress Patinet _et al._ (2016) and near-
neighbor cage dynamics Maestro and Zaccone (2017). Despite their usefulness,
difficulties remain in applying these theories to experiments because of the
need for fitting parameters Bouchbinder and Langer (2011); Siebenbürger _et
al._ (2009) and the use of empirical relationsMaestro and Zaccone (2017) that
are difficult to measure. Moreover, these theories typically do not account
for history-dependent behavior such as material memory, which is necessary to
understand plasticity.

Generally, disordered materials contain memories, i.e., microscopic signatures
related to how the material has been processed Keim and Nagel (2011); Keim _et
al._ (2019); Mukherji _et al._ (2019); Pashine _et al._ (2019); Fiocco _et
al._ (2014); Schwen _et al._ (2020); Mungan _et al._ (2019). Memory of a
previous shearing direction, for example, can be encoded into a material’s
response; once a material is sheared sufficiently in a given direction,
continued shear in that direction requires more force than in the opposite
direction Gadala‐Maria and Acrivos (1980); Keim _et al._ (2019). In jammed
systems, recent experiments and simulations have studied formation of
directional memory at low strain amplitudes, both below and near the yield
transition; far above yield, memories are erased Keim _et al._ (2013); Paulsen
_et al._ (2014); Teich _et al._ (2020). These observations, in turn, raise
important new questions: do memories require elastic storage? Is plasticity
synonymous with erasure? How do these phenomena manifest during yield, e.g.,
in storage and loss moduli?

In this contribution, we utilize excess entropy to quantify material memory
and construct a microstructural model for disordered-material response and
energy dissipation. Experiments and simulations show that three non-
dimensional parameters govern the connections between microstructure and bulk
rheology: packing density, a normalized (non-dimensional) form of the imposed
stress, and an excess entropy (microstructure-related) ratio that quantifies
the material’s ability to retain information about its initial state. Our
results confirm that memory is stored elastically and lost plastically, and
show how yield and the ductile/brittle response emerge from knowledge about
particle configurations at the microscopic scale.

##  II Results

The experiments investigate disordered solids. The solids are colloidal
monolayers of athermal, spherical particles ($\sim$ 40,000) adsorbed at an
oil-water interface (Fig. 1a). The charged particle surfaces generate a
dipole-dipole repulsion between particles. This repulsion is strong enough to
jam the entire material, arresting particle motions. To probe the effects of
disorder, we study both mono-disperse and bi-disperse spherical particle
systems with diameters of 5.6 $\mu$m and 4.1 $\mu$m-5.6 $\mu$m, respectively.
In the bi-disperse system, crystalline domains tend to be much smaller (See
Supplementary Materials). We impose many cycles of sinusoidal stress on these
samples using a custom-made interfacial stress rheometer Keim and Arratia
(2013)that permits measurement of the bulk response of the colloidal monolayer
while simultaneously recording trajectories of individual particles (see
Methods). Cyclic stress is quasi-static, insofar as the time scale for a
completion of a rearrangement ($\sim$0.5s) is much shorter than the shortest
driving period (5s) or largest inverse strain rate (20s).

Figure 1:  Overview of structure, dynamics, and response. We characterize the
disordered solid bulk response to cyclic stress from evolving configurations
of individual constituent particles. (a) Image of $\sim$40,000 particles. Part
of the raw image is shown (left). The scale bar is 200$\mu m$. Detected
particle positions are also shown (right). For illustration, color indicates
$D^{2}_{min,C}$, which quantifies the degree to which a particle has followed
a non-affine returning trajectory (blue), or a non-affine escaping trajectory
(red). The particles in this image are experiencing yield ($\gamma_{0}\sim
15.7\%$). (b) Quantification of the fractions of escaping and returning events
versus total strain amplitude. Returning events rapidly increase near the
yield point ($\gamma_{0}\sim 3.0\%$). (c) The number of particles, $Z(r)$
within a radius, $r$ of a reference particle. The radius is expressed in units
of $a$, the average distance between neighboring particles. Vertical dashed
lines indicate the limit of the first shell of neighboring particles. Inset:
radial distribution function, $g(r)$. (d) The measured strain of the material
versus the imposed stress throughout a cycle. Both stress and strain are
averaged stroboscopically over 25 cycles. The different ellipses correspond to
separate runs at different imposed stress amplitudes. Here, the area enclosed
is a result of the lag between stress and strain, which in turn quantifies the
energy dissipated from the material.

We investigate particle rearrangements by identifying non-affine deformations
within each particle’s neighborhood Falk and Langer (1998); Keim and Arratia
(2013). The degree of non-affinity is quantified by the mean-squared
displacement after subtracting the best fit affine transformation,
$D^{2}_{min}$ (see references Falk and Langer (1998); Keim and Arratia (2013)
for more information). Within cyclically sheared disordered materials, two
types of non-affine events occur (Fig. 1a): those wherein particles return to
their original position at the end of a strain cycle but along different
paths, and those wherein particles escape their nearest neighbors and do not
return Keim and Arratia (2014); Lundberg _et al._ (2008); Möbius and
Heussinger (2014); Regev _et al._ (2013); Galloway _et al._ (2020). For
visualization we define
$D^{2}_{min,C}\equiv\pm\sqrt{(D^{2}_{min,R})^{2}+(D^{2}_{min,E})^{2}}$, where
the subscripts refer to returning (R) and escaping (E) events, respectively,
and sign corresponds to the greater $D^{2}_{min}$. Both types of events
dissipate energy Lundberg _et al._ (2008); Möbius and Heussinger (2014); Keim
and Arratia (2014); Regev _et al._ (2013, 2014); Galloway _et al._ (2020).
Returning non-affine events are known to emerge near the yield point when
elasticity begins to diminish and plasticity starts to increase Keim and
Arratia (2013); Buttinoni _et al._ (2017); escaping events arise well beyond
yield Keim and Arratia (2013) (Fig. 1b). The fraction of particles undergoing
non-affine events is $f_{d}$. By following the rearrangements, we develop
understanding about trajectory dynamics within the microstructure, and we take
steps towards our ultimate goal to relate microstructure to rheology.

To quantify structure, we characterize the inter-particle forces and particle
configurations using the radial distribution function, $g(r)$. Since the
material is jammed, the motion of each particle is arrested by its neighbors
van Hecke (2009); Liu and Nagel (1998); Behringer and Chakraborty (2018); Liu
and Nagel (2010); Das _et al._ (2020); Haxton (2012). This caging, and escape
thereof, provides another lens for the non-affine motions mentioned above;
when enough particles pass each other via small changes in the structure of
their surrounding cage, the material yields Falk and Langer (1998). For
quantitative analysis, we compute $F^{*}$, the sum of the magnitudes of inter-
particle forces acting on the average particle. Specifically:
$F^{*}=2\pi\rho\int_{0}^{r_{N}}{(-\frac{\partial u}{\partial r})g(r)}rdr$;
here $\rho$ is the number density of particles, $r_{N}$ is an upper cutoff
distance below which nearest neighbor particles are found, $u$ is the pair
potential function between any two particles, $\frac{\partial u}{\partial r}$
is the force acting between any two particles, and $g(r)$ is the sample radial
distribution function as a function of separation $r$ (Fig. 1c; Methods). To
determine $r_{N}$, we use the coordination number as a function of radial
distance, $Z(r)$ (Fig. 1c). $Z(r)$ is derived from $g(r)$ and has been studied
van Hecke (2009) and recently used Maestro and Zaccone (2017) to characterize
particle interactions and their effect on bulk materials. In our systems,
neighbor shells are well defined by broad peaks in $g(r)$ separated by troughs
(Fig. 1c-inset). The extent of the nearest neighbor shell is defined as the
radius at which $Z(r)$ begins to increase rapidly for a second time (Fig.
1c-main).

We quantify disorder using excess entropy Rosenfeld (1999), the difference
between the system’s entropy and that of its ideal gas analogue (identical
pressure, temperature, etc.). The two-body approximation of excess entropy,
$s_{2}$, is calculated from $g(r)$ using a formula given in the methods
section (Eq. 4). We calculate $s_{2}$ at discrete time points to characterize
its variation within each shear cycle (more below). Since our systems are
jammed, we interpret the below-yield system $s_{2}$ as ‘frozen in’ excess
entropy.

We seek to relate these microstructural parameters to bulk rheological
properties. Recall that as the yield transition is approached from below, the
strain will begin to lag behind the oscillatory imposed stress by a phase
angle, $\delta$. If $\delta=0[rad]$, then the material is fully elastic. If
$\delta=\pi/2[rad]$, then the material is fully viscous. In between, the
material exhibits both elasticity and plasticity; the phase angle lag
quantifies dissipation (Fig. 1d) and encodes the ratio of the loss
(plasticity) and storage (elasticity) moduli,
$G^{\prime\prime}/G^{\prime}=tan(\delta)$. We will show how
$G^{\prime\prime}/G^{\prime}$ is related to the microstructural and dynamical
quantities described above ($s_{2}$, $F^{*}$, $f_{d}$).

Figure 2:  Memory within microstructure. Microstructural anisotropy reveals
signatures of memory. Below yield, anisotropic orientation remains unchanged
regardless of shear direction. Orientation quantifies stored memory. Above
yield, anisotropic orientation reverses freely to match the direction of
shear, indicating a loss of memory. (a) Radial distribution function, g(x,y,t)
at a time corresponding to one quarter of the way through a shearing cycle. We
fit an ellipse to the first neighbor ring. This ellipse stretches and
reorients over time indicating changes in structural anisotropy of the sample.
Two elliptic fits are shown at two times, t=1.25 (—) and 1.75 [cycles] (\- -
-). (b) Orientation of the sample microstructure over time as a function of
strain amplitude. With increasing strain amplitude, the microstructure
reorients to match the stretching axis. It first reorients completely at the
yield point ($3.2\%$). (c) Elongation quantified by the ratio of ellipse major
and minor axis lengths ($m/n$) over time. Below yield, elongation oscillates
directly with the strain; above yield, elongation oscillates with twice the
frequency of strain perturbation. In b & c data are averaged strobscopically
over 25 cycles.

Next, we examine structural disorder, and its variation as a function of
applied shear. The angle-dependent radial distribution function, $g(x,y)$,
quantifies microstructural order Larson (2010); Denisov _et al._ (2015) (Fig.
2a). Crucially, a nearest-neighbor ring is observable in disordered systems
composed of interacting particles Vermant and Solomon (2005); Parsi and
Gadala‐Maria (1987). In our experiments this ring deforms throughout shear
(Fig. 2a and supplementary video), in agreement with previous observations
Parsi and Gadala‐Maria (1987); Cheng _et al._ (2011); Seth _et al._ (2011);
Ding and Mittal (2015); Maestro and Zaccone (2017); Denisov _et al._ (2015).
Throughout shear, the central ring is ellipsoidal. We can readily track the
orientation and elongation of the ellipse throughout the shear cycle (Fig.
2b&c); ellipse orientation and elongation provide a measure of the sample
anisotropy. Far above yield, as the material is sheared in one direction and
then the other, the microstructural anisotropy switches between two principal
strain axes (oriented at $\pi/4[rad]$ and $-\pi/4[rad]$, counter-clockwise
from horizontal in Fig. 2a); in this situation, microstructural anisotropy is
responsive to the direction of imposed shear (Fig. 2b). Below yield, however,
the microstructural anisotropy remains in its original orientation; shearing
is not sufficient to overcome initial ‘frozen in’ material structure. This
phenomenon is apparent from changes in ring elongation (Fig. 2c) during the
shear cycle. Note that above yield the microstructure elongates twice every
shear cycle, at frequency $2\omega$, but below yield, the microstructure
elongates only once per cycle at $\omega$.

Microstructural anisotropy reveals a memory of the last direction the material
was sheared above yield (Fig. 2). To remove internal stresses, each of our
experiments is pre-sheared well above yield ($\gamma_{0}\sim 50\%$);
nevertheless, this protocol imprints an anisotropy into the sample set by the
last shear direction. Previously it was shown that this type of material
memory is imprinted into $g(x,y)$ Parsi and Gadala‐Maria (1987); Keim _et al._
(2013); Paulsen _et al._ (2014). Here, we find that this memory imprint is
associated with the principal directions of shear (Fig. 2). Once a memory is
stored, the memory is retained as long as the material is sheared elastically.
Precisely when the material yields, all memory is lost, and the microstructure
freely switches between both orientations. Taken together, these results
indicate that materials store and express memories in the elastic regime but
lose them in the plastic regime. Furthermore, recently we showed that
orientational memory is stored most strongly within crystalline domains
wherein particle rearrangements are most intensely suppressed Teich _et al._
(2020).

Figure 3:  Entropy and material memories. Variation of entropy provides means
for predicting system response to a given strain amplitude. (a) Excess
entropy, with the mean value subtracted, follows a sinusoidal response. Below
yield, its oscillation frequency is the shear cycle frequency. At yield, the
excess entropy signal has components at both the driving frequency and twice
the driving frequency: the material is beginning to forget its initial state.
Above yield, the entropy response oscillates almost exclusively at twice the
shear cycle frequency. Black dots indicate experimental data. Red lines are
fits to equation 1 with T as the only fitting parameter. The experimental data
are averaged stroboscopically over 25 cycles. (b) Amplitudes associated with
the first and second harmonics are present within the $s_{2}$ signals. Note,
that the second and first harmonic amplitudes cross each other at the yield
point, $\gamma_{0}=3\%$, designated by the vertical dashed line (- - -).

We now use excess entropy to characterize and relate observations about
imprinted memory to the system microstructure. Above yield, we find that
structural response is independent of the direction of shear (Fig. 3a,
$\gamma_{0}=6.8\%$); when the material is sheared in either direction, the
excess entropy increases and decreases as the shear is reversed. Ostensibly,
the material cannot sustain a memory above yield, because it is continually
forced out of meta-stable states within the energy landscape. Near yield,
however, the direction of shear has an effect on structural response (Fig. 3a,
$\gamma_{0}=2.2\%$). Notice, $s_{2}$ does not increase as the material is
sheared over the second half of a sinusoidal shear cycle. Finally, below
yield, the direction of shear is important; shear in one direction produces an
increase in excess entropy, and shear in the other direction produces a
decrease (Fig. 3a, $\gamma_{0}=0.7\%$).

Figure 4:  Comparisons of imposed force, microstructural excess entropy, and
bulk rheology. a) The imposed force amplitude, $F_{0}$, normalized by the
elastic force capacity, $F^{*}$, is plotted versus the excess entropy harmonic
ratio, $s_{2,h}$ (in both mono-disperse and bi-disperse experiments). A fit of
the data suggests a parabolic relationship (p-value:3.14x$10^{-13}$, and
$r^{2}$:0.989), corroborating equation 2. b) The increase in the ratio of loss
and storage moduli, $(G^{\prime\prime}/G^{\prime})$ versus strain amplitude in
both the mono-disperse and bi-disperse experiments (same legend for mono-
disperse and bi-disperse experiments as panel a). Yield is signaled by the
rapid increase in parameter values at about $0.03$ strain amplitude. Inset:
data from simulations employing Hertzian and Lennard-Jones interaction
potentials. In both cases, markers are measured values and lines are
predictions of equation 3. c) Left and right hand sides of equation 3.
Notably, all parameters are measured. The solid diagonal line (—, slope of
1.0) represents equation 3. The slope of the best fit to the data is 0.981,
p-value:4.43x$10^{-26}$, and $r^{2}$:0.944.

As seen in figure 3b, the $s_{2}$ signals are sinusoidal. The first harmonic
($\omega$) decays and the second harmonic ($2\omega$) grows with increasing
strain amplitude. The first harmonic is dominant below yield, and the second
is dominant above yield. Therefore, the amplitude of the first harmonic of
$s_{2}$ provides quantification of a stored memory, and the amplitude of the
second harmonic characterizes the degree to which memory of the initial state
is lost. Notice, these first and second harmonic amplitudes cross each other
near the yield point.

To build a relationship between excess entropy and bulk rheology, we next
investigate the connection of $s_{2}$ to the other dynamical metrics. For this
comparison, we compute the ratio of the second to first harmonic amplitude,
which we denote as $s_{2,h}$. We can relate $s_{2,h}$ to several quantities in
our system (Fig. 4). For example, $s_{2,h}$ scales with the product of
$F^{*}/F_{0}$ and $f_{d}$ (Fig. 4a), where $F_{0}$ is the amplitude of the
prescribed shear force. This relationship between dimensionless parameters
suggests that when the imposed force on the system grows larger than $F^{*}$,
the microstructure begins to permanently change, losing stored memory. Rapid
variation of $f_{d}$ also signifies the transition. These findings build on
recent work that links excess entropy and non-affine particle dynamics
Galloway _et al._ (2020); Ingebrigtsen and Tanaka (2018). Note that the
scaling in the present case is quadratic because $f_{d}$ varies nearly
linearly with the imposed force, $F_{0}$ (see Supplemental Materials).
Finally, we find that the product of $s^{2}_{2,h}$ and $F_{0}/F^{*}$ scales
linearly with $G^{\prime\prime}/G^{\prime}$ (Fig. 4c). The scaling factor for
this linear relationship is $2\phi/\pi^{2}$; here $\phi=\pi Na^{2}/A$
quantifies the particle spatial density, $a$ is the average nearest neighbor
distance derived from the first peak of $g(r)$ (Fig. 1c: inset), and $A$ is
the total area of the observed sample or simulation.

The yield phenomenology shown in Fig. 4c depends on four dimensionless
parameters: $F_{0}/F^{*}$, $s_{2,h}$, $G^{\prime\prime}/G^{\prime}$, and the
packing density $\phi$. The ratio $F_{0}/F^{*}$ characterizes the shear force
exerted on the material relative to the force required to cause
rearrangements; when $F_{0}/F^{*}\geq 1$ plasticity is non-negligible. The
microstructural quantity $s_{2,h}$ provides a metric for whether a material’s
response is dominated (or not) by memory as it experiences oscillatory strain;
this microstructural property can be interpreted as the degree of plastic
response. Finally, a familiar ratio quantifies the bulk rheological response
of the material: ($G^{\prime\prime}/G^{\prime}$). All experimental (and
simulation) data are collapsed using these dimensionless parameters, and a
direct relationship between rheology, dynamics, and microstructure is
experimentally established in the disordered solid.

Numerical simulations complement the experiments. The simulations enable us to
vary features of the disordered system that are difficult to control
experimentally. In particular, we can test ideas regarding variation of inter-
particle potential. Moreover, unlike the experimental system, which involves a
fluid-fluid interface that gives rise to viscous drag on the particles, the
simulations offer the possibility to study the validity of our new concepts in
disordered materials without viscous drag. Thus, we have conducted shear
simulations without viscous drag and with two different inter-particle
interaction potentials: Lennard-Jones, a model for atomic glass, and Hertzian,
a model for granular systems (see Methods).

The simulations and experiments exhibit remarkably similar behaviors. Across
both the experiments and simulations, a direct and common functional
relationship between excess entropy and rheology is revealed (Fig. 4c). This
relationship does not depend on the details of particle interactions, nor the
amount of disorder. Further, since simulations do not involve a background
fluid, the importance of hydrodynamic effects is ruled out. To test the limits
of applicability of our numerical findings, we introduce varying amounts of
Brownian motion into the Lennard-Jones simulations. At high thermal
temperature, the particles rearrange due to Brownian motion in addition to
shear stress, and memory cannot be formed. However, at low thermal
temperature, the experimentally observed relationship between entropy and
rheology holds (see Fig. 4b&c). Moreover, we find that jamming is required for
the storage of memories in both simulation systems. At low packing densities,
where the system easily un-jams during shear, the relation is violated. The
wide applicability of these ideas suggests the existence of a deeper
theoretical formulation. Thus, in the remainder of this paper we outline how
our results may be derived phenomenologically (for the full derivation see the
Supplemental Materials).

To elucidate the relationship between $s_{2}$ and the material properties
($G^{\prime},G^{\prime\prime}$), we perform a simple energy balance. We start
with the harmonic behavior in $s_{2}$. In this situation, energy is balanced
in terms of accumulation, $T\Delta S_{2}$, reversible (quasi-static) energy
transfer, $F^{*}x/2$, and irreversible dissipation, $f_{d}Fx$:

$T\Delta S_{2}(t)=F^{*}x(t)/2+f_{d}F(t)x(t).$ (1)

Here $x(t)$ is the displacement of the system boundary, $F(t)$ is the imposed
shear force, and $T$ is a parameter (generally different from the thermal
temperature) that converts differences in entropy to differences in energy Bi
_et al._ (2015); Ono _et al._ (2002); Bonnecaze _et al._ (2020); Khabaz and
Bonnecaze (2021). Note that this equation would not apply in a system
dominated by thermal motion, because we do not account for changes in entropy
due to thermal fluctuations. The equation also implicitly reflects the
requirement of jamming via the $F^{*}$ term. With a single fitting parameter,
$T$, the changes in harmonic behavior in excess entropy are reproduced from
below to above yield (Fig. 3a).

The harmonic transition, associated with the excess entropy found in
experiments and simulations, is captured by the first and second terms on the
right-hand-side of equation 1. $s^{2}_{2,h}$ is the ratio of those two terms:

$s^{2}_{2,h}=f_{d}\frac{F_{0}}{F^{*}}.$ (2)

This relation describes the harmonics data remarkably well (Fig. 4a). We next
build on equation 2 by incorporating a finding of shear transformation zone
theory, namely that elastic energy builds up in the microstructure until it is
plastically released via non-affine rearrangement events Falk and Langer
(2011, 1998). Quantitatively, this concept is represented as:
$G^{\prime\prime}\propto Nf_{d}G^{\prime}$, where $N$ is the number of total
particles observed; when substituted into Eq. 2 we obtain:

$\frac{G^{\prime\prime}}{G^{\prime}}=\frac{2\phi}{\pi^{2}}\frac{F^{*}}{F_{0}}s^{2}_{2,h}.$
(3)

Note, that each parameter in this expression is measured and is generally
accessible in many systems. Across strain amplitudes, remarkable agreement is
found between $G^{\prime\prime}/G^{\prime}$ measured in experiments and
simulations, and the predictions by Eq. 3 (see Fig. 4b&c).

##  III Conclusion

ID | Type | Forces | Dispersity | Diameters | $\phi$ | $\Phi[\%]$
---|---|---|---|---|---|---
A | Experiments | Dipole-dipole | Bi-disperse |  4.1, 5.6$\mu m$ | 14.02 |  $\sim$31
B | Experiments | Dipole-dipole | Mono-disperse |  5.6$\mu m$ | 13.99 |  $\sim$35
C | Simulations | Lennard-Jones | Bi-disperse | N/A | 5.03 | N/A
D | Simulations | Hertzian | Bi-disperse | 0.84, 1.16 | 9.68 | 110
E | Simulations | Hertzian | Bi-disperse | 0.84, 1.16 | 10.12 | 120
Table 1:  A summary of the properties of the systems presented, including
variety of inter-particle force, particle dispersity, particle sizes, spatial
density of particles, $\phi$, and simple area fractions of particles, $\Phi$.
We note, particles are point particles in simulations, C; hence, diameters are
not defined in system C.

Our results demonstrate that the yield transition of jammed systems has a
configurational origin rooted in the persistence of material memory. We
investigated the responses of several jammed systems undergoing cyclic shear
deformation, incorporating aspects of STZ theory, excess entropy, and harmonic
analysis into a single framework. The analysis reveals two new dimensionless
parameters and three relations, derived phenomenologically, which connect
particle configurations to bulk rheology. Importantly, the microstructural
information needed, i.e., the radial distribution function, is available in
myriad of scattering/microscopy experiments spanning length scales and
particle types Larson (2010); thus, this analysis is accessible to
experimentalists. In the future, it should be interesting to search for
similar relations for other loading conditions, such as compression or steady
shear, and to explore a wider array of particulate systems in which the
particles are not simple spheres.

We have developed a framework to understand bulk properties of jammed
materials under shear based on microstructural information. The findings hold
potential to predict behavior of a broad range of dynamically arrested
disordered materials including foams, gels, packings of nano- and micro-scale
particles, and atomic/molecular glassy matter. Our findings, perhaps, also
shed light on some deeper questions: in particular, the nature of entropy and
the potential to use entropy ideas in far-from-equilibrium media. While
entropy formulations for non-thermal systems have found utility in modeling
disparate phenomena Shannon (1948); Perunov _et al._ (2016); Jacobson (1995),
its physical interpretation often remains mysterious. Disordered particulate
packings appear to be particularly useful for clarifying this phenomenology,
since their material structure can be interrogated with relatively simple
methods.

##  IV Methods

###  IV.1 Experiments

Using a custom built interfacial stress rheometer (ISR, SI Fig. 1), we
simultaneously measure storage and loss moduli and track particle positions in
2D dense suspensions of athermal, repulsive particles. The ISR measures
rheology by imposing force on a magnetic needle adsorbed at an interface
between oil and waterShahin (1986). A stationary wall is opposite the needle,
so that shear is imposed over a distance visible by a microscope. The
displacement of the rod is measured precisely with the microscope. With
displacement (strain) and imposed force (stress), the storage and loss moduli
are calculated Brooks _et al._ (1999); Reynaert _et al._ (2008). Additionally,
the microscope is used to image the particles ($\sim 40,000$, from wall to
needle) adsorbed at the interface. The particles include charges on their
surfaces, so they exert dipole-dipole repulsive forces on each other Aveyard
_et al._ (2000); Masschaele _et al._ (2010); Park _et al._ (2010). At the
particle densities in these experiments, these forces result in particle
jamming, which we define as full kinematic restraint on each particle by its
neighbors. In all data reported here the systems are in a sinusoidal, steady
state. In the experiments, steady state occurs after five shear cycles.
Twenty-five steady state cycles are used for calculations. For more
information about these experiments and the calculations of $D^{2}_{min}$ see
Refs.Keim and Arratia (2013, 2014, 2015).

An accessible quantity in our experiments is the two-body approximation of
excess entropy, the difference between the system’s entropy and the entropy of
an ideal gas in an equivalent state ($s_{2}\sim s_{sys.}-s_{I.G.}$).
Conveniently, this quantity is calculated from the radial distribution
function, which is available in a wide range of experiments Larson (2010). The
previously derived Baranyai and Evans (1989) formula for excess entropy is:

$s_{2}=-\pi\rho\int_{0}^{\infty}\big{\\{}g(r)ln[g(r)]-[g(r)-1]\big{\\}}rdr$
(4)

where $\rho$ is the particle number density. We implement equation 4 for each
image in our experiments individually to collectively construct an entropy
time signal, $s_{2}(t)$. For specifics of our excess entropy calculations, see
Ref. Galloway _et al._ (2020).

The network force, $F^{*}$ introduced in the paper is calculated based on
inter-particle forces within the average neighborhood of particles. To make
this measurement we estimate the average number of nearest neighbors around a
particle as:

$Z(R_{c})=2\pi\rho\int_{0}^{R_{c}}g(r)rdr$ (5)

where $R_{c}$ values are shown as the horizontal axis in Fig. 1c. We estimate
experimental inter-particle forces based on potentials measured in experiments
and molecular dynamics simulations reported in Ref. Park _et al._ (2010). An
account of our estimate is included in the Supplemental Materials.

###  IV.2 Simulations

The data points for samples C were obtained using LAMMPSPlimpton (1993). At
each strain amplitude, 10 two-dimensional ensembles of 10,000 bi-disperse
Lennard-Jones particlesWidom _et al._ (1987); Patinet _et al._ (2016) were
subjected to sinusoidal shear under periodic boundary conditions at constant
confining pressure. The period of shearing was $100\times$ that of the LJ
time-scale of the particles. Prior to shearing, the samples were dynamically
equilibrated at $1\%$ of the glass-transition temperaturePatinet _et al._
(2016). During strain-controlled shearing LAMMPS’ Nosé-Hoover thermostat was
used to maintain the samples at approximately $1\%$ of the glass-transition
temperature. After 40 cycles of shearing, the shear stress was output for
another 40 cycles for later use in the calculations of the dynamic moduli. We
find that similar calculations at $9\%$ of the glass-transition temperature
begin to violate our assumption of negligible thermal energy.

For simulation samples D and E, we use HOOMD-blue Glaser _et al._ (2015);
Anderson _et al._ (2008) to impose cyclic strain on 10 particle configurations
for each of six strain amplitudes ($1,2,3,4,5,6\%$) at constant confining
pressure. Ensembles are composed of jammed states of 50:50 bidisperse mixtures
of 10,000 Hertzian particles. Ensembles are initialized from a randomly
uniform probability distribution at a packing fraction below jamming, and
subsequently quenched under FIRE minimization Bitzek _et al._ (2006) whilst
increasing the packing fraction until the desired pressure is reached. We then
run a triangle wave shear protocol, imposing a small strain step of
$10^{-4}\%$ and minimizing under FIRE after each step, until a total of 40
cycles have been completed. We calculate dynamic moduli based on the dominant
frequencies of the resulting triangle waves.

##  V Acknowledgements

We thank Doug Durian, Andrea Liu, Rob Riggleman, Ido Regev, and Sèbastien
Kosgodagan Acharige for fruitful discussions. We especially thank Andrea Liu
for the generous contribution of computational resources for our simulations.
This work is partially funded by University of Pennsylvania’s MRSEC NSF-
DMR-1120901 and by ARO W911-NF-16-1-0290. Celia Reina further thanks NSF
Career Award, CMMI-2047506.

a

##  VI Supplemental materials

Figure 1:  Interfacial Stress Rheometer Sketch of interfacial stress rheometer
(ISR). Dense colloidal monolayer sits on water-oil interface that is bounded
by two upright parallel glass walls. A thin magnetic needle cyclically shears
the monolayer using Helmholtz coils. Accurate rheometry is obtained by
tracking needle position as a function of forcing; particle tracking is used
to characterize material microstructure.

###  VI.1 Discussion of systems

Figure 2:  Disorder increasing left to right Crystalline regions visualized
via sixfold bond orientation order, $\Psi_{6}$, measured from particle
positions. Size of crystals decreases from left to right, indicating an
increase in disorder. Colors help to indicate the lattice director
(orientation) as a guide for the eye to help discern ordered and disordered
domains. Dots with large size indicate $|\Psi_{6}|>0.9$, and small dot size
indicates $|\Psi_{6}|<0.9$. (Scale bars: $100\mu m$). a) Mono-disperse,
dipole-dipole, experimental system B. b) Bi-disperse, dipole-dipole,
experimental system A. c) Bi-disperse, Hertzian, simulation system D. d) Bi-
disperse, Lennard-Jones, simulation system C.

Here we make a few comments on the breadth of properties covered by our
systems (as summarized in table 1). First, we note that the systems studied
span a large range of disorder, ranging from crystalline regions of several
hundred particles (Fig. 2a) to merely a few (Fig. 2d). Second we note that,
Lennard-Jones potentials are attractive at long distances. So our results hold
for systems where some of the particles experience attraction, as long as the
system is jammed. Finally, our systems span a wide range of length scales.
Lennard-Jones systems are atomic scale. Dipole-dipole systems are colloidal
scale. Hertzian systems are granular scale. One further difference among our
systems is that our experiments include an intermediary fluid and an
interface, whereas our simulations do not. This serves to explore the role of
specific energy dissipation (viscous drag in the experiments) versus
unspecific dissipation (the simulations).

###  VI.2 Scaling of imposed force and non-affine events

In the main text we identify that the entropic ratio $s_{2,h}$ varies as a
quadratic in figure 4a. Our model, (main body equation 2), makes no prediction
about the form of the scaling unless we know how force amplitude, $F_{0}$, and
the fraction of particles undergoing dissipative events, $f_{d}$, scale with
each other. These values scale linearly with each other in our experiments
(Fig. 3), which results in a quadratic scaling. This is an approximation, as
at very high strain amplitudes it is expected that the fraction of non-affine
events will plateau at unity. But within the linear rheological regime studied
here, this limit is not reached.

###  VI.3  Experimental potentials

In this section we describe our method of estimating the mean inter-particle
potential of our experimental systems A and B. Sulfate latex spheres of
$D_{l}=5.6\mu m$ and $D_{s}=4.1\mu m$ are adsorbed at an interface of decane
and water. The sulfate latex groups cover the surfaces of the particles,
providing a charge. The charges and the presence of the interface cause the
particles to experience dipole-dipole repulsion with each other. The dipole-
dipole form is:

$\frac{u(r)}{k_{B}T}=a_{2}\frac{1}{r^{3}}$ (1)

where $u$ is the potential, $k_{B}$ is Boltzmann’s constant, $T$ is the
thermal temperature, $a_{2}$ is the scaling constant, and $r$ is the center to
center distance of the particles. In our bi-disperse system, the average
separation between small particles is $r_{ss}=7.53\mu m$. Separation between
large particles is $r_{ll}=8.74\mu m$.

This system is often used to study interfacial colloids; Park et al. Park _et
al._ (2010) published a study that precisely measures the form of the
interparticle potential quantitatively using Monte Carlo methods and optical
tweezers. We used particles from the same manufacturer (Invitrogen
Corporation, Carlsbad, CA) as Park et al. and followed the same particle
cleaning procedure. They report that for particles of size $D_{P}=3.1\mu m$,
the mean value of $\langle a_{2,P}\rangle=5.1\pm 2.4\times 10^{-13}m^{3}$,
where $P$ subscripts indicate Park et al.’s values.

Within our bi-disperse system, the osmotic pressure is the same between large-
large and small-small particles. Here the osmotic pressure is
$-\frac{d^{2}u(r)}{dr^{2}}$. This allows us to write:

$\frac{\langle a_{2,s}\rangle}{r_{ss}^{5}}=\frac{\langle
a_{2,l}\rangle}{r_{ll}^{5}}=\frac{\langle a_{2,p}\rangle}{r_{pp}^{5}}$ (2)

These equations are not linearly independent, so we extrapolate from our
diameter-separation information ($[D_{s},r_{ss}]$ and $[D_{l},r_{ll}]$) to
determine $r_{pp}$ using $D_{P}$. We find $r_{pp}=6.72\mu m$, $\langle
a_{2,s}\rangle=9.0\times 10^{-13}m^{3}$, and $\langle a_{2,l}\rangle=1.8\times
10^{-12}m^{3}$.

Figure 3:  Imposed force scales with fraction of dissipative events. Within
the linear rheology regime studied, the fraction of particles undergoing non-
affine, disspative events scales linearly with the imposed force on the
system. The dashed lines (- - -) are added to guide the eye.

For the main body of the text we use the average weighted by particle numbers
of $\langle a_{2,s}\rangle$ and $\langle a_{2,l}\rangle$ for the entire
suspension as $\langle a_{2}\rangle$. Forces are calculated as
$F_{el.}(r)=-\frac{du}{dr}=\frac{3\langle a_{2}\rangle k_{B}T}{r^{4}}$.

###  VI.4 Phenomenological Derivation

Briefly, as discussed in the main text, equation 1 quantifies energy
accumulated in the system, $TdS$ on the left hand side and the contributions
from reversibly transferred and dissipated energy on the right hand side. Our
experiments and simulations indicate that:

$T\Delta S=F^{*}x/2+f_{d}Fx.$ (3)

Here $T$ is a fit constant that converts changes in entropy to changes in
energy. $S$ is the entropy of the entire system. $F^{*}$ is a property of the
material that quantifies how jammed the particles are via the total force
experienced by a typical particle with its neighbors. It is calculated as:
$F^{*}=\rho\int{\int{(-\frac{\partial u}{\partial r})g(x,y)}}dxdy$, where
$\rho$ is the number density of particles, $u$ is the inter-particle
potential, and g(x,y) is the radial distribution function. Here $f_{d}$ is the
fraction of particles undergoing dissipative events detected via non-affine
rearrangements. Non-affine events are detected via $D^{2}_{min}$. See Refs
Keim and Arratia (2013, 2014, 2015) for details on this calculation.
Specifically, $f_{d}=N_{d}/N$, where $N_{d}$ is the number of particles
experiencing non-affine events and $N$ is the number of total particles
observed. The proscribed shear force and resultant displacement of the
shearing surface are $F$ and $x$ respectively. We define the $\Delta$
operation as the difference between entropy at time $t$ and the average
entropy over an entire cycle of shear: $\Delta S=S(t)-\overline{S(t)}$.

We next summarize the specifics of our systems: most notably oscillatory shear
and excess entropy. To apply equation 3 to the oscillatory shear cases
considered in this paper, substitute in the time signals for shear surface
displacement ($x(t)=x_{0}sin(\omega t+\delta)$) and force
($F(t)=F_{0}sin(\omega t)$) on the right-hand-side. On the left-hand-side,
multiply by $Nk_{B}/Nk_{B}$:

$\displaystyle Nk_{B}T\Delta(\frac{S}{Nk_{B}})=$ (4)
$\displaystyle\frac{F^{*}x_{0}}{2}sin(\omega t+\delta)$
$\displaystyle+\frac{N_{d}}{N}F_{0}sin(\omega t)x_{0}sin(\omega t+\delta)$

where $\omega$ is the frequency of the imposed force and $\delta$ is the
resulting time lag between the imposed force and the resulting displacement.
$\delta$ is an important physical parameter in rheology; it helps us to
distinguish between solids, fluids, and everything in between. A fully elastic
material has a $\delta=0[rad]$; stress and strain are in phase as is seen from
Hooke’s law. A fully viscous material has a $\delta=\pi/2[rad]$; stress and
strain are fully out of phase as is seen from Newton’s law of viscosity Larson
(2010).

In our experiments, changes in pressure are negligible. Therefore, changes in
absolute entropy are approximately the same as those for excess entropy
($ds_{2}\sim ds-ds_{I.G.}\sim ds_{total}$); the ideal gas entropy is not
expected to change. Notice, entropy has changed to lower case ’s’ to represent
quantities that are normalized by $N$ and in units of $k_{B}$, which is
convention. In simulations, entropy harmonics are directly calculated on
$TS=E+PV$ because pressure, $P$, volume, $V$, and energy, $E$, are accessible
Ono _et al._ (2002); Sciortino _et al._ (1999); Bonnecaze _et al._ (2020).
Additionally, here we implement the product-to-sum trigonometric identity
($sin(u)sin(v)=(1/2)[cos(u-v)-cos(u+v)]$). Reorganizing gives:

$\displaystyle\Delta s_{2}=$ (5)
$\displaystyle\frac{F^{*}x_{0}}{2Nk_{B}T}sin(\omega t+\delta)$
$\displaystyle+\frac{N_{d}F_{0}x_{0}}{2N^{2}k_{B}T}\\{cos(\delta)-cos(2\omega
t+\delta)\\}.$

Equation 5 describes the evolution of a jammed system as it undergoes
oscillatory shear and is fully non-dimensional. It is now apparent that the
second term on the right-hand side (with $N_{d}$) has the second harmonic of
the forcing frequency $2\omega$; this relation reproduces the frequency shift
of the entropy signals in our simulations and experiments (main body Fig. 3a).
The appearance of the second harmonic in the entropy signal captures well the
transition to plasticity.

We investigate the yield transition further by taking the ratio of the first
and second harmonics within frequency domain of $s_{2}$,
($s_{2,h}\equiv\frac{FFT_{s_{2}}(2\omega)}{FFT_{s_{2}}(\omega)}$), which
follows from equation 4 as:

$s^{2}_{2,h}=\frac{N_{d}}{N}\frac{F_{0}}{F^{*}}.$ (6)

Equation 6 is visualized in figure 4a of the main text. This scaling is
quadratic because $N_{d}$ and $F_{0}$ scale linearly with each other (Fig. 3).
The square of $s_{2,h}$ in equation 6, is included so that linear
relationships are retained throughout.

From here we revisit an idea posited by Falk and Langer (Ref.Falk and Langer
(1998)): relaxation events are due to a local buildups of elastic energy that
suddenly release (i.e. $G^{\prime\prime}\propto N_{d}G^{\prime}$). Recently
quantified for above yield cases in Ref.Keim and Arratia (2014) and here
expanded to below yield, $G^{\prime\prime}=\frac{2a^{2}}{\pi
A}N_{d}G^{\prime}$, where $a$ is the first peak distance of g(r) and $A$ is
the area of observation. Substituting this equation into equation 6 for
$N_{d}$ gives:

$\frac{G^{\prime\prime}}{G^{\prime}}=\frac{2Na^{2}}{\pi
A}\frac{F^{*}}{F_{0}}s^{2}_{2,h}=\frac{2\phi}{\pi^{2}}\frac{F^{*}}{F_{0}}s^{2}_{2,h}$
(7)

which allows us to relate the bulk material response directly to measurable
microstructural properties without the use of fitting parameters. Equation 7
is visualized in Fig. 4b&c. $\phi$ quantifies particle density as $\phi=\pi
Na^{2}/A$, which implicitly takes $a$ as an effective particle diameter. This
relation reveals that the yielding transition of jammed materials is specified
by four dimensionless groups based on imposed force, particle density, a
memory based dimensionless entropy, and the bulk response.

## References

  * Nagel (2017) S. R. Nagel, _Rev. Mod. Phys._ , 2017, 89, 025002.
  * Ioannidou _et al._ (2016) K. Ioannidou, K. J. Krakowiak, M. Bauchy, C. G. Hoover, E. Masoero, S. Yip, F.-J. Ulm, P. Levitz, R. J.-M. Pellenq and E. Del Gado, _PNAS_ , 2016, 113, 2029–2034.
  * Jerolmack and Daniels (2019) D. J. Jerolmack and K. E. Daniels, _Nat. Rev. Phys._ , 2019, 1, 716–730.
  * Nie _et al._ (2020) S. Nie, Q. Jiang, L. Cui and C. Zhang, _Sediment. Geol._ , 2020, 397, 105570.
  * Larson (2010) R. Larson, _The structure and rheology of complex fluids_ , Oxford University Press, 2010.
  * Guazzelli and Pouliquen (2018) E. Guazzelli and O. Pouliquen, _J. Fluid Mech._ , 2018, 852, P1.
  * Chen _et al._ (2010) D. T. Chen, Q. Wen, P. A. Janmey, J. C. Crocker and A. G. Yodh, _Annu. Rev. of Condens. Matter Phys._ , 2010, 1, 301–322.
  * Buttinoni _et al._ (2017) I. Buttinoni, M. Steinacher, H. T. Spanke, J. Pokki, S. Bahmann, B. Nelson, G. Foffi and L. Isa, _Phys. Rev. E_ , 2017, 95, 012610.
  * Pham _et al._ (2008) K. N. Pham, G. Petekidis, D. Vlassopoulos, S. U. Egelhaaf, W. C. K. Poon and P. N. Pusey, _J. Rheol._ , 2008, 52, 649–676.
  * Cipelletti _et al._ (2020) L. Cipelletti, K. Martens and L. Ramos, _Soft Matter_ , 2020, 16, 82–93.
  * Richard _et al._ (2020) D. Richard, M. Ozawa, S. Patinet, E. Stanifer, B. Shang, S. A. Ridout, B. Xu, G. Zhang, P. K. Morse, J. L. Barrat, L. Berthier, M. L. Falk, P. Guan, A. J. Liu, K. Martens, S. Sastry, D. Vandembroucq, E. Lerner and M. L. Manning, _Predicting plasticity in disordered solids from structural indicators_ , 2020\.
  * Galloway _et al._ (2020) K. L. Galloway, X. Ma, N. C. Keim, D. J. Jerolmack, A. G. Yodh and P. E. Arratia, _PNAS_ , 2020.
  * Ingebrigtsen and Tanaka (2018) T. S. Ingebrigtsen and H. Tanaka, _PNAS_ , 2018, 115, 87–92.
  * Bonnecaze _et al._ (2020) R. T. Bonnecaze, F. Khabaz, L. Mohan and M. Cloitre, _J. Rheol._ , 2020, 64, 423–431.
  * Dzugutov (1996) M. Dzugutov, _Nature_ , 1996, 381, 137–139.
  * Rosenfeld (1999) Y. Rosenfeld, _J. Phys. Condens. Matter_ , 1999, 11, 5415–5427.
  * Dyre (2018) J. C. Dyre, _J. Chem. Phys._ , 2018, 149, 210901.
  * Ma _et al._ (2019) X. Ma, J. Liu, Y. Zhang, P. Habdas and A. G. Yodh, _The Journal of Chemical Physics_ , 2019, 150, 144907.
  * Argon (1979) A. Argon, _Acta Metal._ , 1979, 27, 47 – 58.
  * Siebenbürger _et al._ (2009) M. Siebenbürger, M. Fuchs, H. Winter and M. Ballauff, _J. Rheol._ , 2009, 53, 707–726.
  * Fuchs (2010) M. Fuchs, in _Nonlinear Rheological Properties of Dense Colloidal Dispersions Close to a Glass Transition Under Steady Shear_ , ed. M. Cloitre, Springer Berlin Heidelberg, Berlin, Heidelberg, 2010, pp. 55–115.
  * Falk and Langer (2011) M. L. Falk and J. Langer, _Annu. Rev. Condens. Matter Phys._ , 2011, 2, 353–373.
  * Falk and Langer (1998) M. L. Falk and J. S. Langer, _Phys. Rev. E_ , 1998, 57, 7192–7205.
  * Slotterback _et al._ (2012) S. Slotterback, M. Mailman, K. Ronaszegi, M. van Hecke, M. Girvan and W. Losert, _Phys. Rev. E_ , 2012, 85, 021309.
  * Cubuk _et al._ (2017) E. D. Cubuk, R. J. S. Ivancic, S. S. Schoenholz, D. J. Strickland, A. Basu, Z. S. Davidson, J. Fontaine, J. L. Hor, Y.-R. Huang, Y. Jiang, N. C. Keim, K. D. Koshigan, J. A. Lefever, T. Liu, X.-G. Ma, D. J. Magagnosc, E. Morrow, C. P. Ortiz, J. M. Rieser, A. Shavit, T. Still, Y. Xu, Y. Zhang, K. N. Nordstrom, P. E. Arratia, R. W. Carpick, D. J. Durian, Z. Fakhraai, D. J. Jerolmack, D. Lee, J. Li, R. Riggleman, K. T. Turner, A. G. Yodh, D. S. Gianola and A. J. Liu, _Science_ , 2017, 358, 1033–1037.
  * Cubuk _et al._ (2015) E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rottler, D. J. Durian, E. Kaxiras and A. J. Liu, _Phys. Rev. Lett._ , 2015, 114, 108001.
  * Bapst _et al._ (2020) V. Bapst, T. Keck, A. Grabska-Barwińska, C. Donner, E. D. Cubuk, S. S. Schoenholz, A. Obika, A. W. R. Nelson, T. Back, D. Hassabis and P. Kohli, _Nat. Phys._ , 2020, 16, 448–454.
  * Ma _et al._ (2019) X. Ma, Z. S. Davidson, T. Still, R. J. S. Ivancic, S. S. Schoenholz, A. J. Liu and A. G. Yodh, _Phys. Rev. Lett._ , 2019, 122, 028001.
  * Chen _et al._ (2010) K. Chen, W. G. Ellenbroek, Z. Zhang, D. T. N. Chen, P. J. Yunker, S. Henkes, C. Brito, O. Dauchot, W. van Saarloos, A. J. Liu and A. G. Yodh, _Phys. Rev. Lett._ , 2010, 105, 025501.
  * Yunker _et al._ (2011) P. J. Yunker, K. Chen, Z. Zhang and A. G. Yodh, _Phys. Rev. Lett._ , 2011, 106, 225503.
  * Chen _et al._ (2011) K. Chen, M. L. Manning, P. J. Yunker, W. G. Ellenbroek, Z. Zhang, A. J. Liu and A. G. Yodh, _Phys. Rev. Lett._ , 2011, 107, 108301.
  * Seth _et al._ (2011) J. R. Seth, L. Mohan, C. Locatelli-Champagne, M. Cloitre and R. T. Bonnecaze, _Nat. Mater._ , 2011, 10, 838–843.
  * Khabaz _et al._ (2020) F. Khabaz, M. Cloitre and R. T. Bonnecaze, _J. Rheol._ , 2020, 64, 459–468.
  * Xu _et al._ (2007) N. Xu, M. Wyart, A. J. Liu and S. R. Nagel, _Phys. Rev. Lett._ , 2007, 98, 175502.
  * Patinet _et al._ (2016) S. Patinet, D. Vandembroucq and M. L. Falk, _Phys. Rev. Lett._ , 2016, 117, 045501.
  * Maestro and Zaccone (2017) A. Maestro and A. Zaccone, _Nanoscale_ , 2017, 9, 18343–18351.
  * Bouchbinder and Langer (2011) E. Bouchbinder and J. S. Langer, _Phys. Rev. E_ , 2011, 83, 061503\.
  * Keim and Nagel (2011) N. C. Keim and S. R. Nagel, _Phys. Rev. Lett._ , 2011, 107, 010603\.
  * Keim _et al._ (2019) N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry and S. R. Nagel, _Rev. Mod. Phys._ , 2019, 91, 035002.
  * Mukherji _et al._ (2019) S. Mukherji, N. Kandula, A. Sood and R. Ganapathy, _Phys. Rev. Lett._ , 2019, 122, 158001.
  * Pashine _et al._ (2019) N. Pashine, D. Hexner, A. J. Liu and S. R. Nagel, _Sci. Adv._ , 2019, 5, 1–7.
  * Fiocco _et al._ (2014) D. Fiocco, G. Foffi and S. Sastry, _Phys. Rev. Lett._ , 2014, 112, 025702\.
  * Schwen _et al._ (2020) E. M. Schwen, M. Ramaswamy, C.-M. Cheng, L. Jan and I. Cohen, _Soft Matter_ , 2020, 16, 3746–3752.
  * Mungan _et al._ (2019) M. Mungan, S. Sastry, K. Dahmen and I. Regev, _Phys. Rev. Lett._ , 2019, 123, 178002.
  * Gadala‐Maria and Acrivos (1980) F. Gadala‐Maria and A. Acrivos, _J. Rheol._ , 1980, 24, 799–814.
  * Keim _et al._ (2013) N. C. Keim, J. D. Paulsen and S. R. Nagel, _Phys. Rev. E_ , 2013, 88, 032306.
  * Paulsen _et al._ (2014) J. D. Paulsen, N. C. Keim and S. R. Nagel, _Phys. Rev. Lett._ , 2014, 113, 068301.
  * Teich _et al._ (2020) E. G. Teich, K. L. Galloway, P. E. Arratia and D. S. Bassett, 2020.
  * Keim and Arratia (2013) N. C. Keim and P. E. Arratia, _Soft Matter_ , 2013, 9, 6222–6225.
  * Keim and Arratia (2014) N. C. Keim and P. E. Arratia, _Phys. Rev. Lett._ , 2014, 112, 028302\.
  * Lundberg _et al._ (2008) M. Lundberg, K. Krishan, N. Xu, C. S. O’Hern and M. Dennin, _Phys. Rev. E_ , 2008, 77, 041505.
  * Möbius and Heussinger (2014) R. Möbius and C. Heussinger, _Soft Matter_ , 2014, 10, 4806–4812.
  * Regev _et al._ (2013) I. Regev, T. Lookman and C. Reichhardt, _Phys. Rev. E_ , 2013, 88, 062401\.
  * Galloway _et al._ (2020) K. L. Galloway, D. J. Jerolmack and P. E. Arratia, _Soft Matter_ , 2020, 16, 4373–4382.
  * Regev _et al._ (2014) I. Regev, J. Weber, C. Reichhardt, K. Dahmen and T. Lookman, _Nat. Commun._ , 2014, 6, 8805.
  * van Hecke (2009) M. van Hecke, _J. Phys. Condens._ , 2009, 22, 033101.
  * Liu and Nagel (1998) A. Liu and S. Nagel, _Nature_ , 1998, 396, 21–22.
  * Behringer and Chakraborty (2018) R. Behringer and B. Chakraborty, _Rep. Prog. Phys_ , 2018, 82, 012601\.
  * Liu and Nagel (2010) A. J. Liu and S. R. Nagel, _Annu. Rev. Condens. Matter Phys._ , 2010, 1, 347–369.
  * Das _et al._ (2020) P. Das, H. A. Vinutha and S. Sastry, _PNAS_ , 2020, 117, 10203–10209.
  * Haxton (2012) T. K. Haxton, _Phys. Rev. E_ , 2012, 85, 011503.
  * Denisov _et al._ (2015) D. V. Denisov, M. T. Dang, B. Struth, A. Zaccone and G. H. Wegdam, _Sci. Rep._ , 2015, 5, 14359.
  * Vermant and Solomon (2005) J. Vermant and M. J. Solomon, _J. Phys. Condens._ , 2005, 17, R187–R216.
  * Parsi and Gadala‐Maria (1987) F. Parsi and F. Gadala‐Maria, _J. Rheol._ , 1987, 31, 725–732.
  * Cheng _et al._ (2011) X. Cheng, J. H. McCoy, J. N. Israelachvili and I. Cohen, _Science_ , 2011, 333, 1276–1279.
  * Ding and Mittal (2015) Y. Ding and J. Mittal, _Soft Matter_ , 2015, 11, 5274–5281.
  * Bi _et al._ (2015) D. Bi, S. Henkes, K. E. Daniels and B. Chakraborty, _Annu. Rev. of Condens. Matter Phys._ , 2015, 6, 63–83.
  * Ono _et al._ (2002) I. K. Ono, C. S. O’Hern, D. J. Durian, S. A. Langer, A. J. Liu and S. R. Nagel, _Phys. Rev. Lett._ , 2002, 89, 095703.
  * Khabaz and Bonnecaze (2021) F. Khabaz and R. T. Bonnecaze, _Physics of Fluids_ , 2021, 33, 013315\.
  * Shannon (1948) C. E. Shannon, _The Bell System Technical Journal_ , 1948, 27, 379–423.
  * Perunov _et al._ (2016) N. Perunov, R. A. Marsland and J. L. England, _Phys. Rev. X_ , 2016, 6, 021036.
  * Jacobson (1995) T. Jacobson, _Phys. Rev. Lett._ , 1995, 75, 1260–1263.
  * Shahin (1986) G. Shahin, _Ph.D. thesis_ , University of Pennsylvania, 1986.
  * Brooks _et al._ (1999) C. F. Brooks, G. G. Fuller, C. W. Frank and C. R. Robertson, _Langmuir_ , 1999, 15, 2450–2459.
  * Reynaert _et al._ (2008) S. Reynaert, C. F. Brooks, P. Moldenaers, J. Vermant and G. G. Fuller, _J. Rheol._ , 2008, 52, 261–285.
  * Aveyard _et al._ (2000) R. Aveyard, J. H. Clint, D. Nees and V. N. Paunov, _Langmuir_ , 2000, 16, 1969–1979.
  * Masschaele _et al._ (2010) K. Masschaele, B. J. Park, E. M. Furst, J. Fransaer and J. Vermant, _Phys. Rev. Lett._ , 2010, 105, 048303.
  * Park _et al._ (2010) B. J. Park, J. Vermant and E. M. Furst, _Soft Matter_ , 2010, 6, 5327–5333.
  * Keim and Arratia (2015) N. C. Keim and P. E. Arratia, _Soft Matter_ , 2015, 11, 1539–1546.
  * Baranyai and Evans (1989) A. Baranyai and D. J. Evans, _Phys. Rev. A_ , 1989, 40, 3817–3822.
  * Plimpton (1993) S. Plimpton, _Fast parallel algorithms for short-range molecular dynamics_ , Sandia national labs., albuquerque, nm (united states) technical report, 1993.
  * Widom _et al._ (1987) M. Widom, K. J. Strandburg and R. H. Swendsen, _Phys. Rev. Lett._ , 1987, 58, 706.
  * Glaser _et al._ (2015) J. Glaser, T. D. Nguyen, J. A. Anderson, P. Lui, F. Spiga, J. A. Millan, D. C. Morse and S. C. Glotzer, _Comput. Phys. Commun._ , 2015, 192, 97 – 107.
  * Anderson _et al._ (2008) J. A. Anderson, C. D. Lorenz and A. Travesset, _J. Comput. Phys._ , 2008, 227, 5342–5359.
  * Bitzek _et al._ (2006) E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, _Phys. Rev. Lett._ , 2006, 97, 170201.
  * Sciortino _et al._ (1999) F. Sciortino, W. Kob and P. Tartaglia, _Phys. Rev. Lett._ , 1999, 83, 3214–3217.


# Dynamics and entanglement in quantum and quantum-classical systems:
lessons for gravity

Viqar Husain<EMAIL_ADDRESS>Department of Mathematics and Statistics,
University of New Brunswick, Fredericton, NB, Canada E3B 5A3  Irfan Javed
<EMAIL_ADDRESS>Department of Mathematics and Statistics, University of New
Brunswick, Fredericton, NB, Canada E3B 5A3  Suprit Singh<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology Delhi, Hauz Khas, New
Delhi, India 110016

###### Abstract

Motivated by quantum gravity, semi-classical theory, and quantum theory on
curved spacetimes, we study the system of an oscillator coupled to two
spin-1/2 particles. This model provides a prototype for comparing three types
of dynamics: the full quantum theory, the classical oscillator with spin
backreaction, and spins propagating on a fixed oscillator background. From
nonperturbative calculations of oscillator and entanglement entropy dynamics,
we find that entangled tripartite states produce novel oscillator
trajectories, and that the three systems give equivalent dynamics for
sufficiently weak oscillator-spin couplings, but deviate significantly for
intermediate couplings. These results suggest that semiclassical dynamics with
back reaction does not provide a suitable intermediate regime between quantum
gravity and quantum theory on curved spacetime.

A quantum theory of gravity (QG) is expected to provide a unification of
gravity with the other forces of nature (for a recent review, see e.g. Carlip
et al. (2017)). The literature is abound with attempts to quantize gravity or
simplified models of it Misner (1969) with no clear consensus so far on the
approach to a final theory. If QG turns out to be a conventional quantum
theory, it will be a system with a Hilbert space ${\cal H}={\cal
H}_{\rm{gravity}}\otimes{\cal H}_{\rm{matter}}$. The matter component is in
general a “multipartite” system representing several species of matter. This
means that quantum states can have matter-gravity entanglement, and the
corresponding entanglement entropy would be an evolving observable.

If a QG theory were available, there would be several questions to pose. The
Universe we observe is well described by quantum fields on either a background
of an expanding cosmology on large scales or a flat spacetime on smaller
scales. One of the important questions is how such an approximation emerges
dynamically from quantum gravity Kiefer (1994); Padmanabhan (2019). In between
quantum gravity and quantum fields on a classical background spacetime, there
is the intermediate regime of classical gravity coupled to quantum matter with
backreaction. A proposal for this intermediate regime is the much studied
semiclassical Einstein equation Singh and Padmanabhan (1989); Brout et al.
(1995); Anderson (1995):

$G_{ab}(g)=8\pi G\bra{\Psi}\hat{T}_{ab}\left(g,\hat{\phi}\right)\ket{\Psi}.$
(1)

If this equation can be properly defined and solved, it would provide an
association $(g,|\Psi\rangle)$ of a quantum state of matter with a classical
metric $g$ (viewed in the Heisenberg picture). This is a nonperturbative
hybrid classical-quantum equation; it raises many questions, such as what is
the physical interpretation of the metric corresponding to a linear or
entangled combination of matter states, and how exactly the right-hand side is
to be defined if the metric is not known explicitly Isham (1995). There are
other hybrid models of this type: the so-called Newton-Schrodinger equation
Diósi (1984); Bahrami et al. (2014), a Friedmann-Schrodinger generalization to
cosmology and related work Husain and Singh (2019); Bojowald and Ding (2020),
and linear state evolution models using generalizations of the Lindblad
equation Oppenheim (2018).

In a gravity-matter system it is of interest to study and compare three types
of dynamics. These are the full quantum evolution, a suitably defined hybrid
quantum-classical evolution with back reaction, and quantum evolution with no
back reaction on the classical system. In the weak gravity regime, the
gravitational field may be viewed as “heavy” and slowly varying, and weakly
coupled to much lighter and faster moving matter. In this regime of couplings,
it is natural to expect that gravity behaves classically. On the other hand,
in the deep QG regime, matter-gravity coupling would be strong and could
produce highly entangled states. While a study of such comparative dynamics is
technically challenging at the field theoretic level, it is relatively
accessible in simpler models of gravity, such as cosmologies coupled to matter
and non-gravitational systems.

In this paper we study this set of questions in a model that has been a
mainstay for work in atomic physics and quantum optics, the system of an
oscillator coupled to a particle with spin, known as the Jaynes-Cummings
model. We consider a slightly more general model of an oscillator coupled to
two spins-1/2 particles, together with a spin-spin coupling. In addition to
studying the fully quantum case, we utilize the model in a new way by defining
a coupled classical-quantum model with spin back reaction on the oscillator,
and another without back reaction where the two spins propagate on an
“oscillator background.” The former case resembles a Hamiltonian version of
the semi-classical Einstein equation whereas the latter may be viewed as
simple case of quantum theory on curved spacetime. We study the comparative
dynamics numerically for a variety of initial states in the quantum-quantum
(QQ) case and compare it with the dynamics in the semiclassical (SC) and
classical background (CB) cases. (A related hybrid model with different
dynamics is studied in Fratino et al. (2014).) Our main results are that the
oscillator behaves classically even for highly entangled states in the QQ case
for sufficiently weak oscillator-spin couplings, that initial product spin
states can become maximally entangled in the SC and CB cases, and that the SC
case has unusual static solutions not present in the other cases. We discuss
implications of these results for gravitational systems, and for recent
experimental proposals that dynamically generated entanglement of matter
states may provide proof that linearized gravity must be quantum Bose et al.
(2017); Marletto and Vedral (2017); Marshman et al. (2020).

Figure 1: Oscillator coupled to two spin-$\frac{1}{2}$ particles.

The system we consider is shown schematically in Fig. 1. The oscillator takes
the place of gravity, and the spins correspond to matter. The first case is
quantum-quantum, where the entire system is quantized; the second is coupled
classical oscillator-spin system with back reaction; and the third is
(quantum) spin dynamics on a fixed classical oscillator background.
(Generalizations of such models to gravity may be achieved by extending e.g.
the scalar-cosmology case discussed recently Husain and Saeed (2020).)

Quantum oscillator-spin (QQ) The Hilbert space of the model for this case is
the tensor product of the individual Hilbert spaces of the oscillator and the
two spins, $\displaystyle{\cal H}={\cal H}_{o}\otimes{\cal
H}_{\frac{1}{2}}^{(1)}\otimes{\cal H}_{\frac{1}{2}}^{(2)}$, and the
Hamiltonian is

$\displaystyle H=$
$\displaystyle\>\left(\frac{p^{2}}{2m}+\frac{1}{2}m\omega^{2}x^{2}\right)\otimes\left(I^{(1)}\otimes
I^{(2)}\right)+I\otimes\frac{\omega_{S}}{2}\left(\sigma^{(1)}_{z}\otimes
I^{(2)}+I^{(1)}\otimes\sigma_{z}^{(2)}\right)+\frac{g_{1}}{2}\left(a\otimes\sigma^{(1)}_{+}+a^{\dagger}\otimes\sigma_{-}^{(1)}\right)\otimes
I^{(2)}$ (2) $\displaystyle\>+\frac{g_{2}}{2}\left(a\otimes
I^{(1)}\otimes\sigma^{(2)}_{+}+a^{\dagger}\otimes
I^{(1)}\otimes\sigma_{-}^{(2)}\right)+\frac{\lambda}{2}I\otimes\left(\sigma^{(1)}_{+}\otimes\sigma_{-}^{(2)}+\sigma^{(1)}_{-}\otimes\sigma_{+}^{(2)}\right)\equiv
h_{o}+h_{s}+h_{os}+h_{ss},$

where $\sigma_{z},\sigma_{\pm}$ are the Pauli diagonal and ladder operators,
and $a=x\sqrt{m\omega/2}+ip/\sqrt{2m\omega}$. The first two terms in Eq. (2)
are the Hamiltonians of the noninteracting oscillator and spins ($h_{o}$ and
$h_{s}$), the second two are the interactions of the oscillator with each of
the spins with coupling constants $g_{1}/2$ and $g_{2}/2$ ($h_{os}$), and the
third is the spin-spin interaction with coupling $\lambda/2$ ($h_{ss}$).

We restrict our attention to the $4\times 4$ truncation of the oscillator
Hamiltonian and consider initial states that are linear combinations of the
ground and first excited states. This ensures that the coupled quantum
dynamics remains in the $16$-dimensional Hilbert space
$\mathcal{H}=\mathcal{H}_{o}\otimes\mathcal{H}^{(1)}_{\frac{1}{2}}\otimes\mathcal{H}^{(2)}_{\frac{1}{2}}$.
Thus, the time dependent Schrodinger equation (TDSE) in this truncation is a
set of $16$ coupled ODEs.

Semiclassical oscillator-spin (SC) In this case the oscillator is classical
with orbits in the $\mathbb{R}^{2}$ phase space with coordinates $(x,p)$, and
the spin state is given by a vector in the Hilbert space ${\cal
H}_{\frac{1}{2}}^{(1)}\otimes{\cal H}_{\frac{1}{2}}^{(2)}$. The Hamiltonian is

$\displaystyle H=$
$\displaystyle\>\left(\frac{p^{2}}{2m}+\frac{1}{2}m\omega^{2}x^{2}\right)\left(I^{(1)}\otimes
I^{(2)}\right)+\frac{\omega_{S}}{2}\left(\sigma^{(1)}_{z}\otimes
I^{(2)}+I^{(1)}\otimes\sigma_{z}^{(2)}\right)+\frac{g_{1}}{2}\left(a\sigma^{(1)}_{+}+a^{*}\sigma_{-}^{(1)}\right)\otimes
I^{(2)}$ (3)
$\displaystyle\>+I^{(1)}\otimes\frac{g_{2}}{2}\left(a\sigma^{(2)}_{+}+a^{*}\sigma_{-}^{(2)}\right)+\frac{\lambda}{2}\left(\sigma^{(1)}_{+}\otimes\sigma_{-}^{(2)}+\sigma^{(1)}_{-}\otimes\sigma_{+}^{(2)}\right)\equiv
h_{o}^{SC}+h_{s}^{SC}+h_{os}^{SC}+h_{ss}^{SC},$

where each component is defined as for the fully quantum case. However $x$,
$p$, $a$, and $a^{*}$ are now classical variables. We define the coupled
dynamics with the TDSE for the spins and the Hamilton equations for the
oscillator:

$\displaystyle i\frac{d}{dt}\ket{\Psi}$ $\displaystyle=$
$\displaystyle\left(h_{s}^{SC}+h_{os}^{SC}+h_{ss}^{SC}\right)\ket{\Psi}$ (4)
$\displaystyle\dot{q}$ $\displaystyle=$ $\displaystyle\\{q,H_{{\rm
eff}}\\},\quad\dot{p}=\\{p,H_{{\rm eff}}\\},$ (5)

where $\ket{\Psi}$ is a spin state and

$H_{\rm{eff}}(x,p)\equiv\langle\Psi|H|\Psi\rangle.$ (6)

These are a set of $6$ coupled ODEs to be solved with initial data set
$\\{x_{0},p_{0},|\psi\rangle_{0}\\}$. The quantum dynamics is unitary by
definition, and it is readily verified using the evolution equations that the
Hamiltonian $H_{{\rm eff}}$ is a constant of motion.

Spins on classical (oscillator) background (CB): This case is the simplest of
the three. We define it by fixing an classical oscillator “background”
solution $(x_{c}(t),p_{c}(t))$ and
$a_{c}=x_{c}\sqrt{m\omega/2}+ip_{c}/\sqrt{2m\omega}$ with parameters $m$ and
$\omega$, and the time dependent spin Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle\frac{\omega_{S}}{2}\left(\sigma^{(1)}_{z}\otimes
I^{(2)}+I^{(1)}\otimes\sigma_{z}^{(2)}\right)+\ \frac{g_{1}}{2}\
\left(a_{c}(t)\sigma^{(1)}_{+}+a_{c}^{*}(t)\sigma_{-}^{(1)}\right)\otimes
I^{(2)}+I^{(1)}\otimes\frac{g_{2}}{2}\left(a_{c}(t)\sigma^{(2)}_{+}+a_{c}^{*}(t)\sigma_{-}^{(2)}\right)$
(7) $\displaystyle+\ \frac{\lambda}{2}\
\left(\sigma^{(1)}_{+}\otimes\sigma_{-}^{(2)}+\sigma^{(1)}_{-}\otimes\sigma_{+}^{(2)}\right).$

Dynamics is defined solely by the TDSE of the spin state. Thus this is a
system of 4 coupled ODEs for the spin state.

Comparing dynamics All three oscillator-spins cases defined above (QQ, SC, and
SB) have dimensional parameters $m$, $\omega$, $\omega_{S}$, $g_{1}$, $g_{2}$,
and $\lambda$. The oscillator provides “fundamental” time and length scales
$1/\omega$ and $1/\sqrt{m\omega}$, respectively (with $\hbar=1$). We set these
equal to unity and measure the remaining four parameters in these units.

Comparing dynamics in the three systems is accomplished by first fixing
initial data $\\{x(0)$, $p(0),|\Psi\rangle_{s}(0)\\}$ for the SC system, and
then (i) using the same initial spin state $|\Psi\rangle_{s}(0)$ for the QQ
and CB systems, and (ii) matching initial data for the oscillator. The latter
is accomplished for the CB case by using the oscillator solution that goes
through the phase space point $(x(0),p(0))$, and for the QQ case by using the
product oscillator-spin state

$|\Phi\rangle(0)=\left(\cos\left(\frac{\theta}{2}\right)|0\rangle+\sin\left(\frac{\theta}{2}\right)e^{i\phi}|1\rangle\right)\otimes|\Psi\rangle_{s}(0),$
(8)

(where $|0\rangle$ and $|1\rangle$ are respectively the ground and first
excited states of the oscillator), and then fixing $\theta$ and $\phi$ such
that the expectation values of $\hat{x}$ and $\hat{p}$ match:

$\displaystyle x(0)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2m\omega}}\sin(\theta)\cos(\phi)$ (9)
$\displaystyle p(0)$ $\displaystyle=$
$\displaystyle-\sqrt{\frac{m\omega}{2}}\sin(\theta)\sin(\phi).$ (10)

This ensures the closest possible initial data for all the three cases so that
we can compare $(x,p)$ and $(\langle\hat{x}\rangle,\langle\hat{p}\rangle)$,
phase space trajectories, and the evolution of spin entanglement entropy and
energy in the oscillator and spin subsystems. There is also the interesting
possibility of considering the maximally entangled Greenberger–Horne–Zeilinger
(GHZ) states of the QQ system as initial states and computing the phase space
plots $(\langle\hat{x}\rangle,\langle\hat{p}\rangle)$. However, there is no
obvious comparison that can be made with dynamics in the SC and SB systems.

Figure 2: QQ system for initial GHZ state with $m=\omega=1$. The $2\times 2$
square shows phase space plots
$(\langle\hat{p}\rangle,\langle\hat{x}\rangle)$: top left is for
$\omega_{S}=g_{1}=g_{2}=\lambda=1$; top right for $\omega_{S}=g_{1}=g_{2}=1$
and $\lambda=100$; bottom left for $\omega_{S}=0.5$, $g_{1}=g_{2}=1$, and
$\lambda=100$; and bottom right for $\omega_{S}=4$, $g_{1}=g_{2}=0.1$, and
$\lambda=2000$. These illustrate the variety of oscillator phase space
dynamics; the last case is typical of the intricate phase space patterns that
arise for large spin-spin coupling. The frame on the right shows the
corresponding spin entanglement entropy as a function of time; it is constant
at $\log 2$ for cases in the second row, but shows an oscillatory pattern for
the top left case and a slow dip and rise for the top right case. (The black
dot is the initial value of $(\langle\hat{p}\rangle,\langle\hat{x}\rangle)$
fixed by the initial GHZ state; it is the same for all cases shown.)

For all cases we integrated the coupled differential equations numerically for
a variety of initial data. The method we used ensured that probability and
$H_{\rm eff}$ are conserved at least to order $10^{-8}$. Our first set of
results is for the QQ system with the GHZ initial state

$\ket{\Phi}(0)=\frac{1}{\sqrt{2}}\left(|0++\rangle+|1--\rangle\right).$ (11)

Fig. 2 shows phase space plots $(\langle\hat{p}\rangle,\
\langle\hat{x}\rangle)$ and entanglement entropy evolution for a selection of
parameters. The former provide a remarkable characterization of the oscillator
behaviour. These cannot be compared with the other two cases, where the
oscillator remains classical, but may provide unique signatures of the GHZ
state: for instance, the bottom left frame of Fig. 2 is a fuzzy ellipse that
is not centered at the origin, a feature that does not arise in the SC and CB
cases.

We computed several solutions of the three cases with comparable initial data
(as described above) with the aim of studying the parameter ranges where the
oscillator dynamics looks similar. A representative sample is shown in Fig. 3
with parameter values $m=\omega=1$, $\omega_{S}=\lambda=2$, and oscillator-
spin coupling parameters $g\equiv g_{1}=g_{2}=0.0001,0.1,1.5$. The latter are
chosen to highlight how the dynamics in phase space, spin entanglement entropy
$S_{ent}$, and energies in the spin subsystem $E_{\rm ss}=\langle
h_{s}+h_{os}+h_{ss}\rangle$ and oscillator $E_{\rm osc}=\langle H\rangle-
E_{\rm ss}$ change with oscillator-spin couplings.

The initial data for the SC and CB cases in Fig. 3 is
$\\{q(0)=0.1,p(0)=0,|\Psi(0)\rangle=|++\rangle\\}$. The corresponding initial
state for the QQ system is obtained by using Eqs. (9) and (10) with $\phi=0$
and $\sin\theta=0.1\sqrt{2}$; these are very close to the (truncated)
oscillator coherent states.

We highlight the following features evident in Fig. 3: (i) for $g=0.0001$, the
oscillator phase space trajectories and subsystem energies are
indistinguishable in all three cases, and the spin entanglement entropy
remains nearly zero; (ii) as $g$ increases to $0.1$, differences start to
appear in each of the variables plotted, and in particular, spin entanglement
entropy increases from zero, attains its maximum value of $\log 2$, and
oscillates, notably even for the SC and CB cases; (iii) spin entanglement
entropy and subsystem energy oscillations have higher frequencies for larger
$g$ values; and (iv) the SC phase space trajectory is more expansive, with the
range of the oscillator extending an order of magnitude more than that for the
QQ and CB cases. Similar features are evident for other parameter values and
initial data.

Figure 3: Phase space, spin entanglement entropy $S_{\rm ent}$, and subsystem
energies $E_{\rm osc}$ and $E_{\rm ss}$ for $m=\omega=1$ and
$\omega_{S}=\lambda=2$ for the $g$ values indicated. The initial data for the
semiclassical cases is $x=0.1$, $p=0$ (indicated with a black dot), and spin
state $\ket{++}$; the corresponding data for the fully quantum case from (8)
is $(0.99748420879\ket{0}+0.07088902028\ket{1})\ket{++}$. (The phase space
plot axes represent $\langle x\rangle$ and $\langle p\rangle$ for the QQ, and
$x,p$ for the SC and CB cases.)

While it is gratifying to note that all cases approximately agree for
sufficiently small $g$ values, point (iv) above is especially noteworthy: it
provides evidence that the backreaction SC model, which is very similar in
form and spirit to the semiclassical Einstein equation, may not provide a
reasonable transition between the QQ and CB systems.

The SC equations provide a further curious feature not present in the QQ and
CB systems: static solutions for the oscillator. These are obtained by
considering eigenstates of the spin subsystem $h_{s}+h_{os}+h_{ss}$, and
setting $\dot{x}=\dot{p}=0$. The SC equations then reduce to

$\displaystyle\dot{x}$ $\displaystyle=$
$\displaystyle\frac{p}{m}+\\{x,E(a,a^{*},g,\lambda)\\}=0$ (12)
$\displaystyle\dot{p}$ $\displaystyle=$
$\displaystyle-m\omega^{2}x+\\{p,E(a,a^{*},g,\lambda)\\}=0,$ (13)

where $E(a,a^{*},g,\lambda)$ are the corresponding eigenvalues. For
$\lambda=0$ there are particularly simple static solutions: any point on the
circle $x^{2}+p^{2}=g^{2}/2-2\omega_{s}^{2}/g^{2}$ with $m\omega=1$. Thus we
must have $g^{2}>2\omega_{s}$. The physical interpretation of the solution
with $p=0,x\neq 0$ is that the stationary spin state “holds” the stretched
spring of the oscillator. But the physical interpretation of solutions with
$x\neq 0,\ p\neq 0$ are unusual: the spring is held stretched from equilibrium
(since $x\neq 0$) and the mass is in uniform motion since $p\neq 0$! Other
static solutions are readily computed numerically.

Discussion We described in detail three versions of the dynamics of the
oscillator coupled to two spin-$1/2$ particles (QQ, SC, and CB) with a
truncation of the oscillator to a 4-level system. Our aim was to compare the
dynamics of the oscillator, spin entanglement entropy, and subsystem energies
for the same initial conditions. We highlight three results and comment on
some implications:

  * •

For sufficiently small oscillator-spin couplings, the dynamics of the three
systems is identical; this lends support to the idea that similar results
would hold for other systems, including gravity coupled to matter.

  * •

Of particular note is that the SC system gives oscillator trajectories that
are substantially different from the QQ system for larger oscillator-spin
couplings. This may be attributed to the fact that the SC equations are non-
linear in the state, unlike the QQ and CB systems. As the same (but more
consequential) nonlinearity holds for the semiclassical Einstein equation, our
results suggest that the latter does not provide the appropriate transition
between quantum gravity and quantum fields on curved spacetime.

  * •

Spin entanglement is induced in the SC and CB systems for non-zero spin-spin
coupling $\lambda$; e.g. $g=0.1$, $\lambda=2$ in Fig.3. (The CB case is
similar to entanglement generation in Floquet dynamics Claeys et al. (2018).)
The implication for gravity is similar: initial product states of matter can
get entangled thorough the semiclassical Einstein equation, or even by
propagating on a fixed but time dependent background spacetime, provided
matter is self-interacting through any local field.

  * •

In the proposed experiments Bose et al. (2017); Marletto and Vedral (2017);
Marshman et al. (2020) for detecting quantization of linearized gravity
through entanglement generation in mass states, it is posited that interaction
between masses is not action-a-distance (as it is here for spin-spin), but is
instead generated via a mediating quantum gravitational field. However, if the
masses are sufficiently close in a laboratory setting, a point interaction may
be a good approximation, and any entanglement generated through non-
gravitational quantum interactions, whether local or not, could be significant
as demonstrated in the model discussed here. (See Ref. Fragkos et al. (2022)
for related discussion.) In the final analysis, a quantum interaction is of
course necessary to generate entanglement between masses, what ever its
origin, and the spin-spin interaction in the present model is a stand-in for
that.

  * •

There is a curious case for the SC model where for the $|++\rangle$ or
$|--\rangle$ initial states, there is no entanglement induced for $g=0,\
\lambda\neq 0$. Then, increasing $g$ from zero (i.e. turning on the classical
coupling) induces entanglement; this special case is an exception to the proof
in Hall and Reginatto (2018) which covers the $\lambda=0$ case.

Our results suggest several areas for further investigation. These include
considering in the same spirit cosmological and other gravitational models
with scalar and/or spinorial fields by extending the work in Husain and Saeed
(2020), a field theoretic version of the SC model for studying back reaction
induced entanglement in gravity coupled to a scalar field in spherically
symmetric gravity where matter gravity entanglement is a potentially important
feature Husain and Terno (2010), and linear alternatives to the semiclassical
Einstein equation as discussed in Oppenheim (2018) applied to similar model
systems; the latter may address the issue of the significant difference
between the QQ and SC systems for the results presented here.

There are a few broader features of coupled classical-quantum systems where
the quantum part evolves unitarily and the classical one via Hamilton
equations. The SC equations are well defined as a dynamical system whether
physically valid or not, and the corresponding effective Hamiltonian is
conserved; the same holds for the semiclassical Einstein equation at least for
model systems Husain and Singh (2019). This raises questions for the type of
analysis presented in Maudlin et al. (2020), where the question of energy
conservation is raised. Secondly, can classical-quantum systems, the way we
have defined them here, be viewed as “ongoing measurement” of a quantum system
by a classical apparatus (here the oscillator)? Both the oscillator and the
spins undergo continuous smooth evolution, and if the oscillator is
macroscopic, it would not itself be significantly disturbed by backreaction
from the quantum system.

Acknowledgements This work was supported by the Natural Science and
Engineering Research Council of Canada. S.S. is supported in part by the Young
Faculty Incentive Fellowship from IIT Delhi. We thank Stijn De Baerdemacker,
Carlo Rovelli, Mustafa Saeed, Danny Terno, Edward Wilson-Ewing, and Nomaan X
for constructive comments on the manuscript.

## References

  * Carlip et al. (2017) S. Carlip, D.-W. Chiou, W.-T. Ni, and R. Woodard, in _One Hundred Years of General Relativity: From Genesis and Empirical Foundations to Gravitational Waves, Cosmology and Quantum Gravity (Volume 2)_ (World Scientific, Singapore, 2017), vol. Volume 2, pp. 325–347.
  * Misner (1969) C. W. Misner, Phys. Rev. Lett. 22, 1071 (1969).
  * Kiefer (1994) C. Kiefer, Lect. Notes Phys. 434, 170 (1994), eprint gr-qc/9312015.
  * Padmanabhan (2019) T. Padmanabhan, Int. J. Mod. Phys. D 29, 2030001 (2019), eprint 1909.02015.
  * Singh and Padmanabhan (1989) T. Singh and T. Padmanabhan, Annals of Physics 196, 296 (1989), ISSN 0003-4916, URL http://www.sciencedirect.com/science/article/pii/0003491689901802.
  * Brout et al. (1995) R. Brout, S. Massar, R. Parentani, S. Popescu, and P. Spindel, Phys. Rev. D pp. 1119–1133 (1995).
  * Anderson (1995) A. Anderson, Phys. Rev. Lett. 74, 621 (1995), eprint hep-th/9406182.
  * Isham (1995) C. J. Isham, in _General relativity and gravitation. Proceedings, 14th International Conference, Florence, Italy, August 6-12, 1995_ (1995), pp. 167–209, eprint gr-qc/9510063.
  * Diósi (1984) L. Diósi, Phys. Lett. A105, 199 (1984), eprint 1412.0201.
  * Bahrami et al. (2014) M. Bahrami, A. Großardt, S. Donadi, and A. Bassi, New J. Phys. 16, 115007 (2014), eprint 1407.4370.
  * Husain and Singh (2019) V. Husain and S. Singh, Phys. Rev. D 99, 086018 (2019), eprint 1811.03673.
  * Bojowald and Ding (2020) M. Bojowald and D. Ding (2020), eprint 2011.03018.
  * Oppenheim (2018) J. Oppenheim (2018), eprint 1811.03116.
  * Fratino et al. (2014) L. Fratino, A. Lampo, and H. Elze, Physica Scripta 2014, 014005 (2014).
  * Bose et al. (2017) S. Bose, A. Mazumdar, G. W. Morley, H. Ulbricht, M. Toroš, M. Paternostro, A. Geraci, P. Barker, M. S. Kim, and G. Milburn, Phys. Rev. Lett. 119, 240401 (2017), eprint 1707.06050.
  * Marletto and Vedral (2017) C. Marletto and V. Vedral, Phys. Rev. Lett. 119, 240402 (2017), eprint 1707.06036.
  * Marshman et al. (2020) R. J. Marshman, A. Mazumdar, and S. Bose, Phys. Rev. A 101, 052110 (2020), eprint 1907.01568.
  * Husain and Saeed (2020) V. Husain and M. Saeed, Phys. Rev. D 102, 124062 (2020), eprint 2007.06609.
  * Claeys et al. (2018) P. W. Claeys, S. De Baerdemacker, O. El Araby, and J.-S. Caux, Physical Review Letters 121, 080401 (2018).
  * Fragkos et al. (2022) V. Fragkos, M. Kopp, and I. Pikovski (2022), eprint 2206.00558.
  * Hall and Reginatto (2018) M. J. W. Hall and M. Reginatto, J. Phys. A 51, 085303 (2018), eprint 1707.07974.
  * Husain and Terno (2010) V. Husain and D. R. Terno, Phys. Rev. D 81, 044039 (2010), eprint 0903.1471.
  * Maudlin et al. (2020) T. Maudlin, E. Okon, and D. Sudarsky, Stud. Hist. Phil. Sci. B 69, 67 (2020), eprint 1910.06473.