11institutetext: Novosibirsk State University and BLTP, JINR, Dubna

42

Month year

# Anti-stars in the Milky Way and primordial black holes

A.D. Dolgov

###### Abstract

Astronomical data of the several recent years which present an evidence in
favour of abundant antimatter population in our Galaxy, Milky Way, are
analysed. The data include: registration of gamma-rays with energy 0.511 MeV,
which surely originate from electron-positron annihilation at rest, very large
flux of anti-helium nuclei, discovered at AMS, and 14 stars which produce
excessive gamma-rays with energies of several hundred MeV which may be
interpreted as indication that these stars consist of antimatter. Theoretical
predictions of these phenomena, made much earlier ago are described

††issue: 2

##  1 Introduction

The story of the antimatter search in the Galaxy started probably in 1968 from
the attempts initiated by B.P. Konstantinov to search for anti-comets in the
Solar System [1, 2], which were strongly criticised by Ya.B. Zeldovich despite
very good personal relations between them.

Later activity referred mostly to cosmological antimatter but not to
antimatter in our close neighbourhood, Antimatter effects in cosmology was
probably first discussed in 1971 in ref. [3], followed by [4, 5]. Antimatter
domains in the universe were studied in [6]. Reviews on the state of art with
antimatter in cosmology were done in [4, 8]. Antimatter in the Galaxy was
considered only in the last reference [8] based on the theoretical prediction
[9], which was later elaborated in [10].

##  2 Anti-evidence

###  2.1 Cosmic positrons

Observation of intense 0.511 line, a proof of abundant positron population in
the Galaxy. In the central region of the Galaxy electron–positron annihilation
proceeds at a surprisingly high rate, creating the flux [11, 12, 13]:

$\displaystyle\Phi_{511\;{\rm keV}}={1.07\pm 0.03\cdot 10^{-3}}\;{\rm
photons\;cm^{-2}\,s^{-1}}.$ (1)

The width of the line is about 3 keV. Emission mostly goes from the Galactic
bulge and at much lower level from the disk, There are several preceding works
where this phenomenon was observed, the references can be found in the above
quoted papers.

Until recently the commonly accepted explanation was that ${e^{+}}$ [14] are
created in the strong magnetic fields of pulsars but the recent results of AMS
probably exclude this mechanism. The reason is the following. According to the
AMS data the energy spectrum of antiprotons and positrons are exactly the same
both in the form and in the absolute magnitude. This feature implies that the
mechanisms of $\bar{p}$ and $e^{+}$ production are the same and since protons
cannot be produced in the magnetic fields of pulsars, it means that neither
positrons could be created by this mechanism.

###  2.2 Cosmic antinuclei

In 2018 AMS-02 announced possible observation of six $\overline{He}^{3}$ and
two $\overline{He}^{4}$ events [15, 16]. Recent measurements reveal much more
anti-events [14]. According to the data [14] the ratio of fluxes
${\overline{He}/He\sim 10^{-9}}$ is too high, if one assumes that
$\overline{He}$ is produced in the process of cosmic ray collisions. As we see
below, such secondary creation of $\overline{He}^{4}$ is negligibly weak in
comparison with theoretical expectation. A simple possibility to explain so
high value of the observed anti-flux is to assume an existence of primordial
antimatter in the Galaxy. Possibly because of that S. Ting expressed hope to
observe anti-silicon, $\overline{Si}$.

It is not excluded that the flux of anti-helium is even much higher than the
observed one, because the low energy $\overline{He}$ may escape registration
in AMS.

The expected rate of the secondary production of anti-nuclei in cosmic rays
was calculated in ref. [17]. Anti-deuterium can be created in the collisions
${\bar{p}+p}$ or ${\bar{p}+He}$, which would produce the flux of $\bar{D}$
${\sim 10^{-7}/m^{2}/s^{-1}}/$steradian/GeV/neutron, i.e. 5 orders of
magnitude below the observed flux of antiprotons. The fluxes of
${{}^{3}\bar{H}e}$ and ${{}^{4}\bar{H}e}$, which could be created in cosmic
ray collisions are respectively 4 and 8 orders of magnitude smaller than the
flux of anti-D.

After AMS announcement of observations of anti-$He^{4}$ there appeared
theoretical attempts to create anti-$He^{4}$ through dark matter annihilation,
which does not look natural. A recent review on anti-nuclei in cosmic rays can
be found in [18].

###  2.3 Antistars in the Galaxy

A possible striking discovery of antistar population in the Milky Way was
announced recently [19]. Quoting the authors: ”We identify in the catalog 14
antistar candidates not associated with any objects belonging to established
gamma-ray source classes and with a spectrum compatible with baryon-antibaryon
annihilation.”

Figure 1:  Positions and energy flux in the 100 MeV - 100 GeV range of
antistar candidates selected in 4FGL-DR2. Galactic coordinates. The background
image shows the Fermi 5-year all-sky photon counts above 1 GeV

Of course additional confirmation of this result is necessary. A possible new
way to identify antistats was suggested in ref. [20]. It was noticed there
that prior to idirect $p\bar{p}$ contact and annihilation the formation of
atomic-like excited states consisting from proton-antiproton, proton-
antinucleous (or antiproton-nucleous), and at last nucleous-antinucleous could
be formed, which are similar to positronium.

In astrophysically plausible cases formation of such quasi-atoms is possible
in the process of the interaction of neutral atmospheres or winds from
antistars with ionised interstellar gas, These atoms rapidly cascade down to
low levels prior to annihilation giving rise to a series of narrow lines which
can be associated with the hadronic annihilation gamma-ray emission. The most
significant are L (3p-2p) 1.73 keV line (yield more than 90%) from
${p\bar{p}}$ atoms, and M (4-3) 4.86 keV (yield $\sim 60$%) and L (3-2) 11.13
keV (yield about 25%) lines from $He^{4}\bar{p}$ atoms. These lines can be
probed in dedicated observations by forthcoming sensitive X-ray spectroscopic
missions XRISM and Athena and in wide-field X-ray surveys like SRG/eROSITA
all-sky survey.

The search of such objects could be facilitated by some unusual properties of
antistars, e.g. by high velocities or unusual chemical content, see below.
According to the model of antistar-in-the-Galaxy prediction [9, 10] such
peculiarities are typical for antistars and can be good signatures to search
for them.

###  2.4 Observational bounds on galactic antimatter

Possible evidence of admixture of antimatter in the Galaxy i surely beyond the
usual expectations Normally one may expect galaxies consisting purely either
of matter or antimatter. Cosmologically large domains of matter and antimatter
may be formed if C and CP symmetries are broken spontaneously. As a result the
world would be symmeric with respect to matter and antimatter.

From the data on the cosmic gamma rays one can conclude that the nearest anti-
galaxy could not be closer than at $\sim$10 Mpc [21]. Otherwise annihilation
with protons from the common intergalactic cloud consisting of matter would be
too intensive, violating existing observational limits.

The fraction of antimatter in Bullet Cluster should be below ${<3\times
10^{-6}}$ [22] based on the upper bounds to annihilation gamma-rays from
galaxy clusters. Some limits on the fraction of antimatter at large scales can
be obtained from the CMB data which excludes large isocurvature fluctuations
at ${d>10}$ Mpc, and from Big Bang Nucleosynthesis which does not allow large
chemistry fluctuations at $d>1$ Mpc.

According to the ref. [23] the analysis of the intensity of gamma rays created
by the Bondi accretion of interstellar gas to the surface of an antistar would
create the luminosity

$\displaystyle L_{\gamma}\sim 3\cdot 10^{35}(M/M_{\odot})^{2}v_{6}^{-3}.$ (2)

It allows to put a limit on the relative density of antistars in the Solar
neighbourhood: ${N_{\bar{*}}/N_{*}<4\cdot 10^{-5}}$ inside 150 pc from the
Sun.

The presented above bounds are valid, if antimatter makes the same type
objects as the observed matter. For example, compact stellar-like objects
consisting of antimatter may be abundant in the Galaxy but still escape
observations. The bounds on the density of galactic compact antistars are
rather loose, because the annihilation proceeds only on the surface of
antistars which are the objects with short mean free path of protons, as it is
analysed in the papers [24, 25, 26].

##  3 Antistar prediction

Based on the conventional approach no antimatter object is expected to be in
the Galaxy. However, it was predicted in 1993 and elaborated in 2009 that
noticeable amount of antimatter, even antistars might be present in the Galaxy
and in its halo in non-negligible amount [9, 10].

The mechanism the papers [9, 10] was originally dedicated to the formation of
primordial black holes (PBH). It predicts a large population of massive and
supermassive PBH with log-normal mass spectrum in perfect agreement with
observations. This is the only known to us mass spectrum of PBH which is
tested by observational data.

The proposed mechanism [9, 10] also allows to solve multiple problems related
to the observed black holes in the universe in all mass ranges, in particular,
it explains the origin of supermassive BHs and black holes with intermediate
masses, from ${M\sim 10^{2}M_{\odot}}$. up to ${10^{5}M_{\odot}}$, mysterious
otherwise.

As a by-product of this mechanism compact stellar type objects, which are not
massive enough to form BHs, made of matter and antimatter are predicted.

###  3.1 Predicted mass spectrum of PBH

The log-normal mass spectrum is determned by three constant parameters and has
the following very simple form:

$\displaystyle\frac{dN}{dM}=\mu^{2}\exp{[-\gamma\ln^{2}(M/M_{0})].}$ (3)

The parameters $\mu$ and $\gamma$ are determined by an unknown high energy
physics but $M_{0}$ should be equal to the cosmological horizon mass at the
QCD phase transition, so it is predicted to be
${M_{0}\approx(10-20)M_{\odot}}$ [27]. This value perfectly fits the data. In
partociular, the log-normal form of the mass spectrum of PBHs is confirmed by
the chirp mass distribution of the LIGO/Virgo events [28]. The available data
on the chirp mass distribution of the black holes in the coalescing binaries
in O1-O3 LIGO/Virgo/ runs are analysed and compared with theoretical
expectations based on the hypothesis that these black holes are primordial
with log-normal mass spectrum. The inferred best-fit mass spectrum parameters,
$M_{0}=17M_{\odot}$ and ${\gamma=0.9}$, fall within the theoretically expected
range and shows excellent agreement with observations. On the opposite, binary
black hole models based on massive binary star evolution require additional
adjustments to reproduce the observed chirp mass distribution, see figs. 2,3

