# Thixotropic spectra and Ashby-style charts for thixotropy

Samya Sen  Randy H. Ewoldt<EMAIL_ADDRESS>Department of Mechanical
Science and Engineering,
University of Illinois at Urbana-Champaign,
Urbana, IL 61801, USA

###### Abstract

There is no universal model for thixotropy, and comparing thixotropic effects
between different fluids is a subtle yet challenging problem. We introduce a
generalized (model-insensitive) framework for comparing thixotropic properties
based on thixotropic spectra. A superposition of exponential stress modes
distributed over thixotropic timescales is used to quantify buildup and
breakdown times and mode strengths in response to step-change input. This
mathematical framework is tested with several experimental step-shear rate
data on colloidal suspensions. Low-dimensional metrics based on moments of the
distribution reveal characteristic average thixotropic properties which are
visualized on Ashby-style diagrams. This method outlines a framework for
describing thixotropy across a diverse range of microstructures, supporting
scientific studies as well as material selection for engineering design
applications.

††preprint: AIP/JoR/preprint

##  I Introduction

Thixotropic phenomena are observed in fluids with a “structure” that evolves
with time, and such fluids are common both in industry and in daily life.
Common examples include flow batteries [1, 2, 3, 4], crude oils [5], food
materials [6], blood [7, 8], and colloidal suspensions [9, 10, 11, 12, 13],
among numerous others [14, 15, 16, 17]. As defined by Mewis and Wagner [16,
18], thixotropy is “the continuous decrease of viscosity with time when flow
is applied to a sample that has been previously at rest, and the subsequent
recovery of viscosity when the flow is discontinued”. Real materials show both
thixotropy and viscoelasticity, and microstructural buildup can cause an
increase in the viscosity as well as the elastic modulus [19]. In many cases
thixotropy may be distinguished from viscoelasticity based on timescales:
viscoelastic timescales are typically shorter compared to the thixotropic
timescales [20, 21].

Materials with short thixotropic timescales (e.g. $\lesssim\mathcal{O}(1)$ s)
are typically labeled as _not_ thixotropic, but the degree of thixotropy of a
material depends on not only the characteristic timescales, but also the
amount of change in a rheological property of interest, e.g. shear viscosity
or elastic modulus, or more generally a change in state of stress. Stress
changes, and timescales over which these changes occur, must both be
considered to quantify and compare the degree of thixotropy across different
materials. Additionally, thixotropic materials involve a range of
microstructural lengthscales and associated timescales, and this makes
analyzing the rheology non-trivial. Much like viscoelastic materials, these
are characterized by a large diversity of elementary units with multiple
characteristic timescales [22, 23]. There is no universal _predictive_ model
for thixotropy. Therefore, to enable comparison between materials, we
introduce a universal _descriptive_ model. Our objective is to develop
methods, irrespective of the material tested or the underlying predictive
constitutive model, to quantify thixotropy using observed stress changes
($\Delta\sigma$) and timescales ($\tau_{\rm char}$). This concept is
illustrated in Fig. 1.

Figure 1: Proposed Ashby-style mapping for comparing thixotropic responses.
Even with a chosen experimental protocol, the dependencies of two relevant
thixotropic properties on the experimental conditions makes the parameter
space high-dimensional.

The simple idea of using $\Delta\sigma$ and $\tau_{\rm char}$ actually
involves underlying high-dimensionality, noted in Fig. 1. Both thixotropic
breakdown and buildup processes must be considered, and the associated
kinetics warrant the use of a distribution of timescales to describe the data,
eventually reduced to obtain a secondary, average metric of $\tau_{\rm char}$.
The amount of thixotropic change depends on the experimental protocol used;
the four most common ones are shown in Fig. 2. Using shear stress to study the
changes (and hence $\Delta\sigma$) allows us to include information about
viscosity (rate-controlled tests) or modulus (oscillatory tests). Lastly, both
$\Delta\sigma$ and $\tau_{\rm char}$ are functions of the specific
experimental conditions native to the chosen protocol. In step shear rate
tests, $\Delta\sigma$ and $\tau_{\rm char}$ are functions of the initial
($\dot{\gamma}_{\rm i}$) and final ($\dot{\gamma}_{\rm f}$) shear rates, or
initial ($\gamma_{\rm i}$) and final ($\gamma_{\rm f}$) shear strain
amplitudes for oscillatory amplitude step shear (oscillatory) tests. Fig. 1 is
thus a two-dimensional projection of this high-dimensional space onto a plane.
On this plane we can still vary $\dot{\gamma}_{\rm i}$ and $\dot{\gamma}_{\rm
f}$. Other metrics in addition to $\Delta\sigma$ and $\tau_{\rm char}$ may
also exist, as noted in Fig. 1.

Here we introduce the concept of thixotropic spectra, revealed by step
changes, as a description of thixotropic effects. Step changes, rather than
hysteresis loops, allow thixotropic spectra to be revealed. Many other step
tests are available (Fig. 2), but we focus our attention on step shear rate
since this can distinguish thixotropy and viscoelasticity [18]. It may be a
step down ($\dot{\gamma}_{\rm f}<\dot{\gamma}_{\rm i}$) or step up
($\dot{\gamma}_{\rm f}>\dot{\gamma}_{\rm i}$) test, and the transient stress
evolution is used to study thixotropic structure buildup or breakdown. Fig.
2(a) and (b) show typical stress signals for step down (buildup) and step up
(breakdown) rate tests respectively. In step down tests, the stress may first
show a viscoelastic relaxation (decrease) at very short times, which may or
may not be observable [14, 16]. This may be followed by thixotropic recovery
at long times [14, 16]. In step up tests, the stress typically shows a
viscoelastic stress increase at short times, followed by thixotropic decay at
long times [14, 16]. For a given test, the data is analyzed to obtain
necessary information about $\Delta\sigma(\dot{\gamma}_{\rm
i},\dot{\gamma}_{\rm f})$ and $\tau_{\rm char}(\dot{\gamma}_{\rm
i},\dot{\gamma}_{\rm f})$ while keeping in mind the dimensionality of the
problem. The analysis method is dictated by the protocol, both of which are
described in detail in the following section.

Figure 2: Example input protocols where a thixotropic spectrum of timescales
can be applied. (a1) Step down in steady shear ($\dot{\gamma}_{\rm
f}<\dot{\gamma}_{\rm i}$), (b1) step up in steady shear ($\dot{\gamma}_{\rm
f}>\dot{\gamma}_{\rm i}$), (c1) step down in oscillatory shear amplitude
($\gamma_{\rm f}<\gamma_{\rm i}$), (d1) step up in oscillatory shear amplitude
($\gamma_{\rm f}>\gamma_{\rm i}$). Typical response functions (stresses and
moduli) are shown in (a2)-(d2).

It must be noted that our approach only describes the thixotropic component of
the stress response; it requires the viscoelastic component to be sub-dominant
or negligible. As needed, the early-time viscoelastic response (as shown
schematically in Fig. 2(a2), (b2)) can be fit to a viscoelastic spectrum
superimposed on the thixotropic one. This will truncate the range of
observable thixotropic timescales. In our work, we only fit the thixotropic
part following any initial viscoelastic response.

##  II Theory

###  II.1 Flow conditions and rheometry protocols

We propose describing data in step tests as a superposition of stress
components evolving over a spectrum of thixotropic buildup or breakdown
timescales as a universal way to quantify thixotropic dynamics. The transient
thixotropic stress signal $\sigma$ evolving over time $t$, expressed as a
generalized discrete spectrum, is

$\displaystyle\sigma^{+}(t;\mathbb{P})$
$\displaystyle=\sigma_{0}+\sum\limits_{i=1}^{N}\sigma^{+}_{i}\left(1-e^{-t/\tau^{+}_{i}}\right),$
(1a) $\displaystyle\sigma^{-}(t;\mathbb{P})$ $\displaystyle=\sigma_{\rm
ss}+\sum\limits_{i=1}^{N}\sigma^{-}_{i}e^{-t/\tau^{-}_{i}},$ (1b)