Figure 2: Model distribution $F_{PBH}(<M)$ with parameters $M_{0}$ and
${\gamma}$ for two best Kolmogorov-Smirnov tests. EDF= empirical distribution
function. Figure 3: Cumulative distributions $F(<M)$ for some astrophysical
models of binary BH coalescences.

Thus, PBHs with log-normal mass spectrum perfectly agree with the data
obtained by LIGO/Virgo/Kagra, while astrophysical BHs seem to be disfavored.
This agreement presents a very strong support in favour of the theoretical
model [9, 10].

###  3.2 How reliable the galactic antimatter prediction

To summarise: why can we trust the prediction the population of antimatter in
the Milky Way? In short the answer is because the it is a by-product of the
mechanism of PBH formation which pretty well resolves multiple problems with
the properties of observed BH in the early ($Z\sim 10$) and the present day
universe, for a review see [30].

In particular, PBHs formed according to our scenario explain the peculiar
features of the sources of GWs observed by LIGO/Virgo/Kagra [31].

The existence of supermassive black holes (SMBH) discovered in all large and
some small galaxies and even in almost empty environment is explained.
Conventional models of SMBH formation demand time which is not enough by about
two orders of magnitude.

The present day universe is full of SMBH with masses in the range
${M=(10^{6}-10^{10})M_{\odot}}$ as well as the intermediate mass black holes
(IMBH), with $M=(10^{2}-10^{5})M_{\odot}$ in unexpectedly high amount.
Moreover the SMBH are abundant in the the early, ${z=5-10}$, universe. It is
problematic to explain their formation through the canonical astophysical
mechanisms but much more natural is to assume that they are primordial.

Since the predicted features of PBH in all mass ranges well agree with the
data, one may expect that the underlying mechanism of their creation indeed
operated in the early universe. This perfectly working model of PBH creation
also predicts existence of primordial antimatter in our vicinity, so it is
quite natural to expect that antimatter is indeed abundant in the Galaxy.

Theory predicts that the primordial stellar type compact objects have quite
unusual chemical content enriched with metals, It is confirmed in particular,
by observations of extremely old stars, even of a star which formally is older
than the universe. The stars looks too old because their initial chemistry is
enriched by heavy elements, which in the standard approach demand very long
time for their synthesis. Also very high velocity stars can be present in the
Galaxy and they are indeed observed.

##  4 Mechanism of anti-creation

The suggestion of PBH (and antistar) creation proposed in references [9, 10]
is based on the SUSY motivated Affleck-Dine (AD) [32] scenario of
baryogenesis. It is assumed that there exists a scalar field $\chi$ with non-
zero baryonic number, as suggested by supersymmetry. The potential of $\chi$
has, as a rule, the so called flat directions along which the potential does
not change, remaining the same as in the origin at $\chi=0$. All these
properties are generically inherent to high energy supersymmetric models.

The field $\chi$ could reach very large magnitude moving along such flat
directions either due to rising quantum fluctuations of massless fields at de
Sitter (inflationary) stage, as in the original version of the AD baryogenesis
model or, as in the present paper, because of negative effective mass squared
of $\chi$ (higgs-like effect).

Our new input to the model is is an introduction of the interaction of $\chi$
with the inflaton field, such that the mass $m_{eff}^{2}$ stays negative for
some time. Due to this new interaction, the gate to the flat direction were
closed during almost all inflationary epoch, except for a rather short period.
When the gates are closed $\chi$ remained small, sitting near the origin. When
the gates were open $\chi$ might reach large values and when the gates close
again, $\chi$ returned to the origin carrying a large baryonic number but only
inside the bubbles of cosmologically small but possibly astrophysically large
size.

Thus the universe would look like the inverted Swiss cheese. In almost all
space baryon number density is small except for relatively small bubbles with
large $B$. Since the particles which carry baryonic number, i.e, quarks were
massless, the density contrast remained practically unnoticeable. The contrast
became huge after the QCD phase transition when massless quarks turns into
heavy protons and neutrons. Then these High-B Bubbles (HBB) turns into PBH or,
if not massive enough, into compact (anti) stars. Depending upon the direction
of the $\chi$ phase rotation in the complex $\chi$ plane such primordial stars
or antistars. could be created.

When $\chi$ is large, its evolution is governed by the quartic part of the
potential, at smaller $\chi$ the quadratic part dominates and if the flat
directions of quartic and quadratic parts are different $\chi$ starts to
rotate in the complex $\chi$-plane either clock-wise or anticlock-wise. As we
see in what follows it leads either to creation of High-B-Bubble, HBB, or
High-antiB-Bubbles, H${\bar{B}}$B.

For a toy model we assume that the quartic part of the potential has the form:

$\displaystyle U_{\lambda}(\chi)=\lambda|\chi|^{4}\left(1-\cos
4\theta\right),$ (4)

while the mass term can be written as
${U_{m}=m^{2}\chi^{2}+m^{*\,2}\chi^{*\,2}}$, or in a more convenient form:

$\displaystyle
U_{m}(\chi)=m^{2}|\chi|^{2}\left[1-\cos(2\theta+2\alpha)\right],$ (5)

where ${\chi=|\chi|\exp(i\theta)}$ and ${m=|m|e^{\alpha}}$. If ${\alpha\neq
0}$, C and CP would be broken.

In GUT SUSY baryonic number is naturally non-conserved. In our toy model this
non-conservation is achieved due non-invariance of ${U(\chi)}$ with respect to
phase rotation, $\chi\rightarrow e^{i\sigma}\chi$ with a constant phase
$\sigma$.

In the conventional version of the AD-baryogenesis the field $\chi$ after
inflation was away from the origin and, when inflation was over, it started to
evolve down to the equilibrium point, ${\chi=0}$, according to the equation of
Newtonian mechanics with the Hubble friction term:

$\displaystyle\ddot{\chi}+3H\dot{\chi}+U^{\prime}(\chi)=0.$ (6)

Baryonic number of $\chi$, $B_{\chi}=\dot{\theta}|\chi|^{2}$, is analogous to
the mechanical angular momentum of rotation in the complex $\chi$-plane. It is
quite natural that this ”angular momentum” could reach a large value and hence
the baryon asymmetry, that is the ratio of the baryonic number density to the
density of the CMBR photons could be very high $\beta=N_{B}/N_{\gamma}\sim 1$,
much larger than the observed value $\beta_{obs}\approx 10^{-9}$. The
B-conserving decays of ${\chi}$ would transfer its baryonic number to that of
quarks.

Rotation could be either clockwise or anticlockwise for different bubbles with
large $\chi$ depending upon relative locations of the flat directions in
quartic and quadratic parts of the potential in the complex ${\chi}$-plane.
Large $\chi$ lives in quartic valley, but when $\chi$ drops down, it moves to
the quadratic one starting to ”rotate”. So the angular momentum or, in other
words $B$, is generated by possible different directions of the quartic and
quadratic valleys at low ${\chi}$. If CP-odd phase ${\alpha}$ is non-
vanishing, both baryonic and antibaryonic domains might be formed with
possible dominance of one of them. Both matter and antimatter objects may
exist but it is natural to expect a global dominance of one of them, so
${B_{tot}\neq 0}$.

We assume that the Affleck-Dine field $\chi$ has the potential:

$\displaystyle
U={g|\chi|^{2}(\Phi-\Phi_{1})^{2}}+\lambda_{CW}|\chi|^{4}\,\ln\left(\frac{|\chi|^{2}}{\sigma^{2}}\right)+\lambda(\chi^{4}+h.c.)+(m^{2}\chi^{2}+h.c.),$
(7)

where the first term introduced in ref. [9, 10]. is the general renormalizable
coupling of two scalar fields, $\chi$ and the inflaton $\Phi$. The logarithmic
term is is the Coleman-Weinberg potential [33], which originates from summing
up one loop corrections to the quartic potential.

CP would be broken, if the relative phase of ${\lambda}$ and $m$ is non-zero,
otherwise one can “phase rotate” $\chi$ and come to real coefficients,
eliminating thus CP-violation.

The value of $\Phi_{1}$ is chosen so that the $\Phi$ reached the value
$\Phi_{1}$ in the process of inflation but but relatively late, so that the
duration of inflation after that made about 30-40 e-foldings. When the window
to the flat direction is open, near ${\Phi=\Phi_{1}}$, the field ${\chi}$
slowly diffuses to a large value, according to quantum diffusion equation
derived by Starobinsky, generalised to a complex field ${\chi}$. When $\Phi$
turns large, it evolves classically, oscillating near the local minimum of the
potential.

If the window to flat direction, when ${\Phi\approx\Phi_{1}}$ is open only
during a short period, cosmologically small but possibly astronomically large
bubbles with high ${\beta}$ could be created, occupying a small fraction of
the universe, while the rest of the universe has normal ${{\beta\approx 6\cdot
10^{-10}}}$, created by small ${\chi}$.

The mechanism of massive PBH formation quite different from all others. The
fundament of PBH creation is build at inflation by making large isocurvature
fluctuations at relatively small scales, with practically vanishing density
perturbations.

Initial isocurvature perturbations are in chemical content of massless quarks.
Density perturbations are generated rather late after the QCD phase
transition.
density objects occupying a minor fraction of the universe volume.

The main features of the scenario can be summarised as follows:

  * •

$\rho_{\chi}\ll\rho_{\Phi}$, even inside large $\chi$ bubbles.

  * •

Bubbles with large $\chi$ occupy a small fraction of the universe volume.

  * •

When $\Phi<\Phi_{1}$ but inflation still lasts, $\chi$ is large and oscillates
fast. Hence it does not feel shallow valleys of $m^{2}\chi^{2}$. At this stage
baryon asymmetry is not generated.

  * •

Inflation ends and the oscillations of $\Phi$ heats up the universe.

  * •

Ultimately the amplitude of $\chi$ drops down, the field starts to feel
$m^{2}$-valley, and begins to rotate, generating large baryon asymmetry.

  * •

The picture is similar to the original AD-scenario in the universe.

  * •

With the chosen values of couplings and masses the density contrast between
the bubbles and the rest of the world, before the QCD phase transition, could
be rather small, at the per-cent level.

##  5 Summary

Here we summarise some properties of the mechanism [9, 10] of PBH and antistar
formation which are discussed above but not only them.

  * •

The predicted significant antimatter population of Milky Way seems to be
confirmed by observations.

  * •

The Log-normal mass spectrum of PBH is verified by the numerous data and very
well agrees with them.

  * •