where “$+$” is used for recovery (step down) and “$-$” for breakdown (step
up). Subscripts “0” and “ss” refer to initial $(t=0)$ and steady state
$(t\rightarrow\infty)$ values respectively. We assume the basis functions to
be _exponential_ for each individual mode. It is shown here for $N$ discrete
stress modes with mode strengths $\sigma_{i}$ and associated timescales
$\tau_{i}$. These modes can be either for recovery or breakdown, and
superscripts used to distinguish the two as $\tau^{+}_{i}$ or $\tau^{-}_{i}$,
but one would omit the “$\pm$” superscripts for generalized discussions.
$\mathbb{P}$ is the set of experimental parameters for a given protocol (e.g.
$\mathbb{P}=\\{\dot{\gamma}_{\rm i}$,$\dot{\gamma}_{\rm f}\\}$ in step rate
tests, a steady state at $\\{\dot{\gamma}_{\rm i}$ was achieved). The
distribution of $\sigma_{i}$ over $\tau_{i}$ can be reduced further to
generate low-dimensional summarizing metrics, as discussed in § II.4. The key
idea is to represent a signal as the sum (or integral) of multiple component
signals, each with its own characteristic timescale of evolution. It can be
applied to step changes of unidirectional or oscillatory shear forcing; the
input scheduling for each shown in Fig. 2.

Step rate tests are most commonly used for studying thixotropy [16]. The exact
form of Eq. 1 applied to step rate tests is

$\displaystyle\sigma^{+}(t;\dot{\gamma}_{\rm i},\dot{\gamma}_{\rm f})$
$\displaystyle=\sigma_{0}+\sum\limits_{i=1}^{N}\sigma^{+}_{i}\left(1-e^{-t/\tau^{+}_{i}}\right),$
(2a) $\displaystyle\sigma^{-}(t;\dot{\gamma}_{\rm i},\dot{\gamma}_{\rm f})$
$\displaystyle=\sigma_{\rm
ss}+\sum\limits_{i=1}^{N}\sigma^{-}_{i}e^{-t/\tau^{-}_{i}},$ (2b)

for recovery and breakdown respectively. We assume that the system has reached
steady state during $\dot{\gamma}_{\rm i}$, and time at the step from
$\dot{\gamma}_{\rm i}$ to $\dot{\gamma}_{\rm f}$ does not matter; $t=0$ is
measured from the beginning of the $\dot{\gamma}_{\rm f}$ step.

The same idea can be applied to oscillatory shear with step strain amplitude.
From the oscillatory signal, a chosen stress amplitude component evolves as

$\displaystyle\sigma^{+}(t;\omega,\gamma_{\rm i},\gamma_{\rm f})$
$\displaystyle=\gamma_{\rm f}G_{0}(\omega)+\gamma_{\rm
f}\sum\limits_{i=1}^{N}G^{+}_{i}(\omega)\left(1-e^{-t/\tau^{+}_{i}}\right),$
(3a) $\displaystyle\sigma^{-}(t;\omega,\gamma_{\rm i},\gamma_{\rm f})$
$\displaystyle=\gamma_{\rm f}G_{\rm ss}(\omega)+\gamma_{\rm
f}\sum\limits_{i=1}^{N}G^{-}_{i}(\omega)e^{-t/\tau^{-}_{i}},$ (3b)

for recovery and breakdown respectively. The step in strain amplitude occurs
between oscillatory shear inputs $\gamma(t)=\gamma_{\rm i}\sin(\omega t)$ and
$\gamma(t)=\gamma_{\rm f}\sin(\omega t)$. Once again, we assume that the
oscillatory stress has reached steady state during $\gamma_{\rm i}$. The
stress amplitude could be total stress $\sigma$, or in-phase with strain
$\sigma^{\prime}$, or in-phase with strain rate $\sigma^{\prime\prime}$, in
which case the parameters based on $G$ would either be $|G^{*}|$, or
$G^{\prime}$ (as shown in Fig.2), or $G^{\prime\prime}$, respectively. Other
oscillatory measures also fit in this framework (e.g. orthogonal
superposition, where moduli measured in a direction orthogonal to rotational
shear, $G^{\prime}_{\perp}$ and $G^{\prime\prime}_{\perp}$, can be used with
step rate inputs [24, 4]). Other generalizations of Eq. 1 are possible. For
instance, stress in step rate tests can be normalized to define the material
function $\eta_{i}(\tau_{i})\equiv\sigma_{i}(\tau_{i})/\dot{\gamma}_{\rm f}$
and reframe Eq. 1 in terms of $\eta_{i}^{+}$ and $\eta_{i}^{-}$. We use
stresses throughout our work here to keep the analysis general and applicable
to both strain and strain rate controlled tests. The following section
elaborates the mathematical framework for obtaining the spectra.

###  II.2 Thixotropic spectra: discrete and continuous

To derive the continuous spectrum representation, consider thixotropic
recovery as an example, Eq. 1a. The transient contribution can be written in
an integral form, using the properties of the Dirac delta function, as [25]

$\displaystyle\sum\limits_{i=1}^{N}\sigma^{+}_{i}\left(1-e^{-t/\tau^{+}_{i}}\right)=\int\limits_{0}^{\infty}\left[\lim_{N\rightarrow\infty}\sum\limits_{i=1}^{N}\sigma^{+}_{i}\cdot\left(1-e^{-t/\tau^{+}}\right)\cdot\delta\left(\tau^{+}-\tau^{+}_{i}\right)\right]{\rm
d}\tau^{+}.$ (4)

Regrouping the terms, we get

$\displaystyle\sum\limits_{i=1}^{N}\sigma^{+}_{i}\left(1-e^{-t/\tau^{+}_{i}}\right)=\int\limits_{0}^{\infty}\left[\lim_{N\rightarrow\infty}\sum\limits_{i=1}^{N}\sigma^{+}_{i}\cdot\delta\left(\tau^{+}-\tau^{+}_{i}\right)\right]\left(1-e^{-t/\tau^{+}}\right){\rm
d}\tau^{+}.$ (5)

The grouping in square braces, in the limit of the number of modes
$N\rightarrow\infty$, is the continuous spectrum [22, 25]
$X^{+}\left(\tau^{+}\right)$,

$\displaystyle
X^{+}\left(\tau^{+}\right)\equiv\lim_{N\rightarrow\infty}\sum\limits_{i=1}^{N}\sigma^{+}_{i}\cdot\delta\left(\tau^{+}-\tau^{+}_{i}\right),$
(6)

where $\delta(\tau-\tau_{i})$ is the shifted Dirac delta function, with SI
units of $\rm{s}^{-1}$ [22]. $X^{+}\left(\tau^{+}\right)$ is the distribution
of thixotropic recovery modes over a domain of recovery timescales $\tau^{+}$,
each mode of spectral strength $\sigma^{+}_{i}$ and associated timescale
$\tau^{+}_{i}$, i.e. a small magnitude of _recovered_ stress on _recovery_
time. By analogy, from Eq. 1b, the continuous spectrum of breakdown timescales
is

$\displaystyle
X^{-}\left(\tau^{-}\right)\equiv\lim_{N\rightarrow\infty}\sum\limits_{i=1}^{N}\sigma^{-}_{i}\cdot\delta\left(\tau^{-}-\tau^{-}_{i}\right).$
(7)

Integration over the entire domain of $\tau$ yields the total stress from all
modes. This establishes an equivalence between the discrete and continuous
spectra, e.g. for recovery

$\displaystyle\lim_{N\rightarrow\infty}\sum\limits_{i=1}^{N}\sigma^{+}_{i}\left(1-e^{-t/\tau^{+}_{i}}\right)=\int\limits_{0}^{\infty}X^{+}\left(\tau^{+}\right)\cdot\left(1-e^{-t/\tau^{+}}\right){\rm
d}\tau^{+}.$ (8)

$X^{+}\left(\tau^{+}\right)$ is a distribution over timescales, such that an
increment in stress
$\delta\sigma^{+}=X^{+}\left(\tau^{+}\right)\delta\tau^{+}$. $X$ thus has SI
units of Pa s-1. It can also be thought of as a spectrum distributed over
logarithmically-spaced time increments, such that
$\delta\sigma^{+}=\Xi^{+}\left(\tau^{+}\right)\delta\ln\tau^{+}$, which gives

$\displaystyle X^{+}\left(\tau^{+}\right)\delta\tau^{+}$
$\displaystyle=\Xi^{+}\left(\tau^{+}\right)\delta\ln\tau^{+},$ (9a)
$\displaystyle\implies X^{+}\left(\tau^{+}\right)$
$\displaystyle=\frac{\Xi^{+}\left(\tau^{+}\right)}{\tau^{+}},$ (9b)

and $\Xi^{+}\left(\tau^{+}\right)$ has units of Pa. It is useful to recast
$X^{+}\left(\tau^{+}\right)$ as $\Xi^{+}\left(\tau^{+}\right)$ for direct
comparison between the discrete and continuous spectra plotted on log-scale
$\tau^{+}$. Note that such comparison is only possible when the discrete
spectrum $\sigma^{+}_{i}$ is log-spaced in $\tau^{+}_{i}$ [25], which we use
in this work (see the following section for details).

The equivalence between $\sigma^{+}_{i}$ (discrete) and $\Xi^{+}$ (continuous)
spectra, using Eq. 9b in Eq. 8, is

$\displaystyle\lim_{N\rightarrow\infty}\sum\limits_{i=1}^{N}\sigma^{+}_{i}\left(1-e^{-t/\tau^{+}_{i}}\right)$
$\displaystyle=\int\limits_{0}^{\infty}\frac{\Xi^{+}\left(\tau^{+}\right)}{\tau^{+}}\cdot\left(1-e^{-t/\tau^{+}}\right){\rm
d}\tau^{+}$ (10a)
$\displaystyle=\int\limits_{-\infty}^{\infty}\Xi^{+}\left(\tau^{+}\right)\cdot\left(1-e^{-t/\tau^{+}}\right){\rm
d}\ln{\tau^{+}},$ (10b)

noting the change of variables to ${\rm d}\ln\tau^{+}$ changes the lower limit
of integration. Here, we have focused on recovery with basis
$(1-e^{-t/\tau^{+}_{i}})$. Breakdown is represented in the same manner, just
by changing $(1-e^{-t/\tau^{+}_{i}})$ to $e^{-t/\tau^{-}_{i}}$ (see Table 1).

To summarize, we can now generalize the discrete spectrum of Eq. 1 as
continuous spectra. For step down in shear rate, from Eq. 2a, the stress
recovery is

$\displaystyle\sigma^{+}(t;\dot{\gamma}_{\rm i},\dot{\gamma}_{\rm
f})=\sigma_{0}+\int\limits_{-\infty}^{\infty}\Xi^{+}\left(\tau^{+};\dot{\gamma}_{\rm
i},\dot{\gamma}_{\rm f}\right)\cdot\left(1-e^{-t/\tau^{+}}\right){\rm
d}\ln\tau^{+},$ (11)

and from Eq. 2b, the stress breakdown is

$\displaystyle\sigma^{-}(t;\dot{\gamma}_{\rm i},\dot{\gamma}_{\rm
f})=\sigma_{\rm
ss}+\int\limits_{-\infty}^{\infty}\Xi^{-}\left(\tau^{-};\dot{\gamma}_{\rm
i},\dot{\gamma}_{\rm f}\right)\cdot e^{-t/\tau^{-}}{\rm d}\ln\tau^{-}.$ (12)

###  II.3 Parameterized continuous spectra

The particular shape of the continuous distribution $\Xi(\tau)$ is unspecified
in the general theory. Different shapes and parameterized continuous spectra
may be possible, similar to viscoelastic spectra, e.g. BSW, Lorentzian,
fractional Maxwell, Weibull distributions, and their numerous variations, to
name a few [26]. One may fit either discrete or continuous spectra based on
experimental observations (data fitting), or theoretical considerations from
microstructure-based models. In the past, stretched exponential forms of
stress recovery/decay have been widely used to phenomenologically describe
thixotropic transients in step shear tests [27, 16, 20, 28, 17]. Wei _et al._
, in a series of papers, used a structure kinetics based constitutive model to
study thixotropy in aggregating systems, where they employed a stretched
exponential distribution of thixotropic structure parameters over a domain of
characteristic thixotropic rate constants [20, 28]. This idea of a
distribution of structure parameters over characteristic rate constants is
similar to our approach, except that we directly describe the timescales
(rather than assuming a kinetic rate equation) and allow any shape of the
distribution (rather than assuming stretched exponential) of stress modes to
directly quantify thixotropic observations.

One could use the stretched exponential form to fit recovery/breakdown data,
and our framework enables its interpretation as a thixotropic continuous
spectrum [29, 30]. Given this common practice, we also show fits to the
stretched exponential recovery $(+)$ or decay $(-)$, given by

$\displaystyle\sigma^{+}(t)$ $\displaystyle=\sigma_{0}+\sigma^{+}_{\rm
se}\left[1-e^{-\left(t/\tau^{+}_{\rm se}\right)^{\beta}}\right],$ (13a)
$\displaystyle\sigma^{-}(t)$ $\displaystyle=\sigma_{\rm ss}+\sigma^{-}_{\rm
se}e^{-\left(t/\tau^{-}_{\rm se}\right)^{\beta}},$ (13b)

where $\sigma_{\rm se}$ is the total amount of stress change, happening over a
timescale $\tau_{\rm se}$. The stretched exponential function can be related
to an underlying continuous distribution of single-exponential modes [29, 30,
20, 28]. It can therefore be written as a continuous distribution in our
framework, given by [29, 30] (see Supplementary Information § 1 for
derivation)

$\displaystyle\Xi\left(\tau\right)$ $\displaystyle\equiv\sigma_{\rm
se}\frac{1}{\pi}\frac{\tau_{\rm
se}}{\tau}\int\limits_{0}^{\infty}e^{-u^{\beta}\cos\left(\pi\beta/2\right)}\cdot\cos\left[\frac{\tau_{\rm
se}}{\tau}u-u^{\beta}\sin\left(\frac{\pi\beta}{2}\right)\right]{\rm d}u.$ (14)

The mathematical derivation takes care to distinguish between $X(\tau)$ and
$\Xi(\tau)$. We will use $\Xi$ to compare to discrete thixotropic spectra with
log-spaced $\tau_{i}$. Note that using the stretched exponential $\Xi(\tau)$
to describe step shear data directly is _not_ equivalent to the multi-lambda
kinetic rate equation used by Wei _et al._

Among the many possible shapes of $\Xi(\tau)$, we also consider a log-normal
distribution. It is given by [25, 26]

$\displaystyle\Xi(\tau)=\Xi_{\rm
m}\exp\left[-\frac{1}{2}\left(\frac{\ln\tau-\ln\tau_{\rm
m}}{\theta}\right)^{2}\right],$ (15)

where $\Xi_{\rm m}$, $\tau_{\rm m}$, and $\theta$ are parameters of the
distribution pertaining to the peak value (strength of stress change), log-
median relaxation timescale (mean timescale of change), and standard deviation
of the spectrum (breadth of the distribution) respectively. This distribution
captures the key ideas of thixotropic spectra (a dominant timescale and
breadth of distribution) using just three parameters. It will be used to build
intuition for comparing discrete and continuous spectra.

###  II.4 Reduced parameters: moments and timescales

The moments of the discrete $(\tau_{i},\sigma_{i})$ and the continuous
distribution $\Xi\left(\tau\right)$ are defined as [22, 25]

$\displaystyle M_{n}$
$\displaystyle\equiv\sum_{i=1}^{N}\tau^{n}\cdot\sigma_{i},$ (16a)
$\displaystyle M_{n}$
$\displaystyle\equiv\int\limits_{0}^{\infty}\tau^{n}\cdot
X\left(\tau\right){\rm
d}\tau\equiv\int\limits_{-\infty}^{\infty}\tau^{n}\cdot\Xi\left(\tau\right){\rm
d}\ln\tau,$ (16b)

respectively, and so on, for $n\in\mathbb{Z}$. Based on these moments,
meaningful low-dimensional metrics, such as characteristic timescales of
recovery (or breakdown) can be defined. The general definition of the $n$-th
mean timescale of a distribution is [22, 25]

$\displaystyle\tau_{n}$ $\displaystyle\equiv\frac{M_{n}}{M_{n-1}}.$ (17)

We propose the use of Ashby-style diagrams [31] (§ V) for plotting the
quantities defined here. One meaningful quantity is $M_{0}$, the net change in
stress during the recovery/breakdown process. We denote this total change in
stress as $\Delta\sigma$,

$\displaystyle\Delta\sigma$
$\displaystyle=M_{0}\equiv\sum_{i=1}^{N}\sigma_{i},$ (18a)
$\displaystyle\Delta\sigma$
$\displaystyle=M_{0}\equiv\int\limits_{-\infty}^{\infty}\Xi\left(\tau\right){\rm
d}\ln\tau,$ (18b)

respectively, for the discrete and continuous distributions. The average
timescales $\tau_{n}$ are also meaningful. In terms of the spectra, Eq. 17 is

$\displaystyle\tau_{n}=\frac{M_{n}}{M_{n-1}}$
$\displaystyle\equiv\dfrac{\sum\limits_{i=1}^{N}\tau_{i}^{n}\cdot\sigma_{i}}{\sum\limits_{i=1}^{N}\tau_{i}^{n-1}\cdot\sigma_{i}},$
(19a) $\displaystyle\tau_{n}=\frac{M_{n}}{M_{n-1}}$
$\displaystyle\equiv\dfrac{\int\limits_{-\infty}^{\infty}\tau^{n}\cdot\Xi\left(\tau\right){\rm
d}\ln\tau}{\int\limits_{-\infty}^{\infty}\tau^{n-1}\cdot\Xi\left(\tau\right){\rm
d}\ln\tau},$ (19b)

respectively, for the discrete and continuous distributions. From this, we can
define three important quantities: the first and second mean timescales, and a
polydispersity of the timescales, as

$\displaystyle\tau_{1}$ $\displaystyle\equiv\frac{M_{1}}{M_{0}},$ (20a)
$\displaystyle\tau_{2}$ $\displaystyle\equiv\frac{M_{2}}{M_{1}},$ (20b)
$\displaystyle\rm{PDI}$
$\displaystyle\equiv\frac{\tau_{2}}{\tau_{1}}=\frac{M_{2}M_{0}}{M_{1}^{2}}.$
(20c)

There are other ways to define mean timescales. For e.g., one may define mean
timescales in the log space, such as

$\displaystyle\ln\tau_{n,{\rm log}}$
$\displaystyle\equiv\dfrac{\int\limits_{-\infty}^{\infty}(\ln\tau)^{n}\cdot\Xi\left(\tau\right){\rm
d}\ln\tau}{\int\limits_{-\infty}^{\infty}(\ln\tau)^{n-1}\cdot\Xi\left(\tau\right){\rm
d}\ln\tau}.$ (21)

In essence, the definition of mean timescales is a matter of choice. We have
defined the mean timescales in Eq. 20 to mirror the definitions used in linear
viscoelastic spectra [22, 25], and shall be using these to populate Ashby
charts for thixotropy data. Defining the mean timescales in this manner is
akin to a weighted arithmetic average of timescales, while Eq.21 is meaningful
for certain distributions such as the log-normal (Eq. 15), where mean
timescales defined on log space are equal to log-mean timescales of the
specific distribution, $\ln\tau_{n,{\rm log}}=\ln\tau_{\rm m}$.

Eq. 2a, 2b (discrete spectra) and 11, 12 (continuous spectra), in conjunction
with Eq. 18 and 19 to obtain the reduced metrics, convey the core idea of this
work. There are similarities between spectra used in viscoelasticity and the
thixotropic spectra used here. In linear viscoelasticity, the stress
relaxation shear modulus can be expressed in discrete or continuous forms as
$G(t)=\sum\limits_{i=1}^{N}G_{i}\cdot e^{-t/\tau_{i}}$, or
$G(t)=\int\limits_{0}^{\infty}Q(\tau)\cdot e^{-t/\tau}{\rm
d}\tau=\int\limits_{-\infty}^{\infty}H(\tau)\cdot e^{-t/\tau}{\rm d}\ln\tau$,
where $Q(\tau)$ and $H(\tau)$ are the modulus-weighted and viscosity-weighted
viscoelastic relaxation spectra, respectively [25]. A summary of thixotropic
spectra and the relevant basis functions, with parallels to viscoelastic
spectra, is given in Table 1.

Table 1: Thixotropic spectra with parallels to viscoelastic spectra. Definitions of $Q(\tau)$ and $H(\tau)$ following Martinetti _et al._ [22, 25, 26] Property | Data | Basis function | Interpretation | Discrete | Continuous (on ${\rm d}\tau$)
---|---|---|---|---|---
Thixo, $\dot{\gamma}_{\rm f}>\dot{\gamma}_{\rm i}$ | $\sigma^{-}(t;\dot{\gamma}_{\rm i},\dot{\gamma}_{\rm f})$ | $\sigma^{-}_{i}e^{-t/\tau^{-}_{i}}$ | Breakdown | $\sigma^{-}_{i}(\tau^{-}_{i})$ | $X^{-}(\tau^{-})$
Thixo, $\dot{\gamma}_{\rm f}<\dot{\gamma}_{\rm i}$ | $\sigma^{+}(t;\dot{\gamma}_{\rm i},\dot{\gamma}_{\rm f})$ | $\sigma^{+}_{i}\left(1-e^{-t/\tau^{+}_{i}}\right)$ | Recovery | $\sigma^{+}_{i}(\tau^{+}_{i})$ | $X^{+}(\tau^{+})$
Viscoelastic[32] | $G(t;\gamma_{0})$ | $G_{i}e^{-t/\tau_{i}}$ | Decay | $G_{i}(\tau_{i})$ | $Q(\tau)$
(relaxation) | $\eta^{+}(t;\dot{\gamma}_{0})$ | $\eta_{i}\left(1-e^{-t/\tau_{i}}\right)$ | Growth | $\eta_{i}(\tau_{i})=\tau_{i}G_{i}(\tau_{i})$ | $H(\tau)=\tau Q(\tau)$

The reduced parameters from using thixotropic spectra, $\Delta\sigma$ and
$\tau_{n}$, obtained from Eq. 17, 18, can be used to describe any thixotropic
data or constitutive model response. Using $\Delta\sigma$ and $\tau_{n}$, we
can now populate Fig. 1 to report and compare the degree of thixotropy in
different materials. Such results are presented in § V. The application of
this framework to step shear data is illustrated through a theoretical example
in the following section.

###  II.5 Illustration: relation between continuous, discrete, and reduced
metrics

(a)

(b)

Figure 3: Log-normal thixotropic recovery spectrum compared to discrete
spectrum approximations. Simulated data with $\sigma_{0}=10$ Pa, $\Xi_{\rm
m}=10$ Pa, $\tau_{\rm m}=1$ s, $\theta=1$, in Eq. 15. (a) Stress response to
step down in shear rate, (b) co-plot of each spectrum with important reduced
metrics shown.

Using the log-normal $\Xi^{+}(\tau^{+})$ of Eq. 15 with parameters $\Xi_{\rm
m}=10$ Pa, $\tau_{\rm m}=1$ s, $\theta=1$, and $\sigma_{0}=10$ Pa, in Eq. 11,
we can produce the typical stress recovery shape in Fig. 3(a). This is similar
to the schematic in Fig. 2(a2). The shape of $\Xi^{+}(\tau^{+})$ is shown as
the solid line in Fig. 3(b). The log-median timescale $\tau_{\rm m}=1$ s is
associated with the transients shown in Fig. 3(a), while the peak $\Xi_{\rm
m}=10$ Pa is related to the amount of stress change. The PDI of the timescales
is related to $\theta$, the spread of the log-normal spectrum [25]. From
$\Xi^{+}$, we obtain the generalized low-dimensional metrics, namely the
moments and characteristic timescales, using Eq. 18 and 19. The timescales
$\tau^{+}_{1}$ and $\tau^{+}_{2}$ for the spectra are shown on the plot. Note
that $\tau^{+}_{2}>\tau^{+}_{1}$. Also note that $\tau^{+}_{1}=1.6~{}{\rm
s}>\tau_{\rm m}=1~{}{\rm s}$. This is because we calculate mean timescales in
$\tau$ space, and not $\ln\tau$ space. If we were to use Eq. 21, we would get
$\tau^{+}_{1,{\rm log}}=\tau_{\rm m}$. But, as mentioned earlier, this is a
matter of choice, and we use mean timescales defined in $\tau$ space similar
to those in viscoelasticity.

We also show high-density discrete, and low-density discrete distributions,
each co-plotted in Fig. 3(b). The high density (50 modes) and low density (5
modes) distributions $(\sigma^{+}_{i})$ are discrete representations of
$\Xi^{+}$ (obtained by fitting $\sigma^{+}_{i}$ over a set domain of
$\tau^{+}_{i}$ to match the $\sigma^{+}(t)$ in Fig. 3(a) resulting from the
continuous distribution). Also note that $\sigma^{+}_{i}$ lies below $\Xi^{+}$
for a sufficiently large number of modes, whereas for fewer modes the values
of $\sigma^{+}_{i}$ lie above $\Xi^{+}$. This is a consequence of the fact
that each distribution must result in the same stress signal when integrated
(or summed). As more modes are added, the value of each mode must decrease in
order to add up to the same stress response. The plot also shows a single data
point $(\tau^{+}_{1},\Delta\sigma^{+})$, which is the pair of reduced or
summarizing metrics obtained using Eq. 18 and 19, in essence is a single-mode
representation of the polydisperse system. This can be plotted in the
thixotropic Ashby charts. It lies above all other distributions in Fig. 3(b)
as expected. This can be visualized as a single-mode stress response
$\sigma^{+}(t)=\sigma_{0}+\Delta\sigma^{+}\left(1-e^{-t/\tau^{+}_{1}}\right)$,
shown in Fig. 3(a). This clearly shows that $\tau^{+}_{1}\neq\tau_{\rm m}$,
and $\tau^{+}_{1}$ is indeed a mean timescale of the original stress response
since the $(\tau^{+}_{1},\Delta\sigma^{+})$ mode goes through the more complex
$\sigma^{+}(t)$ response. It must be noted that the stress generated from
$(\tau^{+}_{1},\Delta\sigma^{+})$ is not equivalent to a single-mode discrete
spectrum fit of $\sigma^{+}(t)$ data.

The illustrative example in Fig. 3 provides intuition for the breadth of a
thixotropic spectrum. Here, ${\rm PDI}=\tau^{+}_{2}/\tau^{+}_{1}=2.08$,
quantifying the deviation from a single timescale. It shows that discrete
spectra can provide excellent descriptions of $\sigma^{+}(t)$, and that the
log-spaced shape of $\sigma_{i}(\tau_{i})$ will be similar to the shape of any
underlying $\Xi(\tau)$ representation (if $\tau_{i}$ are equally spaced in
$\ln\tau_{i}$). In what follows with experimental data, we use discrete
spectra, log-spaced in $\tau_{i}$, to fit $\sigma^{+}(t)$ data and obtain
reduced metrics on a number of thixotropic materials. The materials used and
rheometry and fitting methods employed are shown in the next section, and the
fitting results are shown in § IV.

##  III Experimental: materials and methods

###  III.1 Materials

We test two common yield-stress fluids: 4.0 wt% suspension of Laponite RD
colloidal clay particles in water, which forms a sparse physical gel when
dispersed, typically identified as a thixotropic material [16]; and 1.0 wt%
suspension of Carbopol 940 microgel particles in water, a model yield-stress
fluid formed of jammed, swollen, crosslinked polyacrylic acid microgel
particles, known to exhibit negligible thixotropy [33, 34]. These suspensions
were prepared using methods from the literature [35, 36]. We also use step
shear data for two model thixotropic yield-stress fluids from the open
literature: 3.23 vol% carbon black suspensions in naftenic oil [27], and 2.9
vol% fumed silica suspension in paraffin oil with PIB [37]. Both form
thixotropic colloidal dispersions when suspended in their respective media. We
show reduced metrics for all the above materials in the Ashby diagrams, but
only show detailed step shear rate data for Laponite in the next section for
the sake of brevity.

###  III.2 Rheometry methods

The data for all materials is either a step down from an initial shear rate of
$\dot{\gamma}_{\rm i}=5.00~{}\rm{s}^{-1}$ to various lower final shear rates
$\dot{\gamma}_{\rm f}=0.25,~{}0.50,~{}1.00,~{}2.50~{}{\rm s}^{-1}$, or a step
up from $\dot{\gamma}_{\rm i}=0.10~{}\rm{s}^{-1}$ to various higher shear
rates $\dot{\gamma}_{\rm f}=0.50,~{}1.00,~{}2.50,~{}5.00~{}{\rm s}^{-1}$ (as
shown in Fig. 2(a) and (b) respectively). Data for fumed silica and carbon
black were used directly from the published literature. For Laponite and
Carbopol, the materials were sheared between 600 grit sandpaper-covered
parallel plate geometries (uncorrected for parallel plate shear rate
inhomogeneity), 25 mm in diameter, at two different gaps (750 and 1000
$\mu{\rm m}$) to verify the absence of wall slip [38], at $25^{\circ}$C, on an
ARES-G2 strain-controlled, separated motor-transducer rheometer (TA
Instruments).

###  III.3 Fitting methods

The model fits were done using iterative linear regression algorithms that in
their most general form minimize a weighted residual sum of squares, ${\rm
RSS}=\sum\limits_{i=1}^{d}\left[y_{i}-f(x_{i})\right]^{2}/w_{i}$, where
$y_{i}$ are the experimental data (shear stress $\sigma(t)$) and $f(x_{i})$
are the model predictions over the control variable (time $t$). The residuals
were weighted by the experimental data, so that $w_{i}=y_{i}$, which
implicitly assumes an error proportional to the data [39]. For the stretched
exponential fits, the data is fit to Eq. 13a, and Eq. 14 is used to interpret
the fit as a continuous spectrum. For the discrete spectrum fits, obtaining
the underlying distribution is an ill-posed, inverse problem, and we employed
Tikhonov regularization [40, 41, 42, 28] to obtain the optimal set of
$\sigma_{i}$ for a pre-distributed, log-spaced set of $\tau_{i}$ values (using
open-source MATLAB scripts available online [43]).

##  IV Spectra from experimental data

We first show results for step up tests, followed by step down tests. We co-
plot the fits to the discrete spectra and the stretched exponential for the
transient stress data, and hence also co-plot the discrete and continuous
spectra (corresponding to the underlying spectra for the stretched
exponential). The fit parameters for the stretched exponential for each
dataset ($\sigma_{\rm se}$ and $\tau_{\rm se}$) are shown in Supplementary
Information, § 2.

###  IV.1 Step up in shear rate

Fig. 4 shows step up in shear rate data for 4.0 wt% Laponite suspension in
water, where (a) shows the transient stress data, and fits to the discrete
spectrum (Eq. 2b, solid lines) and the stretched exponential (Eq. 13a, dashed
lines). The data in (a) suggest that the timescales of breakdown are
$\mathcal{O}(10)~{}{\rm s}$, and the amount of stress change increases from
$\approx 20$ to $100~{}{\rm Pa}$ as $\dot{\gamma}_{\rm f}$ increases. The
effect of thixotropy is most extreme for $\dot{\gamma}_{\rm f}=5.00~{}{\rm
s}^{-1}$.

(a)

(b)

Figure 4: Step up in shear rate data for 4.0 wt% Laponite suspension in water,
with $\dot{\gamma}_{\rm i}=0.1~{}\rm{s}^{-1}$. (a) Transient stress decay fit
to a discrete thixotropic spectrum (solid lines) and to a stretched
exponential (dashed lines) form of breakdown. (b) Co-plots of the discrete
spectrum (circles) and the continuous spectrum interpretation of the stretched
exponential (solid lines) for each value of $\dot{\gamma}_{\rm f}$. The
$\bigstar$ shows the summarizing metric for the discrete spectrum at the most
extreme case of thixotropy, corresponding to $\dot{\gamma}_{\rm
f}=5.00~{}\rm{s}^{-1}$.

Thixotropic spectra help quantify these observations. The existence of
multiple modes of breakdown, evidenced by the changing concavity in the
$\sigma^{-}(t)$ data, is captured well by the discrete spectrum. The multi-
modal nature of the thixotropic breakdown is exemplified in Fig. 4(b), where
the discrete spectrum $\sigma^{-}_{i}(\tau^{-}_{i})$ shows multiple peaks at
distinct timescales, $\tau^{-}_{i}$. This feature is missing from
$\Xi^{-}\left(\tau^{-}\right)$ for the single-mode stretched exponential. The
summarizing metrics of the discrete spectra follow a systematic trend; both
timescales and amount of stress change increase as $\dot{\gamma}_{\rm f}$
increases from 0.50 to 5.00 s-1 ($\tau^{-}_{1}=1.06,~{}2.23,~{}2.56,~{}3.91$
s, and $\Delta\sigma^{-}=18.52,~{}32.64,~{}64.77,~{}102.78$ Pa). The degree of
thixotropy therefore is the most extreme for the highest $\dot{\gamma}_{\rm
f}$, and the corresponding values of $(\tau^{-}_{1},\Delta\sigma^{-})$ are co-
plotted with the spectra in Fig. 4(b). Obtaining such simple yet meaningful
descriptions is the key result of this work.

Similar features are retained for step down in shear data as shown in the next
section, although the evidence is not conclusive if one only looks at the
stress data; one must also look at the underlying spectrum to truly visualize
the distribution of characteristic timescales during transient thixotropic
processes.

###  IV.2 Step down in shear rate

(a)

(b)

Figure 5: Step down in shear rate data for 4.0 wt% Laponite suspension in
water, with $\dot{\gamma}_{\rm i}=5.0~{}\rm{s}^{-1}$. (a) Transient stress
growth fit to a discrete thixotropic spectrum (solid lines) and to a stretched
exponential (dashed lines) form of buildup. (b) Co-plots of the discrete
spectrum (circles) and the continuous spectrum interpretation of the stretched
exponential (solid lines) for each value of $\dot{\gamma}_{\rm f}$. The
$\bigstar$ shows the summarizing metric for the discrete spectrum at the most
extreme case of thixotropy, corresponding to $\dot{\gamma}_{\rm
f}=0.25~{}\rm{s}^{-1}$.

Fig. 5 shows step down in shear data for 4.0 wt% Laponite suspension in water,
where (a) shows the transient stress data, and fits to the discrete spectrum
(Eq. 2a, solid lines) and the stretched exponential (Eq. 13a, dashed lines).
The difference between the discrete _versus_ continuous distributions is not
apparent from the stress data; only from (b) do we see the existence of
multiple modes of recovery, which a single-timescale stretched exponential
function cannot capture. As for the step up data, the multi-modal nature of
the thixotropic recovery is exemplified in Fig. 5(b), which is captured by
$\sigma^{+}_{i}(\tau^{+}_{i})$ but missed by $\Xi^{+}\left(\tau^{+}\right)$.

The effect of $\dot{\gamma}_{\rm f}$ is more prominent for the recovery
spectra, where the peaks in $\sigma^{+}_{i}$ shift toward longer timescales,
indicating that the stress recovery takes longer at lower shear rates. This
trend is also observed in the stretched exponential $\Xi^{+}$. This is
reflected in the summarizing metrics, which again follow a systematic trend;
both timescale and amount of stress recovered increase as $\dot{\gamma}_{\rm
f}$ decreases from 2.50 to 0.25 s-1
($\tau_{1}^{+}=2.84,~{}5.20,~{}5.21,~{}6.14$ s, and
$\Delta\sigma^{+}=18.99,~{}33.63,~{}36.83,~{}42.38$ Pa). The degree of
thixotropy is the most extreme for the lowest $\dot{\gamma}_{\rm f}$, and the
corresponding values of $(\tau_{1}^{+},\Delta\sigma^{+})$ are co-plotted with
the spectra in Fig. 5(b). Such a monotonic trend of increasing recovery
timescales with decrease in shear rate is also predicted by classical
thixotropic constitutive equations [44, 45, 14, 17] (see Supplementary
Information, § 3 for an example), and it is indeed a satisfying if not
surprising result that the same trend exists in this dataset.

Note that during any thixotropic process, breakdown or recovery, both shear-
induced hydrodynamic forces/stresses (which tend to break up particle
aggregates) and Brownian motion due to thermal fluctuations (which helps
particles find minima in the interparticle potential energy landscape and
assist floc growth) contribute to the net resultant transient stress
evolution, i.e. the $\sigma(t)$ data shown in the plots. The relative strength
of these two effects is often compared using the Péclet number [18], defined
as

$\displaystyle{\rm Pe}\equiv\frac{6\pi\eta_{\rm m}\dot{\gamma}_{\rm
c}a^{3}}{k_{\rm B}T}\sim\frac{\text{advective transport rate}}{\text{diffusive
transport rate}},$ (22)

where $\eta_{\rm m}$ is the viscosity of the suspending medium, $a$ is the
primary particle size, and $T$ is the absolute temperature of the system. For
our step rate tests, the characteristic rate of shear in the system,
$\dot{\gamma}_{\rm c}=\dot{\gamma}_{\rm f}$. The strength of Brownian
hydrodynamics is independent of shear rate. At high $\dot{\gamma}_{\rm f}$,
shear-induced stresses dominate, and the effect of intrinsic particle
interactions and Brownian motion is limited. This is the case during step up
tests, and advective hydrodynamic stresses dictate the rate of structure
breakdown. On the other hand, in step down tests, the thixotropic structure
buildup is primarily influenced by interparticle forces and Brownian motion,
and less so by advective hydrodynamics, at very low values of
$\dot{\gamma}_{\rm f}$ such that $\rm Pe\ll 1$. The resultant recovery
processes take much longer and the thixotropic transients are more easily
observable. We therefore choose to compare the thixotropic properties of
Laponite to other fluids by looking at recovery in step down in shear data in
§ V. We show step down in shear data for other materials in the next
subsection.

Also note that this limit of $\rm Pe\ll 1$ for $\dot{\gamma}_{\rm
f}\rightarrow 0$ is a key feature to report, because it represents the most
extreme limit of thixotropic recovery. Using the simplest, classical
thixotropic constitutive equations [44, 45, 14], each timescale mode
$\tau^{+}_{i}$ can be mapped to a kinetic aggregation rate constant in this
limit, which helps connect our descriptions directly with physical quantities
used in thixotropic constitutive modeling. The details are shown in
Supplementary Information, § 3.

###  IV.3 Step down in shear rate for three model materials

(a) carbon black  (b) fumed silica  (c) Carbopol
---

---

Figure 6: Step down in shear data for (a) 3.23 vol% carbon black (from
Dullaert _et al._ [27]), (b) 2.9 vol% fumed silica (from Armstrong _et al._
[37]), (c) 1 wt% Carbopol. The top row shows the transient stress in recovery
for each material, the fit lines correspond to the discrete (solid lines) and
stretched exponential continuous (dashed lines) spectra. The bottom row shows
the corresponding discrete (circles) and stretched exponential continuous
(solid lines) spectra, similar to the data shown for Laponite in Fig. 5. The
$\bigstar$ in each spectrum plot shows the summarizing metric for the most
extreme case of thixotropy, corresponding to $\dot{\gamma}_{\rm
f}=0.25~{}\rm{s}^{-1}$.

In addition to Laponite, we analyzed thixotropic recovery data for three other
materials listed in § III.1. Fig. 6 shows the data for each material as
separate columns, with the transient stress data and fit spectra for each,
similar to Fig. 5. The features observed for both the stretched exponential
$\Xi^{+}$ and discrete $\sigma^{+}_{i}$ in recovery for Laponite are similarly
observed for these materials. The thixotropic timescales range from $<1~{}{\rm
s}$ (e.g. carbon black at $\dot{\gamma}_{\rm f}=2.5~{}\rm{s}^{-1}$) to
$>10~{}{\rm s}$ (e.g. fumed silica at $\dot{\gamma}_{\rm
f}=0.25~{}\rm{s}^{-1}$). For carbon black and silica, the amount of stress
change relative to the steady state value increases as $\dot{\gamma}_{\rm f}$
decreases, while this ratio is much smaller for Carbopol at all shear rates;
this is quantified with $\Delta\sigma^{+}/\sigma_{\infty}$, shown in Fig. 7.
We also observe the existence of distinct peaks in the spectra, in addition to
a spread in timescales, for Laponite, carbon black, and Carbopol, whereas we
observe a broad spectrum and an absence of distinct, separate peaks in the
fumed silica data. The polydispersity of thixotropic modes is captured in the
spectral method developed here, including the ability to see distinct groups
of modes via peaks in the distributions. The reduced summarizing metrics for
Fig. 5, 6 are shown in the following section.

##  V Low-dimensional descriptions from spectra

Fig. 7 is an Ashby-style co-plot [31] of reduced, low-dimensional metrics of
thixotropic recovery spectra from Figs. 5 and 6. This allows direct comparison
between the four different materials using $\Delta\sigma^{+}/\sigma_{\infty}$
_versus_ $\tau^{+}_{1}$, where $\Delta\sigma^{+}$ is the total stress
recovered under $\dot{\gamma}_{\rm f}$, and is obtained from Eq. 18. The
recovered stress is normalized by the steady state stress at long times,
defined as $\sigma_{\infty}\equiv\lim_{t\rightarrow\infty}\sigma(t)$. Since
$\sigma_{\infty}$ is used as the reference stress and the recovery is
quantified as a percentage of this value, we get

$\displaystyle 0\leq\frac{\Delta\sigma^{+}}{\sigma_{\infty}}\leq 1,$ (23)

and this allows for comparison of thixotropic recovery for different
materials. $\tau^{+}_{1}$ is the first mean timescale of recovery, obtained by
using $n=1$ in Eq. 19.

(a)

(b)

Figure 7: Ashby co-plots for thixotropic recovery, based on data in Figs. 5,
6. (a) Reduced metrics $(\tau^{+}_{1},\Delta\sigma^{+}/\sigma_{\infty})$ are
obtained from discrete spectrum fits to step down in shear data for different
fluids. The symbols get smaller as the final shear rate decreases
($\dot{\gamma}_{\rm f}=2.50,~{}1.00,~{}0.50,~{}0.25~{}{\rm s}^{-1}$) for each
material. (b) The mean timescales $\tau^{+}_{1},\tau^{+}_{2}$ (visualizing
${\rm PDI}=\tau^{+}_{2}/\tau^{+}_{1}$ for the most thixotropic response for
each material (observed at $\dot{\gamma}_{\rm f}=0.25~{}{\rm s}^{-1}$) plotted
against $\Delta\sigma^{+}/\sigma_{\infty}$.

The closer the value of $\Delta\sigma^{+}/\sigma_{\infty}$ is to 1, the
greater is the change in stress due to thixotropic recovery relative to the
final state of the sample, and the greater the significance of thixotropic
effects in the material for a given set of $\left(\dot{\gamma}_{\rm
i},\dot{\gamma}_{\rm f}\right)$. Among these data, there is a single value of
$\dot{\gamma}_{\rm i}=5~{}{\rm s}^{-1}$, and only $\dot{\gamma}_{\rm f}$ is
varied. Accordingly, we see the effect of changing the shear history in Fig.
7(a), where smaller symbols correspond to a lower $\dot{\gamma}_{\rm f}$,
approaching the Brownian-dominated recovery limit. As $\dot{\gamma}_{\rm f}$
is decreased, $\Delta\sigma^{+}/\sigma_{\infty}$ increases, as observed in
Fig. 7, and hence the degree of thixotropy also increases in terms of the
stress recovered. But as mentioned in § I, one must consider both stresses and
timescales to judge the degree of thixotropy in a material.

The longer the timescale $\tau^{+}_{1}$, the longer the thixotropic transients
last, and the greater the significance of thixotropic effects in the material,
once again for a given set of $\left(\dot{\gamma}_{\rm i},\dot{\gamma}_{\rm
f}\right)$. The same trend with changing $\dot{\gamma}_{\rm f}$ is observed in
$\tau^{+}_{1}$ as with $\Delta\sigma^{+}/\sigma_{\infty}$: as
$\dot{\gamma}_{\rm f}$ is decreased, $\tau^{+}_{1}$ increases, and thus the
degree of thixotropy is most significant at the lowest $\dot{\gamma}_{\rm f}$
(Fig. 5, Fig. 6). This is observed for all four materials shown in Fig. 7(a).
We see that the timescales are comparable across all materials ($\approx 1-10$
s), but the fractional stress recovery is very different. Considering the
lowest $\dot{\gamma}_{\rm f}$ for each, Carbopol has
$\Delta\sigma^{+}/\sigma_{\infty}=0.09$, and is therefore the least
thixotropic, even though its $\tau^{+}_{1}=5.41$ s is similar to that of the
other materials. Laponite is more thixotropic than Carbopol, while carbon
black and fumed silica are the most thixotropic.

The trends in Fig. 7(a) form a compact summary and are consistent with basic
thixotropic constitutive models. As mentioned in § IV, the monotonic trend in
$\tau^{+}_{1}$ with $\dot{\gamma}_{\rm f}$ is in accordance to classical
structural kinetics-based thixotropic models in the literature [16, 17] (see
Supplementary Material, § 3). Structure recovery through Brownian motion
proceeds slowly, and at lower $\dot{\gamma}_{\rm f}$, the effect of Brownian
dynamics is more prominent. The increase in $\Delta\sigma^{+}/\sigma_{\infty}$
with decreasing $\dot{\gamma}_{\rm f}$ also is physically meaningful; as
$\dot{\gamma}_{\rm f}$ decreases, $\rm Pe$ decreases, and Brownian dynamics
dominate over advective hydrodynamics, and more structure can build up (or
equivalently, less structure is broken down due to shear-induced stresses). It
is interesting that the two quantities used here qualitatively depend on
$\dot{\gamma}_{\rm f}$ in the same way. The Ashby plots are simple to
interpret and use: the further away a material lies from the origin, the more
significant is the effect of thixotropy in the material (moving from the gray
to the white region).

Fig. 7(a) does not indicate the breadth of underlying thixotropic spectra, but
showing both $\tau^{+}_{1}$ and $\tau^{+}_{2}$, as in Fig. 7(b) for the
minimum $\dot{\gamma}_{\rm f}$ data, reveals this. The more distinct
$\tau^{+}_{1}$ and $\tau^{+}_{2}$ are, the more polydisperse the thixotropic
recovery spectrum, ${\rm PDI}\equiv\tau^{+}_{2}/\tau^{+}_{1}$ from Eq. 20c
(${\rm PDI}=3.31,~{}1.77,~{}1.49,~{}1.41$ for Carbopol, Laponite, carbon
black, and fumed silica respectively). This representation in Fig. 7(b) is
general, a description applicable to any material, even those with different
underlying constitutive behavior. It is therefore a powerful tool that is
useful for comparison and selection across different microstructural classes
of thixotropic materials. This technique has already been applied by the
authors in studying drop impact dynamics of thixotropic yield-stress fluids
[36, 46], where reduced thixotropic metrics obtained from step shear data were
used to select a material (Laponite) with a long recovery timescale and
significant stress changes which has important consequences in their drop
impact stick-splash outcomes [46].

##  VI Conclusions

In this work, we have developed a general framework for thixotropic spectra to
quantitatively describe thixotropy using a distribution of recovery/breakdown
timescales in step tests. The moments and means of the distributions are used
to generate low-dimensional metrics of quantifying the degree of thixotropy in
the material, and these descriptions are both model-independent and material
agnostic. Alternative reduced metrics for quantifying thixotropic timescales
have been used by Divoux _et al._ [34] and Jamali _et al._ [47] using areas of
thixotropic hysteresis loops. However, the methods developed in our work
employ faster experimental protocols, largely isolate the effect of
viscoelastic transients from thixotropic ones, can be applied to not only
unidirectional but also oscillatory shear tests, and most importantly,
separately describe recovery and breakdown and reveal spectral distributions.
We mention the analogy of the methods developed here with viscoelastic
relaxation spectra in Table 1. Similarly, an analogy may be drawn between
stress controlled step tests where the strain rate response is expressed
analogous to retardation spectra distributed over retardation timescales.

The methodology proposed in this work is not free from limitations,
experimental or analytical. Rheometer motors that impose strains/strain rates
on samples via angular displacements have a finite response time to ramp up or
down $(t_{\rm mot})$, and step changes are not instantaneous. The shortest
duration of change between $\dot{\gamma}_{\rm i}$ and $\dot{\gamma}_{\rm f}$
is limited by the inertia and response time of the rotor-geometry assembly.
The quickest step changes are essentially ramps that at best occur over
$t_{\rm mot}\sim\mathcal{O}(10)$ ms, and therefore transients in
$\sigma(t;\dot{\gamma}_{\rm i},\dot{\gamma}_{\rm f})$ for $t\lesssim t_{\rm
mot}$ cannot be observed. The credibility of short-time data is also limited
by viscoelastic waves in the sample [38]. Consequently, the range of $\tau$
and $\tau_{i}$ used for fitting spectra to stress data has to be truncated.
The temporal location and value of $\sigma_{0}$ is thus also influenced by
$t_{\rm mot}$. This in turn makes the value of $\Delta\sigma/\sigma_{\infty}$
and therefore the reduced descriptions used in the Ashby charts susceptible to
instrument limitations used to collect the data.

Specific details of the shape of the distribution are lost upon reduction into
its moments. It is possible to generate completely different distributions and
stress responses from the same set of moments, one prominent example being
Anscombe’s quartet [48]. This also typifies the ill-posed nature of the
fitting problem. Multiple distributions can give the same stress response, and
therefore the same means and moments. Inverting the summation in Eq. 1 to
obtain $\sigma_{i}$ is an ill-posed problem, similar to those faced with
viscoelastic spectra. The approach of using a specific analytically defined
continuous spectrum to fit to stress data and obtain $\Xi(\tau)$ may be a
workaround, especially when employed with Bayesian inference methods and
credibility metrics, as has been applied to LVE [49], MAOS [25], and
thixotropic constitutive model selection [37].

Future work will benefit from this methodology of thixotropic spectra and
Ashby diagrams. The high-density discrete spectra $\sigma_{i}$ can be employed
to obtain dominant timescale modes, and a low-density superposition of
continuous spectra centered around these timescales could be used to fit the
data to very high accuracy. The method developed here is not predictive, but
only descriptive. Yet, the value of descriptive methods cannot be overstated,
e.g. material properties which are extremely useful in rheological analysis
despite not being predictive. Our method also gives insight into thixotropic
properties of fluids that might inform predictive models; for instance,
revealing or suggesting shapes of underlying spectra for multi-mode
constitutive equations, which can lead to better predictions. Shapes of the
discrete spectra can also be used to inform the shape of parameterized
continuous distributions of $\Xi(\tau)$. These make the fitting process less
computationally intensive. The concept of thixotropic spectra also provides
interpretation to commonly-used fit functions, such as the stretched
exponential (Eq. 13a, 14). Lastly, the methodology developed here can be
applied to any material system that evolves with time to obtain summarizing
metrics that describe the data in a simple yet meaningful way. We now have a
way to judge “how thixotropic” a material is in terms of both timescale and
amount of change, which is important for dimensionless groups involving
thixotropic and aging systems in general (e.g. the Mutation number [50, 51,
52]), and for rheological design requirements involving thixotropy [53] in
applications such as 3D printing [54], spray coating, fire suppression [55,
56], flow batteries [1, 2, 3], food processing, and other uses not yet
imagined.

##  VII Supplementary Material

See Supplementary Information for more information about the continuous
spectrum for the stretched exponential and the parameters for the step down
data fits for all materials to the stretched exponential model.

##  VIII Acknowledgements

This project was funded by the National Science Foundation CAREER Award,
CBET-1351342, and the Joint Center for Energy Storage Research (JCESR), U.S.
Department of Energy. S.S. thanks Y. Wang, N. Ramlawi, and Dr. C. Saengow for
helpful discussions.

## References

  * Fan  _et al._ [2014] F. Y. Fan, W. H. Woodford, Z. Li, N. Baram, K. C. Smith, A. Helal, G. H. McKinley, W. C. Carter, and Y. M. Chiang, “Polysulfide flow batteries enabled by percolating nanoscale conductor networks,” Nano Lett. 14, 2210–2218 (2014).
  * Helal, Divoux, and McKinley [2016] A. Helal, T. Divoux, and G. H. McKinley, “Simultaneous rheoelectric measurements of strongly conductive complex fluids,” Phys. Rev. Appl. 6, 1–19 (2016).
  * Narayanan, Mugele, and Duits [2017] A. Narayanan, F. Mugele, and M. H. G. Duits, “Mechanical history dependence in carbon black suspensions for flow batteries: A rheo-impedance study,” Langmuir 33, 1629–1638 (2017).
  * [4] Y. Wang and R. H. Ewoldt, “New insights on carbon black suspension rheology - anisotropic thixotropy and anti-thixotropy,” Journal of Rheology, submitted .
  * Dimitriou and Mckinley [2014] C. J. Dimitriou and G. H. Mckinley, “A comprehensive constitutive law for waxy crude oil: A thixotropic yield stress fluid,” Soft Matter 10, 6619 (2014).
  * Glicerina  _et al._ [2015] V. Glicerina, F. Balestra, M. D. Rosa, and S. Romani, “Microstructural and rheological properties of white chocolate during processing,” Food Bioprocess Technol. 8, 770–776 (2015).
  * Jin  _et al._ [2011] L. Jin, Z. Wu, T. Wei, J. Zhai, and X. Zhang, “Dye-sensitized solar cell based on blood mimetic thixotropy sol–gel electrolyte,” Chem. Commun. 47, 997–999 (2011).
  * Armstrong and Tussing [2020] M. J. Armstrong and J. Tussing, “A methodology for adding thixotropy to Oldroyd-8 family of viscoelastic models for characterization of human blood,” Phys. Fluids 32, 094111 (2020).
  * Dullaert and Mewis [2005] K. Dullaert and J. Mewis, “A model system for thixotropy studies,” Rheol. Acta 45, 23–32 (2005).
  * Alessandrini, Lapasin, and Sturzi [1982] A. Alessandrini, R. Lapasin, and F. Sturzi, “The kinetics of thixotropic behaviour in clay/kaolin aqueous suspensions,” Chem. Eng. Commun. 17, 13–22 (1982).
  * Beris, Stiakakis, and Vlassopoulos [2008] A. N. Beris, E. Stiakakis, and D. Vlassopoulos, “A thermodynamically consistent model for the thixotropic behavior of concentrated star polymer suspensions,” J. Non-Newtonian Fluid Mech. 152, 76–85 (2008).
  * Burgos, Alexandrou, and Entov [2001] G. R. Burgos, A. N. Alexandrou, and V. Entov, “Thixotropic rheology of semisolid metal suspensions,” J. Mater. Process. Technol. 110, 164–176 (2001).
  * Kelessidis [2008] V. Kelessidis, “Investigations on the thixotropy of bentonite suspensions,” Energy Sources A 30, 1729–1746 (2008).
  * Mewis [1979] J. Mewis, “Thixotropy – a general review,” J. Non-Newtonian Fluid Mech. 6, 1–20 (1979).
  * Barnes [1997] H. A. Barnes, “Thixotropy - a review,” J. Non-Newtonian Fluid Mech. 70, 1–33 (1997).
  * Mewis and Wagner [2009] J. Mewis and N. J. Wagner, “Thixotropy,” Adv. Colloid Interface Sci. 147–148, 214–227 (2009).
  * Larson and Wei [2019] R. G. Larson and Y. Wei, “A review of thixotropy and its rheological modeling,” J. Rheol. 63, 477–501 (2019).
  * Mewis and Wagner [2012] J. Mewis and N. J. Wagner,  _Colloidal Suspension Rheology_ (Cambridge University Press, 2012).
  * Joshi [2021] Y. M. Joshi, “Thixotropy, nonmonotonic stress relaxation, and the second law of thermodynamics,” J. Rheol. 66, 111–123 (2021).
  * Wei, Solomon, and Larson [2016] Y. Wei, M. J. Solomon, and R. G. Larson, “Quantitative nonlinear thixotropic model with stretched exponential response in transient shear flows,” J. Rheol. 60, 1301 (2016).
  * Larson [2015] R. G. Larson, “Constitutive equations for thixotropic fluids,” J. Rheol. 59, 595 (2015).
  * Tschoegl [1989] N. W. Tschoegl,  _The Phenomenological Theory of Linear Viscoelastic Behavior_ (Springer-Verlag, Berlin, 1989).
  * Vlad  _et al._ [1996] M. O. Vlad, R. Metzler, T. F. Nonnenmacher, and M. C. Mackey, “Universality classes for asymptotic behavior of relaxation processes in systems with dynamical disorder: Dynamical generalizations of stretched exponential,” Energy Sources A 37, 2279 (1996).
  * Vermant  _et al._ [1998] J. Vermant, P. Walker, P. Moldenaers, and J. Mewis, “Orthogonal versus parallel superposition measurements.” J. Non-Newtonian Fluid Mech. 79, 173–189 (1998).
  * Martinetti, Soulages, and Ewoldt [2018] L. Martinetti, J. M. Soulages, and R. H. Ewoldt, “Continuous relaxation spectra for constitutive models in medium-amplitude oscillatory shear,” J. Rheol. 62, 1271–1298 (2018).
  * Martinetti and Ewoldt [2019] L. Martinetti and R. H. Ewoldt, “Time-strain separability in medium amplitude oscillatory shear,” Phys. Fluids 31, 021213 (2019).
  * Dullaert and Mewis [2006] K. Dullaert and J. Mewis, “A structural kinetics model for thixotropy,” J. Non-Newtonian Fluid Mech. 139, 21–30 (2006).
  * Wei, Solomon, and Larson [2018] Y. Wei, M. J. Solomon, and R. G. Larson, “A multimode structural kinetics constitutive equation for the transient rheology of thixotropic elasto-viscoplastic fluids,” J. Rheol. 62, 321–342 (2018).
  * Johnston [2006] D. C. Johnston, “Stretched exponential relaxation arising from a continuous sum of exponential decays,” Phys. Rev. B 74, 184430 (2006).
  * Berberan-Santos, Bodunov, and Valeur [2005] M. N. Berberan-Santos, E. N. Bodunov, and B. Valeur, “Mathematical functions for the analysis of luminescence decays with underlying distributions 1. Kohlrausch decay function (stretched exponential),” J. Chem. Phys. 315, 171–182 (2005).
  * Ashby [2011] M. F. Ashby,  _Materials Selection in Mechanical Design_ , 4th ed. (Elsevier, 2011).
  * Bird, Armstrong, and Hassager [1987] R. B. Bird, R. C. Armstrong, and O. Hassager,  _Dynamics of Polymeric Liquids_ , Vol. 1 (Wiley, 1987).
  * Piau [2007] J. M. Piau, “Carbopol gels: Elastoviscoplastic and slippery glasses made of individual swollen sponges Meso- and macroscopic properties, constitutive equations and scaling laws,” J. Non-Newtonian Fluid Mech. 144, 1–29 (2007).
  * Divoux, Grenard, and Manneville [2013] T. Divoux, V. Grenard, and S. Manneville, “Rheological hysteresis in soft glassy materials,” Phys. Rev. Lett. 110, 018304 (2013).
  * Blackwell  _et al._ [2015a] B. C. Blackwell, M. E. Deetjen, J. E. Gaudio, and R. H. Ewoldt, “Sticking and splashing in yield-stress fluid drop impacts on coated surfaces,” Phys. Fluids 27, 043101 (2015a).
  * Sen, Morales, and Ewoldt [2020] S. Sen, A. G. Morales, and R. H. Ewoldt, “Viscoplastic drop impact on thin films,” J. Fluid Mech. 891, A27 (2020).
  * Armstrong  _et al._ [2016] M. J. Armstrong, A. N. Beris, S. A. Rogers, and S. A. Wagner, “Dynamic shear rheology of a thixotropic suspension: Comparison of an improved structure-based model with large amplitude oscillatory shear experiments,” J. Rheol. 60, 433–450 (2016).
  * Ewoldt, Johnston, and Caretta [2015] R. H. Ewoldt, M. T. Johnston, and L. M. Caretta, “Experimental challenges of shear rheology: how to avoid bad data,” in  _Complex Fluids in Biological Systems_ , edited by S. Spagnolie (Springer Biological Engineering Series, 2015) pp. 207–241.
  * Singh, Soulages, and Ewoldt [2019] P. K. Singh, J. M. Soulages, and R. H. Ewoldt, “On fitting data for parameter estimates: residual weighting and data representation,” Rheol. Acta 58, 341–359 (2019).
  * Tikhonov and Arsenin [1979] A. N. Tikhonov and V. Y. Arsenin, “Solutions of ill-posed problems,” SIAM Rev. 21, 266–267 (1979).
  * Kontogiorgos [2010a] V. Kontogiorgos, “Calculation of relaxation spectra from mechanical spectra in MATLAB,” Polym. Test. 29, 1021–1025 (2010a).
  * Kontogiorgos [2010b] V. Kontogiorgos, “Calculation of relaxation spectra from stress relaxation measurements,” in  _Biopolymers_ , edited by M. Elnashar (IntechOpen, Rijeka, 2010) Chap. 25.
  * [43] “MATLAB regtools by Per Christian Hansen,” https://www.mathworks.com/matlabcentral/fileexchange/52-regtools.
  * Goodeve and Whitfield [1938] C. Goodeve and G. Whitfield, “The measurement of thixotropy in absolute units,” Trans. Faraday Soc. 34, 511–520 (1938).
  * Moore [1959] F. Moore, “The rheology of ceramic slips and bodies,” Trans. Br. Ceram. Soc. 58, 470–494 (1959).
  * Sen, Morales, and Ewoldt [2021] S. Sen, A. G. Morales, and R. H. Ewoldt, “Thixotropy in viscoplastic drop impact on thin films,” Phys. Rev. Fluids 6, 043301 (2021).
  * Jamali, Armstrong, and McKinley [2019] S. Jamali, R. C. Armstrong, and G. H. McKinley, “Multiscale nature of thixotropy and rheological hysteresis in attractive colloidal suspensions under shear,” Phys. Rev. Lett. 123, 248003 (2019).
  * Anscombe [1973] F. J. Anscombe, “Graphs in statistical analysis,” Am. Stat. 27, 17–21 (1973).
  * Freund and Ewoldt [2015] J. B. Freund and R. H. Ewoldt, “Quantitative rheological model selection: good fits versus credible models using Bayesian inference,” J. Rheol. 59, 667–701 (2015).
  * Mours and Winter [1994] M. Mours and H. H. Winter, “Time-resolved rheometry,” Rheol. Acta 33, 385–397 (1994).
  * Laukkanen, Winter, and Seppälä [2018] O.-V. Laukkanen, H. H. Winter, and J. Seppälä, “Characterization of physical aging by time-resolved rheometry: fundamentals and application to bituminous binders,” Rheol. Acta 57, 745–756 (2018).
  * Geri  _et al._ [2018] M. Geri, B. Keshavarz, T. Divoux, C. Clasen, D. J. Curtis, and G. H. McKinley, “Time-resolved mechanical spectroscopy of soft materials via optimally windowed chirps,” Phys. Rev. X 8, 041042 (2018).
  * Ewoldt and Saengow [2022] R. H. Ewoldt and C. Saengow, “Designing complex fluids,” Annu. Rev. Fluid Mech. 54, 413–441 (2022).
  * Corker  _et al._ [2019] A. Corker, H. C.-H. Ng, R. J. Poole, and E. García-Tuñón, “3D printing with 2D colloids: designing rheology protocols to predict ‘printability’ of soft-materials,” Soft Matter 15, 1444–1456 (2019).
  * Blackwell  _et al._ [2015b] B. C. Blackwell, A. E. Nadhan, A. Wu, and R. H. Ewoldt, “Fire suppression with yield-stress fluids,” Gallery of Fluid Motion, APS Division of Fluid Dynamics 68th Annual Meeting (2015b).
  * Sen and Ewoldt [2018] S. Sen and R. H. Ewoldt, “Stick or Splash? Thixotropy Decides!” Gallery of Fluid Motion, APS Division of Fluid Dynamics 71st Annual Meeting (2018).