PBHs formed through the mechanism [9, 10] explain the peculiar features of the
sources of GWs observed by LIGO/Virgo/Kaga.

  * •

The considered mechanism solves the numerous mysteries of ${z\sim 10}$
universe: supermassive black holes,, early created gamma-bursters and
supernovae, early bright galaxies, and evolved chemistry including a high
level of dust.

  * •

Inverted picture of the galaxy formation is proposed when a supermassive PBH
is first formed and subsequently seeds galaxy formation.

  * •

It is predicted that globular clusters (GB) and dwarf galaxies are seeded by
IMBHs, so there should be an intermediate mass BH in the center with mass
about 2000 $M_{\odot}$ in GBs and typically somewhat higher in dwarfs. An
estimate of the density of GCs is made and agrees with the data.

  * •

A large number of the recently observed intermediate mass black holes was
predicted.

  * •

An existence of supermassive black holes observed in all large and some small
galaxies and in almost empty environment is explained.

  * •

It is claimed that a large fraction of the cosmological dark matter up to 100%
can consist of PBHs.

  * •

Strange stars with unusual chemistry, enriched with metals and high velocities
are predicted and observed. Extremely old stars are predicted and discovered.
Even, an ”older than the universe” star is found; the old age is mimicked by
the unusual initial chemistry.

  * •

In refs [9, 10] the new ideas of invoking inflation and Affleck-Dime
baryogenesis for PBH creation were proposed, which were later repeated in a
number of subsequent works.

## References

  * [1] B.P. Konstantinov, M.M. Bredov, A.I. Belyaevskij, et al, Cosmic Research, 4, 66 (1968)
  * [2] B.P. Konstantinov, M.M. Bredov, S.V. Gellentskij, et al Bulletin of the Academy of Sciences of the USSR. Physical series, 33, №,11, 1820 (1969)
  * [3] F. W. Stecker, D. L. Morgan, Jr., J. Bredekamp, Possible ”Evidence for the Existence of Antimatter on a Cosmological Scale in the Universe”, Phys. Rev. Letters 27, 1469 (1971).)
  * [4] F. W. Stecker, ”Grand Unification and possible matter-antimatter domain structure in the the universe”. Tenth Texas Symposium on Relativistic Astrophysics, p. 69 (1981).
  * [5] V.M. Chechetkin, M.Yu. Khlopov, M.G. Sapozhnikov, Ya.B. Zeldovich, ”Astrophysical aspects of anti-proton interaction with He-4 (antimatter in the universe). Phys. Lett. B 118, 329 (1982).
  * [6] M.Yu. Khlopov, S.G. Rubin, A.S. Sakharov, ”Possible origin of antimatter regions in the baryon dominated universe”, Phys. Rev. D, 62, 083505 (2000).
  * [7] F. W. Stecker, ”The Matter-Antimatter Asymmetry of the Universe (keynote address for XIVth Rencontres de Blois)” arXiv:hep-ph/0207323.
  * [8] A.D. Dolgov, ”Cosmological matter antimatter asymmetry and antimatter in the universe”, keynote lecture at 14th Rencontres de Blois on Matter - Anti-matter Asymmetry • e-Print: hep-ph/0211260.
  * [9] A. Dolgov, J. Silk, ”Baryon isocurvature fluctuations at small scale and baryonic dark matter”, Phys. Rev. D 47, (1993) 4244,
  * [10] A.Dolgov, M. Kawasaki, N. Kevlishvili, ”Inhomogeneous baryogenesis, cosmic antimatter, and dark matter”, Nucl. Phys. B 807 (2009) 229,
  * [11] G. Weidenspointner et al., ”The sky distribution of positronium annihilation continuum emission measured with SPI/INTEGRAL”, Astron. Astrophys. 450, 1013 (2006); astro-ph/0601673.
  * [12] J. Knodlseder et al., ”The sky distribution of positronium annihilation continuum emission measured with SPI/INTEGRAL”, Astron. Astrophys. 441, 513 (2005); [arXiv:astro-ph/0506026].
  * [13] P. Jean et al. ”Spectral analysis of the Galactic $e^{+}e-$ \- annihilation emission”, Astron. Astrophys. 445, 579 (2006); [arXiv:astro-ph/0509298].
  * [14] S. Ting, L’Aquila Joint Astroparticle Colloquium, 10th November, 2021.
  * [15] A. Choutko, AMS-02 Collaboration, “AMS Days at La Palma, La Palma, Canary Islands, Spain,” (2018).
  * [16] S. Ting, Latest Results from the AMS Experiment on the International Space Station. Colloquium at CERN, May, 2018.
  * [17] R. Duperray, B. Baret, D. Maurin et al, ”Flux of light antimatter nuclei near Earth, induced by cosmic rays in the Galaxy and in the atmosphere”, Phys. Rev. D 71, 083013 (2005),
  * [18] P. von Doetinchem, K. Perez, T. Aramaki, et al, ”Cosmic-ray antinuclei as messengers of new physics: status and outlook for the new decade”, JCAP 08, 035 (2020): e-Print: 2002.04163 [astro-ph.HE].
  * [19] S. Dupourqué, L. Tibaldo and P. von Ballmoos, ”Constraints on the antistar fraction in the Solar System neighborhood from the 10-year Fermi Large Area Telescope gamma-ray source catalog”, Phys Rev D 103, 083016 (2021).
  * [20] A.E. Bondar, S.I. Blinnikov, A.M. Bykov, A.D. Dolgov, K.A. Postnov, ”X-ray signature of antistars in the Galaxy”, JCAP 03, 009 (2022); e-Print: 2109.12699 [astro-ph.HE]
  * [21] G. Steigman, ”Observational tests of antimatter cosmologies”, Ann. Rev. Astron. Astrophys. 14, 339 (1976).
  * [22] G. Steigman, ”When Clusters Collide: Constraints On Antimatter On The Largest Scales”, JCAP 10, 001 (2008); e-Print: 0808.1122 [astro-ph].
  * [23] P. von Ballmoos, ”Antimatter in the Universe : Constraints from Gamma-Ray Astronomy”, Hyperfine Interact. ]bf 228, 1, 91 (2014); Contribution to: LEAP2013, 91-100; e-Print: 1401.7258 [astro-ph.HE].
  * [24] C.Bambi, A.D. Dolgov, ”Antimatter in the Milky Way”, Nucl. Phys.B 784, 132 (2007); astro-ph 0702350.
  * [25] A.D. Dolgov, S.I. Blinnikov, ”Stars and Black Holes from the very Early Universe”, Phys. Rev. D 89 (2014) 2, 021301 (2014); astro-ph 1309.3395.
  * [26] S.I.Blinnikov, A.D., K.A.Postnov, ”Antimatter and antistars in the universe and in the Galaxy”, Phys. Rev. ]D 92, 023516 (2015); astro-ph 1409.5736.
  * [27] A.Dolgov, K.Postnov, ”Why the mean mass of primordial black hole distribution is close to $10M_{\odot}$”, JCAP 07, 063 (2020); e-Print: 2004.11669 .
  * [28] A.D. Dolgov, A.G. Kuranov, N.A. Mitichkin, S. Porey, K.A. Postnov, O.S. Sazhina, I.V. Simkin, ”On mass distribution of coalescing black holes,” JCAP 12, 017 (2020); e-Print: 2005.00892.
  * [29] Author Book title (Publisher, city year) page numbers
  * [30] A.D. Dolgov ”Massive and supermassive black holes in the contemporary and early Universe and problems in cosmology and astrophysics,” Usp. Fiz. Nauk 188, 121 (2018); Phys. Uspekhi. 61, 115 (2018).
  * [31] S. Blinnikov, A. Dolgov, N. Porayko, K. Postnov, ”Solving puzzles of GW150914 by primordial black holes”, JCAP 11, 036 (2016).
  * [32] I. Affleck, M. Dine, ”A new mechanism for baryogenesis”. Nuclear Physics B249, 361, (1985).
  * [33] S. Coleman, E. Weinberg ”Radiative Corrections as the Origin of Spontaneous Symmetry Breaking”. Phys. Rev. D7 1888 (1973).


# The emergence of soft-glassy mechanics in simulated foams

Amruthesh Thirumalaiswamy  Robert A. Riggleman<EMAIL_ADDRESS>John C.
Crocker<EMAIL_ADDRESS>Department of Chemical and Biomolecular
Engineering, University of Pennsylvania, Philadelphia, Pennsylvania

###### Abstract

Several seemingly different soft materials, including foams, cells, and many
complex fluids, exhibit remarkably similar rheological properties and
microscopic dynamics, termed soft glassy mechanics. Here, we show that such
behavior emerges from a simple model of a damped ripening foam, for
sufficiently weak damping. In particular, we observe intermittent avalanchey
dynamics, bubble super-diffusion, and power-law rheology that vary as the
damping factor is changed. In the limit of weak damping, the dynamics are
determined by the tortuous low-lying portions of the energy landscape, as
described in a recent study. For strong damping the viscous stresses cause the
system configuration to evolve along higher energy paths, washing out small-
scale tortuosity and producing motion with an increasingly ballistic
character. Using a microrheological approach, the linear viscoelastic response
of the model can be efficiently calculated. This resembles the power-law
rheology expected for soft glassy mechanics, but unexpectedly, is only weakly
sensitive to the damping parameter. Lastly, we study the reported memory
effect in foams after large perturbations and find that the timescale of the
memory goes to zero as the damping parameter vanishes, suggesting that the
effect is due to viscous stress relaxation rather than slow structural changes
stabilized by the energy landscape.

††preprint: APS/123-QED

##  I Introduction

Soft glassy materials (SGMs) [1, 2, 3] such as foams and emulsions exhibit
complex physical and rheological properties that continue to defy explanation.
Moreover, the similarity of soft glassy mechanics to that of living cells [4,
5, 6] and glassy materials [7] has long been noted. Previous experimental and
theoretical models have captured different aspects of such systems while
falling short of a complete physical picture. For foams, rheological
experiments have shown conflicting results — showing weak [8, 9] or no [10]
power-law frequency dependence of the dynamic shear modulus. Modeling efforts
have largely focused on the now canonical ‘bubble model’ [11, 12], but the
dynamic shear modulus of this model has not been reported. While a more recent
study did report power-law rheology [13] it used a simplified system without
damping. Further, experiments have shown memory effects [14, 15] in which a
deformed foam shows perturbed mechanics which relaxes back to the unperturbed
trend after a long time. The physical origin of this memory effect remains
poorly understood.

Here, we study the soft glassy mechanics and rheology of foams, as well as
their recovery from mechanical perturbation using a 3-D bubble model [11, 12,
16] with a simple damping law [12, 7], driven by simulated Ostwald ripening
[17]. Previous stress-strain simulations [18, 16] of a 2-D bubble model
without ripening have indicated a transition to avalanchey dynamics with
reduced applied strain rate. We look for a similar effect in our ripening foam
model by changing the damping parameter $\xi$, effectively changing the
relative rates of ripening and viscous relaxation. This however requires the
computationally expensive integration of the bubble model’s equation of motion
at low $\xi$. We find that for sufficiently low damping (or equivalently slow
ripening), the system dynamics are determined by the tortuous character of the
energy landscape, as observed in a damping-free model [13], leading to
avalanches in energy, super-diffusive bubble motion, and fractal
configuration-space paths. For stronger damping, this behavior disappears,
being replaced by a more continuous motion having a ballistic character. We
use a microrheological approach to determine the dynamic shear modulus of our
model from its intrinsic, non-thermal fluctuations [19, 20], and find that it
generically has power-law rheology resembling recent experimental measurements
[8, 9]. The rheology exponent is, unexpectedly, only a weak function of
damping, providing new insights into the origin of power-law rheology in SGMs.
Lastly, we study foam’s recovery from mechanical perturbation by randomly
scrambling the locations of bubbles in our model, finding that scrambling
leads to perturbed mechanics that slowly return to the (average) unperturbed
baseline, resembling experimental reports of mechanical memory in foams [14,
15]. The foam recovers to the baseline more quickly as the damping factor is
reduced, and does so immediately when damping is removed, indicating that the
memory effect is controlled by viscous stress relaxation, and not due to
activation between energy minima.

##  II Damped SGM Model

###  II.1 Coarsening bubble dynamics

We model a coarsening foam using the bubble model [11, 12] with a simplified
damping rule and simulated Ostwald ripening. While the bubble model has been
traditionally used to simulate foams [11, 18], it also serves as an effective
model for many other SGMs [21, 7, 13]. The constituent bubbles of foam in this
model are treated as soft-sphere particles that can overlap and interact via a
pairwise repulsive potential when overlapping:

$V(\mathbf{r}_{ij})=\begin{cases}\frac{\epsilon}{2}{\left(1-\frac{\lVert{\mathbf{r}_{ij}}\rVert}{a_{i}+a_{j}}\right)}^{2},&\text{if
$\lVert\mathbf{r}_{ij}\rVert<a_{i}+a_{j}$}\\\ 0,&\text{otherwise}.\end{cases}$
(1)

with $r_{ij}$ being the distance between two bubbles of radii $a_{i}$ and
$a_{j}$, and $V(\mathbf{r}_{ij})$ being the corresponding potential.

The positions of the bubbles interacting via the pairwise potential in Eq. 1
are evolved using an (over-damped) equation of motion [11]. Notably, we
consider a simplified version of the viscous force,
$\mathbf{F}_{i}=\xi\mathbf{v}_{i}$ [12, 7], on each bubble to reduce
computational overload while preserving relevant model physics.

$\begin{split}\xi\frac{d\mathbf{r}_{i}}{dt}=-{{\sum}}_{j}^{nn}\frac{\partial
V(r_{ij})}{\partial r_{i}}\end{split}$ (2)

with the right-hand side representing a summation over neighboring bubbles
that contribute to the force on bubble $i$. Meanwhile, the left side is the
damping force with $\xi$ being the effective viscous damping factor.

To simulate the mass exchange between bubbles due to Ostwald ripening [17],
the bubble radii are allowed to evolve while keeping total bubble volume
constant (preserving notional mass). Ripening causes larger bubbles to grow
and smaller ones to shrink over time via a pairwise mass flux. We model this
process with a flow rate that depends on the degree of overlap between
neighboring bubbles (over the overlap cross-section), along with a mean-field
flux that flows through the connected phase medium.

$\begin{split}Q_{i}=\underbrace{-\alpha_{1}\rho\sum_{j}^{nn}\left(\frac{1}{a_{i}}-\frac{1}{a_{j}}\right)A_{\text{overlap}}}_{\text{neighbor-
neighbor}}\\\
\underbrace{-\alpha_{2}\rho\left(\frac{1}{a_{i}}-\frac{1}{\left<a\right>}\right)a_{i}}_{\text{mean-
field}}\end{split}$ (3)

As indicated, the first term represents the pair-wise mass-flux between
neighboring bubbles over the cross-section of overlap and the second
represents the mean field contribution. The values for $\alpha_{1}(=0.05)$ and
$\alpha_{2}(=0.002)$ were chose akin to our previous study [13]. When a
bubble’s volume turns negative over the course of the simulation, we remove it
from the simulation box, while ensuring that the mass of the deleted bubble
and its neighbors is conserved.

To realize the overall evolution of the system, we independently evolve the
system using the equation of motion (Eq. 2) and ripening rules (Eq. 3) for all
bubbles at every time step. We use a simple explicit Euler scheme with small
$dt$ values, to numerically integrate the equations of motion (more details in
Section A.2. The use of other integrators like a second-order Runge-Kutta
numerical discretization led to similar results. While we note that the
equations of motion can be physically unstable at very high energies (when
there is a significant overlap between bubbles), we verify that such large
overlaps are not present at the energy levels presented here. Further, while
using small step-sizes within the range of numerical stability (more details
in Section A.2), we ensure the simulation has converged by cross-validating
with smaller step-sizes ($dt$). It may be noted here that smaller step-sizes
($dt$) are required for lower values of $\xi$, and thus are computationally
more expensive to simulate.

Bubble model simulations [13] when initialized randomly, eventually reach a
dynamic steady state with a characteristic bubble size distribution and
various systematic trends in properties such as total energy or mean bubbles
size, as with experiments [22]. We initialize a system of $N\sim 1000$ bubbles
at a volume fraction of $\phi=0.75$ (just above its jamming volume fraction
[13, 23]) with a Weibull bubble radius distribution,
$\textbf{P}(a)\sim(k/\lambda)(x/\lambda)^{k-1}$, where $k=1.75,\lambda=0.73$,
and let the system evolve as a function of time. This distribution is
representative of the steady state bubble size distribution that a Gaussian
initialized system reaches when evolved in the quasi-static limit ($\xi\to 0$)
[13]. Using this as a starting point for all our simulations, we model over a
range of damping factors ($\xi$) and calculate various physical quantities of
interest. It must be noted that we end our simulations when bubbles grow
considerably large leading to multiple ($>1$) overlaps between $2$ bubbles.
This happens much earlier in larger $\xi$ systems thus producing shorter
simulation trajectories overall.

(a)(a)(a)

(a)

(a)

(a)

(b)

(a)

(c)

Figure 1: (a) Traces of relative energy differences $\Delta U/{U(t)}$ (for
$\Delta t=1$ or simulation points spaced by $1$ time unit) are sensitive to
intermittent dynamics. For lower damping, $\xi\lesssim 0.01$, the relative
change in energy shows abrupt peaks characteristic of intermittent motion. (b)
Lower $\xi$ simulations show a heavy-tailed probability distribution of energy
fluctuations, typical of an avalanchey system. As the system becomes more
viscous, the energy fluctuations become more Gaussian. (c) The average system
coordination number $\langle z\rangle$ remains low for lower $\xi$
simulations, characteristic of lower energy configurations close to potential
energy minima on the landscape. $z_{C}$ is the critical coordination for
jamming, with $z_{C}=6$ in $3-D$ [21].

Fluctuations in the system’s total potential energy change significantly for
different simulated viscosities. Two distinct limits are observed, as shown in
Fig. 1. Low viscosity simulations ($\xi\leq 0.001$), produce large
fluctuations in $\Delta U(\Delta t=1)/{U(t)}$ (see Fig. 1a, indicative of
avalanchey, intermittent dynamics. These are suggestive of the system
following the ‘bumpy’ lower levels of the energy landscape. Conversely, with
higher $\xi$ values, the system no longer moves from minimum to minimum of the
underlying energy landscape but evolves in a dynamic force balance between the
larger interaction forces and viscous stresses. This allows the system to fly
over the barriers and rugged features of the energy landscape, with a higher
time average potential energy for the system (see Fig. 6a). This change in the
fluctuations is shown more clearly in the distribution of energy drops, Fig.
1b, which becomes more heavy-tailed at lower $\xi$. Further, a similar trend
can also be seen in Fig. 1c, where the average coordination number (over a
system configuration) is higher at higher viscosities, indicating that viscous
stresses are shifting the foam structure away from the minimum energy states,
and farther from jamming, defined as coordination with $\langle z\rangle\simeq
z_{c}$ [21]. Thus foam configurations formed at low damping explore lower and
more tortuous portions of the energy landscape and ones with higher damping
cruise through higher and apparently smoother portions of the energy
landscape.

(a)(a)

(a)

(a)

(a)

(b)

Figure 2: Analysis of bubble motion in 3-$N$ dimensional space and real space
shows a mixture of fractal and ballistic motion. (a) The different simulation
points in high dimensional ($3N$) space are analyzed for end-to-end distances
$\Delta R^{2}$ and contour lengths $\Delta s$ to study fractal scaling over
different length scales. Simulations with larger $\xi$ values give almost
ballistic scaling in hyperspace. However, lower $\xi$ lead to a more tortuous
trajectory characteristic of a fractal path, leading to super-diffusive
scaling. The grey dashed line is a reference with ballistic scaling $\Delta
R^{2}\sim\Delta s$ throughout. All data points above represent values pooled
over $4$ simulations and log-bin-averaged over contour distances. (b) Time and
ensemble ($4$) averaged mean-squared displacement for an ensemble of bubbles
that remain finite sized throughout our simulations, plotted for the different
$\xi$ values shows ballistic motion rolling over to a super-diffusive form for
lower $\xi$ simulations. In comparison, more viscous simulations show more
ballistic behavior over a larger range of $\tau$.

To characterize the system’s high dimensional motion over the energy
landscape, we look at the path traversed by the system through configuration
space for the range of viscosities studied. The different time points on a
simulation trajectory in configuration space are analyzed for end-to-end
distances $(\Delta R^{2})$ and path contour lengths $(\Delta s)$. This serves
as a measure of the tortuosity of the $3N$ dimensional configurational
trajectory taken by the system over time. As expected from our conclusions
above, we observe that lower viscosities yield fractal, self-similar scaling
at large lengthscales with a fractal dimension of $D_{f}\sim 2/1.38\text{\
(slope at large distances)}\simeq 1.45$ (see Fig. 2a) – capturing the
intrinsic fractal physics of the landscape [13]. Simulations with higher
damping show almost no large lengthscale fractal character, indicative of
their ability to avoid lower energy portions of the energy landscape.
Alternatively, the slight bends on Fig. 2 may be interpreted as a shift in the
lengthscale (as a function of $\xi$) over which a fractal slope would be
observed. This, however, is further evidence for the self-similar fractal
nature of the landscape and indicates that one would have to examine
considerable lag times (or configuration distances traveled by the system) to
observe soft-glassy mechanics for systems with larger damping. Here, it must
be noted that when particles shrink to zero size as a result of ripening, we
fix their positions in space, thus conserving the number of dimensions ($3N$)
used to calculate $\Delta R^{2}$. In Fig. 2b, we compute the ensemble and
time-averaged mean-squared displacement as a function of lag time ($\tau$).
These curves show a functional form that is similar to that of the $\Delta
R^{2}$ above because the mean-squared displacement is a projection of those
curves to 3-D space; the slight difference in the exponent is due to the
calculation being done on a slightly different ensemble of bubbles that remain
finite sized throughout the simulation.

Here, one may also identify a dimensionless group of interest called the
Deborah number $De$, which can be expressed as the ratio of time scales
associated with relaxation and the mode of driving—the two relevant dynamic
processes for this system. Here, that would be the damped relaxation time from
Eq. 2 ($\tau_{R}=\xi{\langle a\rangle}^{2}/\epsilon$) and timescale associated
with changing bubble radii $(a)$ imparted by the ripening process, Eq. 3
($\tau_{C}={\langle a\rangle}^{2}/\alpha_{1}$ when $\alpha_{1}>\alpha_{2}$).
This gives us a ripening Deborah number,
($De_{\alpha}=\xi\alpha_{1}/\epsilon$) which is a ratio of the relaxation
($\tau_{R}$) and coarsening ($\tau_{C}$) times (typically ranging between
$10^{-3}-10^{-6}$ for our simulations). This dimensionless group presumably
depends on the system’s volume fraction $\phi$ and its proximity to the
jamming volume fraction $\phi_{J}$ [11, 18].

This dimensionless group formalism can be a useful way to explain many
previous experimental and simulation results [12, 18, 16]. We begin by noting
that the avalanchey dynamics and intermittent rearrangements observed in our
simulations resemble previous studies of similar systems [11] driven by shear
strain instead of coarsening. Various comprehensive studies [18, 16] using
$2D$ shear strain point out a similar transition to avalanchey rearrangement
events below a certain shear strain rate. Thus, our results can be interpreted
as a transition in landscape physics as a function of $De_{\alpha}$ while
shear simulation results [16] can be explained using a corresponding shear
Deborah number, $De_{\gamma}$.

For a foam experiment, we note that the energy scale and damping factor vary
as the system evolves: $\epsilon\simeq\sigma{\langle a\rangle}^{2}$ [11] and
$\xi\propto{\langle a\rangle}$, while $\alpha_{1}$ is effectively independent
of $\langle a\rangle$. Thus experimentally, $De_{\alpha}\propto{\langle
a\rangle}^{3}$ changes for dynamically aging foam where ${\langle a\rangle}$
increases as a function of time [13, 22](see Fig. 5). This keeps pushing the
aging system away from the landscape-dominated regime, potentially explaining
the issue associated with the shifting cut-off [1], and tending to produce
behavior akin to high $\xi$ simulations.

###  II.2 Rheology of SGMs

The rheology of soft-glassy systems is typically found to be weakly frequency
dependent (solid-like), often with a power-law form, while different
experiments on foams [10, 8] yield apparently conflicting results.
Computationally, capturing low-frequency responses to applied strains can be
very expensive, making the determination of rheology difficult [13]. Here, we
provide a numerical procedure that derives its essentials from a
microrheological approach [24, 19] that computes the power spectra of the
active, fluctuating shear strain and stress from the particle motions, and
computes the dynamic shear modulus from their ratio.

We begin by noting that one can relate the stress ($\mathbf{\sigma}(t)$) and
strain ($\mathbf{\gamma}(t)$) to the creep compliance ($\mathbf{J}(t)$) using
the theory of linear response [25, 26] and the Boltzmann superposition
principle, relating them through a convolution:

$\begin{split}\mathbf{\mathbf{J}}(t)\circledast\dot{\mathbf{\sigma}}(t)=&\mathbf{\gamma}(t)\\\
\int_{-\infty}^{t}{{\mathbf{J}}(t-t^{\prime})\dot{\mathbf{\sigma}}(t^{\prime})dt^{\prime}}=&\mathbf{\gamma}(t)\end{split}$
(4)

While this basic constitutive equation represents the relation between the
macroscopic stress and strain for a linear material, we extend this formalism
to its micro-rheological version wherein each bubble/ particle can be treated
as a tracer moving in a homogeneous viscoelastic continuum (formed by all the
other bubbles) driven by active fluctuating stresses. Thus, typically one can
use the bubbles’ positional vectors describing their motion in the effective
medium to describe the local, time-dependent strain in the effective medium
[24]. Similarly, the local fluctuating active stress acting on each bubble in
the system can be computed as follows [27, 28]:

$\sigma(\mathbf{r}_{i})=-\left(\sum_{j}^{nn}\mathbf{r}_{ij}\otimes\mathbf{F}_{ij}\right)\delta(\mathbf{r}-\mathbf{r}_{i})$
(5)

where $r_{ij}$ and $F_{ij}$ represent the inter-particle displacements and
forces between particles $i$ and $j$.

Applying the above equations directly to the data would be impractical because
the $\sigma(t)$ and $\gamma(t)$ signals for each bubble are random functions
of time. Instead, we transform the equation described in Section B.1 to a
relation between the ensemble-averaged mean squared differences (MSD), the
stress, and strain. The stress MSD is calculated by considering the squared
difference between the three off-diagonal elements of the bubble-wise
symmetric tensor (see Eq. 5). Further, we consider the ensemble average over
all bubbles in our system and over similar lag times to get a statistically
consistent MSD. Meanwhile, the strain MSD can be estimated using the
positional MSD or mean-squared displacement introduced earlier (see Fig. 2a).
These quantities can further be related using the modified Fourier transformed
(FT) version of the above equation [19, 20, 13]:

${\lvert{G^{*}(\omega)}\rvert}^{2}\simeq\frac{\widetilde{\Delta\sigma^{2}}(\omega)}{3\pi\left<a\right>\widetilde{\Delta\mathbf{r}^{2}}(\omega)}$
(6)

To avoid assumptions and approximations related to computing Fourier
transforms of these MSDs over a finite range of lag times, [20, 13], we
consider the exact convolutional relation described in Section B.1. This
equation can be further modified using the Wiener-Khinchin theorem and the
relationship between autocorrelation and MSD for the stress and strain, giving
us the following equation:

$\displaystyle\begin{split}2J^{2}(0)\langle\sigma^{2}\rangle+\int_{0}^{\tau_{i}}{f(\tau_{i}-t^{\prime})(\langle\sigma^{2}\rangle-\langle\Delta\sigma^{2}\rangle(t^{\prime})/2)dt^{\prime}}\\\
=(\langle\gamma^{2}\rangle-\langle\Delta\gamma^{2}\rangle(\tau_{i})/2)\\\
\simeq 3\pi\langle
a\rangle(\langle\mathbf{r}^{2}\rangle-\langle\Delta\mathbf{r}^{2}\rangle(\tau_{i})/2)/{\langle
a\rangle}^{3}\\\ \text{where $f(\tau_{i})$ is defined as follows,}\\\
f(\tau_{i})=\left(\int_{0}^{\tau_{i}}{\dot{J}(\tau_{i}-t")\dot{J}(t")dt"}+2J(0)\dot{J}(\tau_{i})\right)\\\
\end{split}$ (7)

where ${\Delta\sigma^{2}(\tau)}$ and ${\Delta\mathbf{\gamma}^{2}(\tau)}$
represent the time-averaged, mean-squared difference of the bubbles’ stress
and strains, in our analyses. We approximate the strain using the position
vector, $\mathbf{r}$ [24, 19] as discussed above. Further, we ensemble average
our MSDs over $4$ simulation runs. Finally, to represent the creep compliance,
we use a modified version of the model suggested by Lavergne and co-authors in
Ref. [8]: $J(t)=1/G_{\infty}+k_{D}/G_{\infty}[(1+{t/\tau_{0}})^{\beta}-1]$
(more details in Section B.2). Using this as a model for the viscoelastic
rheology for the foam, we undertake a simultaneous fitting operation for the
parameters of the model, i.e., $G_{\infty},k_{D}$ and $\beta$, at various lag
times or $\tau_{i}$ in the convolutional integral equation shown above (Eq.
7). Using the optimal parameters from the fit gives us the creep compliance
and, subsequently, the complex modulus $G^{*}(\omega)$ using the relation,
$G^{*}(\omega)J(\omega)=1/{i\mkern 1.0mu}\omega$. Further details of the
derivation and mathematics of the numerical procedure are provided in Section
B.1. It may be further noted that attempts to model the rheology using a
Maxwell model produced inferior solutions to Eq. 7, with the power-law model
cited above providing significantly better fits.

(a)(a)

(a)

(a)

(a)

(b)

Figure 3: We compute the viscoelastic moduli for the dynamic viscous
simulations considered from the fluctuating stresses and displacements of
bubbles in the simulation, as described in the text. $\xi$ values with data
over a significant $\tau$ range were considered for the calculation. The
dotted grey lines indicate the $\tau$ range of the MSD data used for the above
calculation. (a) Fitting the model explained in the text to simulation data
gives us suitable fits with a family of curves with power-law behavior. The
creep compliance scales as $J(\tau)\sim t^{\beta}$ in the lag time range shown
above. (b) $G^{*}(\omega)$ obtained from $J(\tau)$, gives us power-law
rheology in $\omega$ i.e., $G^{*}(\omega)\sim\omega^{\beta}$. This behavior is
observed at all $\xi$ values calculated above. (inset) The predicted $\beta$
values, indicative of the log-slope for the curves in (b), hover consistently
in the range $\sim 0.15-0.2$, similar to previously observed values in
simulation [13], and experiments [8, 9].

The results from the computed creep compliance and dynamic shear moduli are
summarized in Fig. 3. $G^{*}(\omega)$ exhibits a power-law regime over the
$\omega$ range of interest and is characteristic of behavior predicted in
theory [1], simulations [13] and observed in experiments [10, 8, 9]. Recent
experiments and our simulation results here (see Fig. 3), evaluated with a
robust numerical approach provide clear evidence in support of the existence
of power-law rheology in SGMs. Fig. 3a, shows the fits for the $J(t)$ model
described above, with a family of curves with similar power-low exponents.
Considering the semi-analytical FT to obtain $G^{*}(\omega)$ gives us the
viscoelastic moduli with a power regime defined by
$G^{*}(\omega)\sim\omega^{\beta}$, with weak dependence of the exponent
$\beta$ on damping, showing that this is a universal feature for foams,
regardless of damping.

###  II.3 Memory and recovery in perturbed SGMs

(a)(a)

(a)

(a)

(a)

(b)

Figure 4: Scrambling a quasi-static system shows an immediate return to trend.
(a) We scramble the configurational positions of a system in steady state (at
$t=400$) in a quasi-static ($\xi=0$) simulation. Surprisingly, the system
always finds a ’new’ steady state right away, as indicated by the coordination
number ($z-z_{C}$) measured here, and continues to evolve with similar dynamic
properties. The dark and light symbols represent the scrambled and unscrambled
simulation, respectively. (b) Running multiple $(\sim 100)$ such scrambles at
$t=400$, gives us a Gaussian distribution of $\langle z\rangle$ as shown
above. This overlaps well with $\langle z\rangle$ values obtained at $t=400$
for $10$ different realizations of the same simulation as indicated by the
mean $\pm$ standard deviations. This tells us that the scrambled simulation
returns to the newly found steady state instantaneously. (inset) Moreover, the
temporal autocorrelations for these $z$ ensembles - scrambled and unscrambled
- provide similar decorrelation times. These findings indicate similar dynamic
properties for the scrambled and unscrambled simulation.

The SGM system shows a significant downhill descent in energy as the largest
bubbles coarsen and grow. As this downward trend continues, the system reaches
a dynamical scaling steady state [13]. While it is unclear whether
configurations in this regime form an ergodic ensemble over some
characteristic time, the bubbles show stable trends in various structural
quantities like average coordination number, mean bubble radius, normalized
radial distribution, etc. Bubbles initially move around to reach the steady
state, defined by the dynamical scaling ‘attractor’ on the energy landscape,
and then continue to evolve in this steady state ensemble. Any perturbation
away from the attractor would thus lead the system back to a ’new’ steady
state as defined by the structural and dynamical properties of the attractor
and the system landscape. Experiments have observed [14, 15] that a strain-
perturbed foam relaxed back to its unperturbed steady state after an
unexpectedly long waiting time, and have described this as a memory phenomenon
or measure of history dependence. The consensus [14, 15] on the origin of this
memory is that coarsening mediated excitations are needed to enable the system
to overcome local minima that the perturbed system relaxes into. Thus, the
long waiting time has been considered a result of slow coarsening.

To study this phenomenon’s structural and dynamical significance
computationally, we run a set of simulations using our modified damped model
over various $\xi$ values. We consider the theoretical extreme of a
perturbation by introducing positional scrambles in our system. To do so, we
begin with a typical steady-state system and randomly scramble the various
$3N$ positions of the bubbles. This scramble randomly assigns a point in
hyperspace for the system of soft spheres, providing a random structural
perturbation. We then continue with the relaxation-coarsening procedure
described previously in Section II.1. It must be noted here that for the
quasi-static case when $\xi=0$, we relax the system to its first energy
minimum (i.e. mechanical equilibrium) using FIRE [29] instead of using Eq. 2.

Figure 5: Scrambling a quasi-static system shows no change to ripening
evolution. Here, we look at the structure through the radial distribution
formed at steady state for a scrambled and unscrambled system. As previously
in Fig. 4, the system instantaneously continues in steady state. As can be
noticed, the slope changes for the scrambled simulation at $t=400$ (indicated
by arrows), indicative of a new foam initiation time [13]; however, the trend
remains linear, consistent with dynamic scaling state behavior ${\langle
a\rangle}^{2}\sim t_{\text{age}}$.

For the quasi-static case, we see that the system, upon one (or even multiple)
scrambles, returns to the earlier dynamical scaling steady trend (see Fig. 4a)
immediately. Indicators like $\langle z\rangle$ and $\langle a\rangle^{2}$
show no significant change from steady-state behavior, as can be seen in Fig.
4 which plots the scrambled (at $t=400$) and unscrambled average coordination
number as a function of time. Here the scrambled system experiences no
barriers to reaching this ’new’ steady state with FIRE traversing the large
configurational distance on a relatively smooth portion of the energy
landscape (at higher energies) to find the nearest (primary) minima. It may be
noted that the scramble moves the system to a random $3N$-dimensional
configuration on the energy landscape. The system thus evolves in the
particular meta-basin corresponding to the scrambled positions going forward.
So while a single simulation might not erogodically explore all portions in
configuration space, these different hyperspaces on the energy landscape have
similar structural properties. Small changes in moving ensemble averages
indicate the slight variations in different regions or metabasins of the
energy landscape. Effectively, all these primary minima that the minimizer
finds belong to the ’steady-state ensemble’ of the particular foam radii
distribution realization at $t=400$.

To test whether the scrambled simulation actually returns to the ’steady-state
ensemble’ for a similar foam at the same age, we compare the average
coordination numbers at $t=400$ for $10$ different quasi-static simulations
(different positional initializations at $t=0$) in Fig. 4b; with a pool of
coordination numbers obtained from scrambling the same test simulation (from
Fig. 4a) configuration at $t=400$, a $100$ different ways and quenching them
using FIRE. As seen in Fig. 4b, the distribution of average coordination
numbers for a quenched minima from the test simulation lies in the range of
expected steady-state $\langle z\rangle$ values for a similar foam simulation
at the same age. Lastly, looking at the inherent temporal correlations in the
coordination number further tells us that similar correlations get rebuilt
after the system evolves on the randomly chosen ’metabasin’ on the energy
landscape Fig. 4a(inset). Thus we conclude that while the system doesn’t
return to the exact configurational hyperspace on the energy landscape, the
foam exhibits a ’memory’ effect in various physical and dynamical properties.

Interestingly, this also doesn’t affect the coarsening mediated bubble size
distribution reached at dynamic scaling as seen in Fig. 5. This can be seen in
the average system radii measured over time in Fig. 5. While there is a
noticeable change in the rate of radial change or the slope, the trajectory
remains in steady state as indicated by the linear gradient. Additionally, the
scrambled simulation continues to evolve with similar moments of the radii
distribution (see Fig. 4a(inset)). Overall, this shows that while any steady
state structure built in by the dynamics and coarsening before the scramble
gets ruined by the perturbation, it gets restored immediately by quenching to
the nearest minima (using FIRE).

(a)(a)

(a)

(a)

(a)

(b)

Figure 6: Scrambled foams with damping show very slow relaxation towards their
prior trends in energy and coordination. Finite $\xi$ simulations are
scrambled at $t=400$ and evolved using the dynamical equation Eq. 2. (a) The
potential energies post scrambling show a progressive trend in reaching the
steady state. Fitting the energies (b) The mean coordination number $\langle
z\rangle$ of the system shows similar $\xi$ dependent trends. Interestingly,
the initial dynamics directs the system to lower $\langle z\rangle$
configurations before relaxing to the appropriate steady-state value.

Repeating the same computational experiment at finite $\xi$ provides insight
into the mechanism of the memory phenomenon. In agreement with previous
experiments [10, 15], we see that the system requires a surprisingly long time
to recover to its former steady state trend (see Fig. 6). However, unlike
previous suggestions of this time-scale being coarsening-mediated, we observe
a $\xi$ dependent phenomenon. This viscous time scale dictates the time the
system takes to relax any energetic stress built in by the overlaps caused by
the positional scramble. Larger $\xi$ leads to a longer time for the system to
relax these unstable overlaps. Fig. 6 shows the progression of energy and
average coordination number towards equilibrium after a positional scramble.
Interestingly, the $z$ values shoot below the steady state line post scramble
before trending back to steady state, much like previous experiments measuring
rearrangement rate [14]. Finally, one may note here that while the waiting
times seem to be damping dependent, they are, however, much larger than
$\tau_{R}=\xi\langle a\rangle/\epsilon\sim\mathcal{O}(\xi)$. This discrepancy
might be due to the extreme nature of positional perturbation introduced in
our simulation – which introduces many large and small perturbations for all
$N$ bubbles away from the nearest steady-state ensemble configuration. Thus,
the waiting time is a compounded sum of all these different distances that the
bubbles must traverse to reach the ’new’ steady state.

##  III Conclusions

We have shown that the ‘bubble model’ with simple damping and simulated
ripening recreates many of the exciting phenomena reported for soft-glassy
materials. Specifically, this model exhibits avalanchey, intermittent dynamics
at low viscosities with non-Brownian super-diffusive motion. Considered in
high-dimensional configuration space, such motion occurs along a fractal
configurational path that is constrained to the lowest energy portions of the
potential energy landscape. We find that the energy minima of the bubble model
are clustered together in configuration space, scattered along the
configuration path followed by the model (at low viscosity). In the practical
absence of viscous stresses, the system hops from each (ripening destabilized)
energy minimum to a nearby, adjacent energy minimum. This landscape-dominated
motion subsequently produces the observed super-diffusive motion and stress
and strain fluctuations corresponding to power-law rheology. As the simulated
viscosity is increased, the system shows progressively smoother dynamics and
motion with a more ballistic character. In this case, viscous stresses cause
the configuration never to explore the true potential energy minima but
instead evolve along a path that is adjacent to the cluster of energy minima.
The system effectively stays at higher potential energy, and the finer details
of the fractal configuration path seen at lower viscosity are washed out,
leading to a straighter path and more ballistic motion. Between these two
limits, one can find a gradation of properties, where the system displays
increasing characteristic length and time scales above which the low viscosity
behavior may still be observed.

This model also successfully generates power-law rheology, which previously
has not been reported for a damped bubble model. However, unlike other
properties, power-law rheology seems to be a consistent feature in SGMs over a
wide range of viscosity values. This suggests an extended fractal nature for
the energy landscape, consistently producing power-law rheology even when the
configurations are at energy somewhat above the energy minima. Further, our
microrheology-based approach provides a robust and reliable way to compute
viscoelastic moduli from force and strain fluctuations measurements of
constituent particles and is free of any systematic truncation errors
associated with earlier microrheological methods.

Lastly, we investigate the ‘memory’ of the ripening bubble model for
mechanical perturbations by randomly scrambling the bubble positions. We find
that scrambled configurations (effectively a random point in configuration
space) must relax a long configuration distance before reaching their first
potential energy minimum. Moreover, those first energy minima are
indistinguishable (statistically) from the ensemble of configurations explored
by other ripening simulations of the same age. This is most clearly shown by
the quasi-static simulation that immediately recovers its earlier (ensemble
averaged) baseline properties when it reaches its first energy minimum. For
finite viscosity, the system can take a long time to traverse the required
configuration distance to return to the vicinity of the energy minima cluster
and recover its earlier mechanical properties. In viscous systems, this
recovery time (and effective ‘memory’ time) is proportional to the viscous
relaxation time of the model. This viscosity-mediated recovery process is
contrary to previous experimental inferences and a consequence of the barrier-
free potential energy landscape at higher energies that the perturbed system
must traverse.

Future work would include developing a model that captures other long-time
characteristics of these SGMs. We hope that such a description will provide a
more complete and practical model for SGMs. Modeling the physical properties
of the many materials categorized under SGMs could be of potential use in
fields ranging from material science (foams and complex fluids) to biology
(living cells [4, 5]).

## Author Contributions

A.T., R.A.R., and J.C.C. designed research; A.T. performed research and
analyzed data; A.T., R.A.R., and J.C.C. wrote the paper. R.A.R. and J.C.C.
contributed equally to this work.

###### Acknowledgements.

We are grateful for valuable conversations with Douglas Durian, François
Lavergne, Andrea Liu, Talid Sinno, and Véronique Trappe. This work was
supported by NSF-DMR 1609525 and 1720530 and computational resources provided
by XSEDE through TG-DMR150034.

##  Appendix A Damped SGM Model

###  A.1 Underdamped limit of an overdamped equation

The equation used for the simulations, as described in the main text, is:

$\displaystyle\begin{split}\xi\frac{d\mathbf{r}_{i}}{dt}&=\mathbf{F}_{i}\\\
&=-{{\sum}}_{j}^{nn}\frac{\partial V(r_{ij})}{\partial r_{i}}\end{split}$ (8)

It can be seen that the equation is similar to an overdamped equation of
motion. However, it must be noted that though the equation resembles and has
the characteristics of an overdamped equation of motion, the same is not due
to a large viscosity but rather the non-inertial nature of the constituent
particles considered in the system. One may recall that dynamics for a mass
attached to a damped spring are mediated by the damping factor
$\zeta=b/(2\sqrt{km})$. That can be evaluated for our system of interest as
follows $\zeta\simeq\xi/\sqrt{\epsilon\rho\langle a\rangle}$. Since overdamped
dynamics is achieved when $\zeta\geq 1$, we see that the non-inertial
particles $(\rho\to 0)$ in our case give rise to the so-called overdamped
equation of motion. Meanwhile, we continue to operate with a finite value of
$\xi$.

###  A.2 Integration and stability

The simulation can be summarized as a numerical integration of the two
equations – Eq. 2 and eq. 3 using a numerical integration technique. Due to
the stiff nature of Eq. 2 (especially at small $\xi$ values, one needs to
choose appropriate $dt$ values to ensure any error perturbations don’t diverge
as the simulation proceeds and that the solution is a converged one. We use a
simple Explicit Euler scheme to perform our integration here. We note that
other methods, like implicit Euler and second-order Runge Kutta scheme,
provide more extensive stability regimes for $dt$ and are more accurate but
can have more significant computational overload associated with the
integration scheme. Below, we perform a simple numerical stability test.

We start by considering Eq. 2 for all $N$ particles or $3N$ degrees of
freedom, i.e., $i\in\\{1,2,...3N\\}$, which can be expressed in terms of the
Hessian for using a Taylor expansion as follows: and

$\displaystyle\begin{split}\xi\frac{d\mathbf{r}}{dt}&=\mathbf{F}\\\
&=\mathbf{F}_{0}-\mathbf{H}\mathbf{r}\end{split}$ (9)

where $\mathbf{r}$, $\mathbf{F}$ are $3N$ dimensional vectors and $\mathbf{H}$
is a $3N\times 3N$ matrix or the Hessian of the potential field. We may note
here that for most $\xi$ simulations, the system configurations are close to
mechanical equilibrium, so for our stability analysis, we may approximate this
using $\mathbf{F}_{0}\simeq 0$. Further one may note that the any error
$\epsilon_{i}$ would propagate via an equation similar to Eq. 9:

$\displaystyle\begin{split}\xi\frac{d\mathbf{\epsilon}}{dt}&=-\mathbf{H}\mathbf{\epsilon}\end{split}$
(10)

Now, using the Explicit Euler formalism, for time steps $n+1$ and $n$, we get:

$\displaystyle\begin{split}\xi\frac{\mathbf{\epsilon_{n+1}}-\mathbf{\epsilon_{n}}}{dt}&=-\mathbf{H}\mathbf{\epsilon}_{n}\\\
&=-\mathbf{\lambda}\mathbf{\epsilon}_{n}\\\
\frac{||\mathbf{\epsilon_{n+1}}||}{||\mathbf{\epsilon_{n}}||}&=||\mathbf{I}-\mathbf{\lambda}||dt/\xi\end{split}$
(11)

where $\mathbf{\lambda}$ is a matrix containing all eigenvalues of
$\mathbf{H}$. Enforcing the criteria of stability on the equation above we
have,

$\displaystyle\begin{split}\frac{||\mathbf{\epsilon_{n+1}}||}{||\mathbf{\epsilon_{n}}||}&\leq
1\\\ ||\mathbf{I}-\mathbf{\lambda}||dt/\xi&\leq 1\\\
0\leq||\lambda_{max}||dt/\xi&\leq 2\end{split}$ (12)

where $\lambda_{max}$ is the largest eigenvalue for $\mathbf{H}$. Since all
eigenvalues would be real for this physical system, we now have,

$\displaystyle\begin{split}0\leq||\lambda_{max}||dt/\xi&\leq 2\\\ dt&\leq
2\xi/\lambda_{max}\end{split}$ (13)

For most configurations, explored in our system simulation the $\lambda_{max}$
varies around $\sim 1-10$. This gives us that $dt\leqslant\xi/5$ is the
condition for stability. Here, we choose $dt=\xi/10$, as the step size for all
simulations reported in this study. Since we have Eq. 3, which also controls
overall dynamics, this choice of time-step was validated for convergence. We
have verified that our explicit Euler scheme was converged by checking other
smaller values of $dt$. Other schemes like the RK-2 also produced similar
results. Further, for $\xi>0.01$, we stuck with the use of $dt=0.001$, as the
system moves further away from mechanically stable states and the above
approximation in Eq. 10 fails to strictly hold.

###  A.3 Dimensionless group analysis: Deborah number

Apart from evaluating the Deborah number $De$ as the ratio of the damped
relaxation time from Eq. 2 ($\tau_{R}=\xi{\langle a\rangle}^{2}/\epsilon$) and
probing time associated with changing bubble radii $(a)$ imparted by the
coarsening process, Eq. 3 ($\tau_{C}={\langle a\rangle}^{2}/\alpha_{1}$, we
can do a simple Buckingham Pi analysis to determine the relevant $\Pi$ group.
Below, we present the analysis to derive the Deborah number as a $\Pi$ group.

One can re-model the system through an experimental lens and pose the problem
statement as measuring the average radii $\langle a\rangle$ as a function of
time. Intuitively, this might be influenced by system properties like
$\epsilon$, $\rho$, $\alpha_{1}$, $\xi$. These $4$ quantities along with
$\langle a\rangle$, are comprised of the dimensions $M$, $L$ and $T$. Thus $2\
\Pi$ groups can be made using these variables for every combination of $3$
repeating variables being chosen. Here, we choose $\rho$, $\alpha$, and $\xi$
as are repeating variables.

$\displaystyle\begin{split}\Pi_{1}&=f(\epsilon,\rho,\alpha_{1},\xi)\\\
&=\epsilon\rho^{x}\alpha_{1}^{y}\xi^{z}\end{split}$ (14)

Solving for $x$, $y$, and $z$ so that $\Pi_{1}$ is dimensionless, we get
$\Pi_{1}=\epsilon/(\alpha_{1}\xi)$ or $De=\xi\alpha_{1}/\epsilon$.

##  Appendix B Rheology

###  B.1 Analytical Derivation

Here, we provide a derivation for the integral equation Eq. 7, which we used
to compute the viscoelastic moduli for our simulation. We start by noting that
the theory of viscoelasticity for linear materials [25, 26] shows that the
creep compliance $J$, can be related to the stress $\sigma$ and strain
$\gamma$ as follows:

$\displaystyle\begin{split}\int_{-\infty}^{t}{J(t-t^{\prime})\dot{\sigma}(t^{\prime})dt^{\prime}=\gamma(t)}\end{split}$
(15)

Here we may note that $J$, $\sigma$ and $\gamma$ are $=0\ \forall
t\in(-\infty,0)$ and $\geq 0\ \forall t\in[0,\infty)$. Thus, we can extend the
integral limits by doing the following:

$\displaystyle\begin{split}\int_{-\infty}^{\infty}{J(t-t^{\prime})\dot{\sigma}(t-t^{\prime})dt^{\prime}}&=\gamma(t)\\\
\int_{-\infty}^{\infty}{\dot{J}(t-t^{\prime}){\sigma}(t-t^{\prime})dt^{\prime}}&=\gamma(t)\\\
\text{ using the product rule}\\\
\int_{0}^{t}{\dot{J}(t-t^{\prime}){\sigma}(t-t^{\prime})dt^{\prime}}+J(0)\sigma(t)&=\gamma(t)\end{split}$
(16)

Taking the Fourier Transform of the non-decomposed equation above and applying
the convolution theorem gives us,

$\displaystyle\begin{split}{\widetilde{\dot{J}}}\ \widetilde{{{\sigma}}}&=\
\widetilde{\gamma}\\\ \end{split}$ (17)

While one could potentially work with Eq. 16 or Eq. 17, the numerical
inaccuracies associated with an FT [13, 20] and the statistical noise in a
trajectory function like $\sigma(t)$ or $\gamma(t)$, would make the procedure
more difficult. Thus, we use the Wiener–Khinchin theorem and further transform
the auto-correlation into its mean squared version as follows:

$\displaystyle\begin{split}{\widetilde{\dot{J}}}\
{\widetilde{\dot{J}}}&=\frac{{||\widetilde{\gamma}||}^{2}}{{||\widetilde{\sigma}||}^{2}}\\\
&=\frac{\widetilde{R_{\gamma\gamma}}}{\widetilde{R_{\sigma\sigma}}}\\\
&=\frac{\widetilde{\langle\gamma^{2}\rangle-\langle\Delta\gamma^{2}\rangle/2}}{\widetilde{\langle\sigma^{2}\rangle-\langle\Delta\sigma^{2}\rangle/2}}\end{split}$
(18)

Reshuffling this equation and taking the inverse FT yields an integral
equation. We decompose the limits to stay between $0$ and $t$, which adds a
few boundary terms for the step function jump in $J$ and $\sigma$ at $t=0$.
Further, we change our notation for $t$ to $\tau$, to be consistent with the
MSDs which are calculated as averages over lag times.

$\displaystyle\begin{split}2J^{2}(0)\langle\sigma^{2}\rangle+\int_{0}^{\tau}{f(\tau-t^{\prime})(\langle\sigma^{2}\rangle-\langle\Delta\sigma^{2}\rangle(t^{\prime})/2)dt^{\prime}}\\\
=(\langle\gamma^{2}\rangle-\langle\Delta\gamma^{2}\rangle(\tau)/2)\\\
\text{where $f(\tau)$ is defined as follows,}\\\
f(\tau)=\left(\int_{0}^{\tau}{\dot{J}(\tau-t")\dot{J}(t")dt"}+2J(0)\dot{J}(\tau)\right)\\\
\end{split}$ (19)

This equation can be approximated, using similar mathematical approximations
as used earlier in Ref. [20].

$\displaystyle\begin{split}\int_{0}^{\tau}{g(\tau-t^{\prime})\langle\Delta\sigma^{2}\rangle(t^{\prime})dt^{\prime}}=\langle\Delta\gamma^{2}\rangle(\tau)\\\
\text{where $g(\tau)$ is defined as follows,}\\\
g(\tau)=\left(\int_{0}^{\tau}{\dot{J}(\tau-t")\dot{J}(t")dt"}\right)\\\
\end{split}$ (20)

We approximate the right-hand side of this equation using the bubble portions
$\mathbf{r}$ [24, 19], giving: $\simeq 3\pi\langle
a\rangle(\langle\mathbf{r}^{2}\rangle-\langle\Delta\mathbf{r}^{2}\rangle(\tau_{i})/2)/{\langle
a\rangle}^{3}$. However, it may be noted that this equation is not well
defined at $\tau=0$. Thus we evaluate this only for lag times greater than
zero. To get an accurate solution, we consider the above equation at various
finite lag time values or $\tau_{i}$ and solve a set of simultaneous equations
to find the appropriate creep compliance, $J(t)$. Specifically, we choose
$\tau_{i}\in\\{\tau_{1},\tau_{2},\tau_{3}...\tau_{max}\\}$. Here, $\tau_{1}$
can be as small as $dt$. We report here results for $\tau_{1}=1$. This choice,
however, brings in some numerical error due to the integrals going from
$0\to\tau$. It may be noted that this equation is mathematically exact for
$\forall\tau>0$ and that the upper limit of our observation – $\tau_{max}$,
does not affect the numerical procedure, effectively avoiding a source of
truncation error present in many earlier approaches.

###  B.2 Choice of Fitting Model

One may notice that solving Eq. 19 or Eq. 20 requires a model for $J(t)$. Here
we choose a modified version of the model suggested in Ref. [8]. The original
model put forth in the above study has a terminal mode of relaxation at long
times, given by $t/\eta_{R}$, and has been observed previously in experiments
[10]. In our simulations, we, however, do not observe any terminal relaxation
and thus ignore the additional term mentioned above. We considered a modified
version of the model given as follows.

$\displaystyle\begin{split}J(t)=1/G_{\infty}+k_{D}/G_{\infty}[(1+{t/\tau_{0}})^{\beta}-1]\end{split}$
(21)

## References

  * Sollich  _et al._ [1997] P. Sollich, F. Lequeux, P. Hébraud, and M. E. Cates, Rheology of soft glassy materials, Physical Review Letters 78, 2020–2023 (1997).
  * Sollich [1998] P. Sollich, Rheological constitutive equation for a model of soft glassy materials, Physical Review E 58, 738–759 (1998).
  * Hébraud and Lequeux [1998] P. Hébraud and F. Lequeux, Mode-coupling theory for the pasty rheology of soft glassy materials, Physical Review Letters 81, 2934–2937 (1998).
  * Bursac  _et al._ [2005] P. Bursac, G. Lenormand, B. Fabry, M. Oliver, D. A. Weitz, V. Viasnoff, J. P. Butler, and J. J. Fredberg, Cytoskeletal remodelling and slow dynamics in the living cell, Nature Materials 4, 557–561 (2005).
  * Hoffman  _et al._ [2006] B. D. Hoffman, G. Massiera, K. M. V. Citters, and J. C. Crocker, The consensus mechanics of cultured mammalian cells, Proceedings of the National Academy of Sciences 103, 10259–10264 (2006).
  * Hoffman and Crocker [2009] B. D. Hoffman and J. C. Crocker, Cell mechanics: Dissecting the physical responses of cells to force, Annual Review of Biomedical Engineering 11, 259–288 (2009).
  * Ikeda  _et al._ [2012] A. Ikeda, L. Berthier, and P. Sollich, Unified study of glass and jamming rheology in soft particle systems, Phys. Rev. Lett. 109, 018301 (2012).
  * Lavergne  _et al._ [2022] F. A. Lavergne, P. Sollich, and V. Trappe, Delayed elastic contributions to the viscoelastic response of foams, The Journal of Chemical Physics 156, 154901 (2022), https://doi.org/10.1063/5.0085773 .
  * Rodriguez-Cruz  _et al._ [2022] C. Rodriguez-Cruz, M. Molaei, A. Thirumalaiswamy, K. Feitosa, V. N. Manoharan, S. Sivarajan, D. H. Reich, R. A. Riggleman, and J. C. Crocker, Fractal landscape dynamics in dense emulsions and stock prices (2022).
  * Gopal and Durian [2003] A. D. Gopal and D. J. Durian, Relaxing in foam, Phys. Rev. Lett. 91, 188303 (2003).
  * Durian [1995] D. J. Durian, Foam mechanics at the bubble scale, Physical Review Letters 75, 4780–4783 (1995).
  * Durian [1997] D. J. Durian, Bubble-scale model of foam mechanics: Melting, nonlinear behavior, and avalanches, Physical Review E 55, 1739–1751 (1997).
  * Hwang  _et al._ [2016] H. J. Hwang, R. A. Riggleman, and J. C. Crocker, Understanding soft glassy materials using an energy landscape approach, Nature Materials 15, 1031–1036 (2016).
  * Gopal and Durian [1995] A. Gopal and D. J. Durian, Nonlinear bubble dynamics in a slowly driven foam, Physical Review Letters 75, 2610 (1995).
  * Höhler  _et al._ [1999] R. Höhler, S. Cohen-Addad, and A. Asnacios, Rheological memory effect in aqueous foam, Europhysics Letters (EPL) 48, 93 (1999).
  * Ono  _et al._ [2003] I. K. Ono, S. Tewari, S. A. Langer, and A. J. Liu, Velocity fluctuations in a steadily sheared model foam, Phys. Rev. E 67, 061503 (2003).
  * Stevenson [2010] P. Stevenson, Inter-bubble gas diffusion in liquid foam, Current Opinion in Colloid & Interface Science 15, 374 (2010).
  * Tewari  _et al._ [1999] S. Tewari, D. Schiemann, D. J. Durian, C. M. Knobler, S. A. Langer, and A. J. Liu, Statistics of shear-induced rearrangements in a two-dimensional model foam, Phys. Rev. E 60, 4385 (1999).
  * Lau  _et al._ [2003] A. W. C. Lau, B. D. Hoffman, A. Davies, J. C. Crocker, and T. C. Lubensky, Microrheology, stress fluctuations, and active behavior of living cells, Phys. Rev. Lett. 91, 198101 (2003).
  * Crocker and Hoffman [2007] J. C. Crocker and B. D. Hoffman, Multiple‐particle tracking and two‐point microrheology in cells, Methods in Cell Biology Cell Mechanics , 141–178 (2007).
  * O’Hern  _et al._ [2002] C. S. O’Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, Random packings of frictionless particles, Phys. Rev. Lett. 88, 075507 (2002).
  * Feitosa  _et al._ [2006] K. Feitosa, O. L. Halt, R. D. Kamien, and D. J. Durian, Bubble kinetics in a steady-state column of aqueous foam, Europhysics Letters 76, 683 (2006).
  * Liu and Nagel [1998] A. J. Liu and S. R. Nagel, Jamming is not just cool any more, Nature 396, 21–22 (1998).
  * Mason [2000] T. G. Mason, Estimating the viscoelastic moduli of complex fluids using the generalized stokes-einstein equation, Rheologica Acta 39, 371 (2000).
  * Volterra [1909] V. Volterra, Sulle equazioni integro-differenziali della theoria dell’ elasticita, Atti Reale Accad. naz. Lincei. Rend. Cl. sci. fis., mat. e natur. 18, 295 (1909).
  * Volterra and Joseph [1913] V. Volterra and P. Joseph,  _Leçons sur les fonctions de lignes: Professés a la sorbone en 1912_ (Gauthier-Villars, 1913).
  * Rossi and Testa [2010] G. C. Rossi and M. Testa, The stress tensor in thermodynamics and statistical mechanics, The Journal of Chemical Physics 132, 074902 (2010).
  * Zhou [2003] M. Zhou, A new look at the atomic level virial stress: on continuum-molecular system equivalence, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 459, 2347 (2003).
  * Bitzek  _et al._ [2006] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Structural relaxation made simple, Phys. Rev. Lett. 97, 170201 (2006).