Anticipating gelation and vitrification with medium amplitude parallel superposition (MAPS) rheology and artificial neural networks

Anticipating qualitative changes in the rheological response of complex fluids (e.g., a gelation or vitrification transition) is an important capability for processing operations that utilize such materials in real-world environments. One class of complex fluids that exhibits distinct rheological states are soft glassy materials such as colloidal gels and clay dispersions, which can be well characterized by the soft glassy rheology (SGR) model. We first solve the model equations for the time-dependent, weakly nonlinear response of the SGR model. With this analytical solution, we show that the weak nonlinearities measured via medium amplitude parallel superposition (MAPS) rheology can be used to anticipate the rheological aging transitions in the linear response of soft glassy materials. This is a rheological version of a technique called structural health monitoring used widely in civil and aerospace engineering. We design and train artificial neural networks (ANNs) that are capable of quickly inferring the parameters of the SGR model from the results of sequential MAPS experiments. The combination of these data-rich experiments and machine learning tools to provide a surrogate for computationally expensive viscoelastic constitutive equations allows for rapid experimental characterization of the rheological state of soft glassy materials. We apply this technique to an aging dispersion of Laponite® clay particles approaching the gel point and demonstrate that a trained ANN can provide real-time detection of transitions in the nonlinear response well in advance of incipient changes in the linear viscoelastic response of the system.


Introduction
Identifying incipient changes in the state of engineering systems is a critical task.This process, which typically comprises recognizing states of the system characteristic of sudden qualitative change in performance, is sometimes referred to as "fault detection" or "health monitoring."A prominent example of engineering systems in which incipient changes of state are identified are the mechanical structures that constitute civil infrastructure, such as bridges and buildings (Doebling et al. 1998).Structures fatigue over time, a physical change, and can suddenly lose rigidity resulting in a total mechanical failure (Bathias and Pineau 2010).
This study is dedicated to James W. Swan who passed away during the preparation of this study.
B Gareth H. McKinley gareth@mit.eduExtended author information available on the last page of the article Structural health monitoring can be accomplished in different ways, but in the most general sense, these techniques amount to methods for nonlinear system identification (NLSI) (Ewins 2009;Doyle et al. 2002).That is, a non-destructive probe of a system's structural mechanics is used to construct a nonlinear mathematical model for its structural response, and the nonlinear components of the model are subsequently monitored for changes that are empirically correlated with different modes of failure (Dwivedi et al. 2018).
One of the most common modes of probing structures is familiar to rheologists: the propagation of acoustic waves through the structure and the measurement of nonlinear interactions among the waves (Solodov et al. 2004).These techniques measure the weak nonlinear response through the use of special intermodulation tones whose nonlinear interactions can be easily measured (Cheng et al. 2017;Boyd et al. 1983;Chua and Liao 1989).Although there is limited rigorous proof that failure modes can be identified first in the nonlinear response, this kind of structural health monitor-ing sees wide applicability in the civil (Brownjohn 2007) and aerospace (Giurgiutiu 2015) industries.Experimental studies have demonstrated that harmonic analysis also provides a useful characterization of fatigue dynamics in soft materials; for instance, one study applied stress-controlled oscillatory shear tests for time-resolved monitoring of yielding and fluidization in a colloidal gel (Perge et al. 2014).More recent work has suggested that harmonic analysis may also provide sensitive indications of fatigue in brittle polymers and visco-hyper-elastic elastomers (Hirschberg et al. 2018(Hirschberg et al. , 2019(Hirschberg et al. , 2021)).
For complex fluids, methods of monitoring certain aspects of the state of a fluid that are critical for engineering applications, such as its linear rheology, both online (Covas et al. 2004;Konigsberg et al. 2013;Luger and Miethlinger 2019) and offline (Mours and Winter 1994) are robust.However, methods of measurement, either on-or offline, for anticipating dramatic changes in the linear rheology do not exist.Such transitions might come about through physical or chemical changes in the material and can be categorized in terms of many well-known phenomena including gelation and vitrification.When certain features of the linear rheology, or features of the nonlinear rheology that can be correlated with the linear response (such as the steady flow curve (Cox and Merz 1958)), are crucial for engineering applications, anticipating these transitions is essential to the design and maintenance of products and processes.In this work, we apply the rheological version of the intermodulation technique used in structural health monitoring, termed medium amplitude parallel superposition (MAPS) rheology (Lennon et al. 2020b), to study incipient changes in a fluid's rheological state (e.g., if it is in a "sol" or "gel" state) as characterized by its linear viscoelastic response.
Although empirical correlations can be used for this task, we seek a first principles foundation for such a measurement.To do so, we investigate a well understood model that is capable of describing the linear and nonlinear rheology of a broad range evolving and aging viscoelastic materialsthe soft glassy rheology (SGR) model (Sollich et al. 1997;Sollich 1998).This model is appropriate for describing materials whose micro-and meso-scale elements must yield in order to induce flow, such as aqueous dispersions of clay platelets, hard sphere colloidal glasses, and dispersions of polymer microgels (Bhattacharyya et al. 2023).It has been applied to describe dense suspensions (Purnomo et al. 2008), foams (Höhler and Cohen-Addad 2005), colloidal gels (Yin and Solomon 2008), emulsions (Mason et al. 1995), and dispersions of soft particles (Bonn et al. 2002;Suman and Joshi 2020).Qualitatively, the linear viscoelastic response of this nonlinear constitutive model exhibits three distinct regimes: single mode Maxwell behavior typical of simple viscoelastic fluids (Tschoegl 1989), power-law or fractional Maxwelllike behavior characteristic of gelling materials (Palade et al. 1996), and age-dependent fractional Kelvin-Voigt behavior typical of glassy materials (Papoulia et al. 2010).Which class of behavior is exhibited by the model is controlled entirely by a single parameter called the "noise temperature."The progressive transition or mutation in time of a given material from one regime to another will change the dynamical characteristics of any flow system employing such a fluid.Thus, a method of monitoring the "rheological health" of the fluid for incipient changes in its material response is critical to any engineering problem involving the rheomechanical characteristics of these fluids.
The paper is organized as follows.In Sect.2, we review the basic formulation for MAPS rheology to provide a foundation for intermodulation techniques in rheology.In Sect.3, we solve for the weak nonlinearities in both the scalar form and a tensorial form of the SGR model analytically, and qualitatively illustrate how these nonlinearities may be used to anticipate changes in the linear rheology of an aging or timeevolving fluid that can be described by the SGR model.In Sect.4, we introduce a machine learning method for "inverting" the SGR model, which enables real-time identification of incipient changes in the state of the fluid.Finally, in Sect.5, we demonstrate an implementation of this methodology on a gelling aqueous dispersion of Laponite ® clay platelets-a material with practical engineering uses that span oilfield and cosmetic applications (Cummins 2007;Xiong et al. 2019).

Background on MAPS rheology
A version of the intermodulation technique commonly applied to structural health monitoring of material systems has recently been developed for measurements of the weakly nonlinear mechanical response of complex fluids, and we refer to this as MAPS rheology (Lennon et al. 2020b).In a MAPS experiment, multiple sinusoidal shear waves are superimposed in a rheometer, and the time-resolved shear stress response is measured (Lennon et al. 2020a).For small strains, the material is homogeneously deformed, and the state of the fluid can be described in terms of the timedependent shear stress and shear strain alone.The analysis in this scenario is best performed in terms of the Fourier transformations of the shear stress and shear strain: In the limit of relatively small strain amplitudes, the Fourier transformation of shear stress is well approximated by a truncated polynomial expansion in the Fourier transformation of the strain: This is the frequency-space representation of a Volterra series expansion for the shear stress in terms of the strain (Lennon et al. 2020b).The quantity G * (ω) is the familiar complex modulus measured in ordinary oscillatory rheometry at very small strains.Due to time-reversability, the leading order nonlinearity in the shear stress during homogeneous deformation arises at third order in the strain, and the response function, G * 3 (ω 1 , ω 2 , ω 3 ), is called the third-order complex modulus with SI units of Pa.This material function encodes the weakly nonlinear response of the shear stress to a homogeneous deformation with arbitrary time dependence.In previous work (Lennon et al. 2020b), we have analyzed this function for several constitutive models and showed that the frequency dependence of G * 3 (ω 1 , ω 2 , ω 3 ) is very sensitive to the details of the general functional relationship between stress and deformation in each model.
Many soft materials are fragile and can only sustain small deformations before the imposed flow drives significant evolution of the microstructure (Butler 1999).In such a case, measurements of weak nonlinearities under strain control can be difficult.However, the weakly nonlinear strain response under an imposed shear stress history can be measured instead.The Volterra series expansion of the strain in terms of the shear stress in frequency space is as follows: This stress-controlled mode of experimentation provides measurements of the linear compliance: for the linear response, along with another Volterra kernel, called the third-order complex compliance, J * 3 (ω 1 , ω 2 , ω 3 ).Because the functional relationship between shear strain and stress is invertible, this material function can be related to the third-order modulus directly through analysis of Eqs. 1 and 2: Thus, the third-order complex compliance inherits the same model sensitivity present in the complex modulus.
In a typical MAPS experiment with characteristic strain amplitude γ 0 1, we use a superposition of three particular sinusoids (Lennon et al. 2020a): The particular value of the amplitude γ 0 should be selected based on an analysis of experimental variance (i.e., noise) and the relative bias generated by higher order effects in the material under study (Lennon et al. 2020a).The set of integers {n 1 , n 2 , n 3 } should be chosen carefully, so that all sums and differences of triplets selected from the set are unique.In this case, the third-order complex modulus can be measured directly as distinct harmonics of ω 0 in the shear stress.The data returned by a single MAPS experiment with three distinct input tones consists of both the linear response at frequencies: n 1 ω 0 , n 2 ω 0 , n 3 ω 0 and the nonlinear response function at the frequency triplets listed in Table 1.Sweeps of the fundamental frequency, ω 0 , and different sets of integer triplets, {n 1 , n 2 , n 3 }, can allow one to probe the third-order MAPS response function in its entirety.Just as a three-tone strain waveform can be used to measure G * (ω) and G * 3 (ω 1 , ω 2 , ω 3 ), a three-tone shear stress waveform can be used to measure J * (ω) and J * 3 (ω 1 , ω 2 , ω 3 ).
Table 1 The harmonics of the base tone ω 0 at which the Fourier transformation of the stress shows a spectral contribution to the three-dimensional function G * 3 (ω 1 , ω 2 , ω 3 ) arising from a three-tone, strain-controlled MAPS experiment

Harmonic
Measured stress response / iπγ 3 0 /4 Previous studies have demonstrated that three-tone MAPS experiments can provide data that is more sensitive to distinct aspects of a material's nonlinear rheology than other spectral and transient rheometric methods (Lennon et al. 2021b) and that this sensitivity leads to more precise estimates of constitutive model parameters when fitting the model predictions to real experimental data (Lennon et al. 2021a).This sensitivity is critical for health monitoring, because the monitoring system strives to identify small changes in the state of the material.Three-tone MAPS experiments also possess an advantage in data density over other experimental techniques, as evidenced by Table 1.The acquisition time of such experiments is set by the fundamental frequency ω 0 -as is also the case for single-tone small and medium amplitude oscillatory shear experiments (SAOS and MAOS, respectively)-but the number of independent measurements is set by the number of input tones, with the three-tone experiment resulting in nearly ten times more data than its single-tone counterpart.It has recently been shown that the data throughput of multi-tone MAPS experiments can be further increased by implementing signals with more than three tones, though these experiments require a more complicated experimental design that will not be considered here (Lennon 2023).As we will demonstrate in Sect.4, this data richness enables accurate health monitoring from a single MAPS experiment, which can facilitate the development of online rheology health monitoring tools.

Weakly nonlinear response of the SGR model
The scalar soft glassy rheology (SGR) model (Sollich et al. 1997;Sollich 1998) describes the flow of an amorphous material due to local yielding events having a distribution of yield energies E, denoted ρ(E).For a material subjected to an oscillating strain history, the relationship in periodic steady state (the alternance state) between the shear stress, σ (t), and the strain, γ (t), in the SGR model can be written as follows: with auxiliary equations defining the rate at which yielding events occur, (t), and an effective time interval, Z (t, t ): and The function G(z) is a relaxation modulus, defined as follows: with τ = e E/x representing a dimensionless relaxation time, which depends on the yield energy E, and the average in angle brackets being taken over the equilibrium distribution of yield energies: eq ρ(E)e E/x .The parameter x is a socalled "noise temperature" setting the scale for the depth of the local energy minima trapping the yielding elements.In this model, time is made dimensionless on the inverse of the attempt frequency for yielding, f −1 , and stress is made dimensionless on the characteristic stiffness of the yielding elements, k.We note that throughout this section, we use the symbols σ and G to denote the dimensionless shear stress and modulus, respectively, and the symbols t and ω to denote dimensionless time and frequency, respectively.Although Eqs. 5-8 are an exact solution to the SGR model with periodic straining, computing response functions with these equations analytically may be impossible in general.Computing these functions numerically is also quite laborious because of the integral equation describing how the yielding rate changes in time under flow (Eq.6) (Radhakrishnan and Fielding 2016Fielding , 2018;;Park and Rogers 2018).However, asymptotic solutions to this set of equations can be constructed through straightforward methods (Fielding et al. 2000), and these are all that are needed to extract the weakly nonlinear response probed in a MAPS experiment.To O(γ 2 ), the time interval can be written as follows: Similarly, to O(γ 2 ), the relaxation modulus G(Z (t, t )) is as follows: As a final ingredient in this asymptotic analysis, we supply the ansatz that the yielding rate to leading order is (t) = eq (1 + ˜ 2 (t)) + O(γ 3 ), where ˜ 2 (t) ∼ O(γ 2 ).Substituting these expansions for G(Z (t, t )) and (t) into Eq.6 and matching terms up to O(γ 2 ) shows that which can be solved to give the perturbation to the equilibrium yield rate, ˜ 2 (t).With these same expansions, the stress to leading order in the strain can be written as follows: where the first line gives the linear response and the remainder is the leading-order nonlinearity.
To find a closed form for the weakly nonlinear, timedependent response of the SGR model, it is convenient to work with the stress and strain in Fourier space: σ (ω), γ (ω).Then, Eq. 12 can be written as follows: . (14) The Fourier transformation of Eq. 11 gives the perturbation to the yielding rate as follows: Finally, substituting for ˆ 2 (ω) in the transformed stress and integrating over the frequency variable ω yields the following: with the third order complex modulus defined as follows: The leading-order term in the Fourier transformation of the stress can be identified as the linear response, with the transfer function being the complex modulus of the material.The leading nonlinearities are identified with the remaining terms in Eq. 17.This nonlinear response is third order and controlled by the Volterra kernel, G * 3 (ω 1 , ω 2 , ω 3 ).As with other third-order Volterra kernels, the expression in Eq. 18 can be made symmetric without altering the structure of the response by permuting the variables ω 1 , ω 2 , ω 3 and averaging the permutations (Lennon et al. 2020b, a).Appendix 1 contains a non-symmetric formula for G * 3 (ω 1 , ω 2 , ω 3 ) written in terms of confluent hypergeometric functions.Figure 1 plots the real and imaginary parts of for noise temperatures of x = 4 and x = 5 and a few distinctive values from the set of integer triplets: {n 1 , n 2 , n 3 }.The third-order moduli exhibit particular power-law scaling in the low-frequency limit.The real and imaginary parts of the third-order moduli also exhibit sign changes with changes in frequency-a phenomenon that is known from MAOS studies of other nonlinear viscoelastic models (Ewoldt and Bharadwaj 2013).
In the SGR model, there are postulates about the distribution of yield energies necessary to produce a glass transition (Bouchaud 1992).It is generally assumed that ρ(E) ∼ exp(−E), so that the equilibrium distribution of yield energies, eq ρ(E)e E/x , remains integrable up to x = 1, Fig. 1 The third-order complex modulus: , of the SGR model with x = 4 and x = 5 as a function of the dimensionless fundamental frequency ω 0 for {n 1 , n 2 , n 3 } = {1, 4, 16} (black symbols), {5, 6, 9} (red symbols), and {1, 1, 1} (blue symbols).The left panel shows scale-free power-law behavior at low frequency, while the right shows the expected behavior at third order for a viscoelastic fluid with a compact spectrum of relaxation times where the glass transition occurs.This exponential distribution also gives rise to regimes in which the predicted material response can be described by a fractional viscoelastic constitutive equation (Jaishankar and McKinley 2013;Rathinaraj et al. 2021).For the present purposes, we take the distribution to be exp(−E), so that eq = (x − 1)/x and evaluate the averages in the expansion of the shear stress.As noted by Sollich (1998), the linear response produced by this distribution of relaxation times is as follows: which has the well-known (Bonn et al. 2002) asymptotic scaling in the low-frequency limit: The noise temperature in the SGR model thus controls a set of rheological transitions in the material.For x > 3, the linear viscoelastic response in both G and G is that of a simple Maxwell fluid.For 1 < x < 3, the linear response is of the fractional Maxwell type (Bonn et al. 2002;Papoulia et al. 2010).Below the glass transition, 0 < x < 1, the linear response takes on a fractional Kelvin-Voigt form, which is also age-dependent (Fielding et al. 2000).We do not explore this age-dependent regime here, but refer the reader to the work of Fielding et al. (2000).
The approach to the glass transition at x = 1 is often analyzed in terms of the phase angle δ(ω) of the complex modulus, defined by tan δ(ω) = G (ω)/G (ω).For large x In this last regime, the phase angle becomes independent of frequency (at low frequencies).This is often called the critical gel state (Winter and Chambon 1986), in which the relative viscoelastic characteristic (i.e., the phase angle given by the ratio of G to G ) is scale-free (i.e., independent of the observation or deformation timescale, in the limit of long time-scales or small frequencies).Observing this scale-free response in which tan δ(ω) is constant can be a sensitive probe of the gelation time in an aging soft glass (Keshavarz et al. 2021).
The third-order modulus can be evaluated in a similar manner so that Eq. 18 yields a complex expression involving hypergeometric functions.We omit the full expression here as it is not easily simplified or interpreted (see Eqs. A.1-A.4).However, in the limit of low dimensionless frequencies, 1, this expression has the asymptotic form: Perhaps unsurprisingly, in the regime just above the glass transition 1 < x < 2, the third-order modulus is scale-free just as the linear response.That is, the real and imaginary parts of the third-order modulus G * 3 (ω 1 , ω 2 , ω 3 ) exhibit the same power-law dependence, or logarithmic slope, at low frequency.The ratio of regime.This is in contrast to a third-order fluid with a Maxwell-like linear response, such as the SGR model with x > 5, for which this same ratio grows linearly with ω (Lennon et al. 2020b).One surprising characteristic of the 1, as a function of noise temperature x.There exists a distinct regime (shaded) in which the logarithmic slopes of G 3 and G 3 serve as antecedents of a transition in the character of the linear response third-order modulus in the SGR model is that it becomes scale-free further from the glass transition at x = 1 than is the case with the linear response.That is, scale-free behavior appears in the linear response when 1 < x < 2, while it appears in the third-order complex modulus when 1 < x < 4, which suggests that the weak nonlinearities in the SGR model are more sensitive to the incipient glass transition that is parameterized by a time-varying and slowly decreasing noise temperature x(t).
Being able to establish and monitor the rheological state of a soft material-whether it is presently fluid-, gel-, or glass-like in characteristic-is critical to use of that fluid in applications that require a particular rheological response.Figure 2 illustrates how one could potentially use the log-arithmic slopes of the third-order modulus, d ln G 3 /d ln ω, as advanced warning of a crossover to a scale-free linear response.If that scale-free behavior were to correspond to the process of gelation and x is decreasing along a reaction coordinate for the sol-gel transition, then the observation of scale-free response in the third-order modulus (e.g., in the range 3 < x < 4) could be used to anticipate the transition to the gelled regime well in advance of any qualitative change in the linear response (1 < x < 2).
As with the moduli, the linear and third-order compliance of the SGR model exhibit distinct power law behaviors at small frequencies much like the moduli: In Fig. 3, we plot the real and imaginary parts of for the SGR model with x = 4 and x = 5 obtained from Eqs. 18 and 3, with the same values of the set {n 1 , n 2 , n 3 } as in Fig. 1.The predicted asymptotic limits are confirmed.
From these values of the third-order compliance, the logarithmic slopes of the real and imaginary parts can be extracted and the response tested for scale-free behavior.When the nonlinear compliance becomes scale-free at low frequency as indicated in Fig. 4, the material is in the same pre-vitrified state as shown in Fig. 2 and is about to undergo a rheological transition.

Fig. 3
The third-order complex compliance, J * 3 (n 1 ω 0 , n 2 ω 0 , n 3 ω 0 ), of the SGR model with x = 4 and x = 5 as a function of the dimensionless fundamental frequency ω 0 for {n 1 , n 2 , n 3 } = {1, 4, 16} (black symbols), {5, 6, 9} (red symbols), and {1, 1, 1} (blue symbols).The left panel shows scale-free power-law behavior at low frequency, while the right shows the expected behavior at third order for a viscoelastic fluid with a compact spectrum of relaxation times Fig. 4 Logarithmic slopes of the linear and weakly nonlinear response functions of the SGR model in a stress-controlled experimentcompliances: 1, as a function of noise temperature x.There exists a distinct regime (shaded) in which the logarithmic slopes of J 3 and J 3 serve as antecedents of a transition in the character of the linear response One practical complication to this picture is that within these asymptotic power law regimes, the nonlinear response functions can exhibit changes in sign.Unlike the linear response, there is no physical constraint that forces the nonlinear response function to be positive over any given frequency range.This means that unlike examinations of processes like critical gelation, the power law slopes cannot necessarily be read off a log-log plot of the nonlinear response function data.In light of limited data sets, a more sophisticated analysis of the data is required and will be examined in the following section.
An additional consideration in this sort of analysis is whether the weak nonlinearities that are predicted in the shear stress of a tensorial SGR model differ significantly from the scalar model.In the tensorial SGR model, the stress tensor is written as follows (Cates and Sollich 2004): with the effective time increment where is the accumulated strain tensor between t and t , and R and Q are functions of the accumulated strain tensor.The tensorial SGR formulation also gives predictions of normal stresses and the response of soft glassy materials under other kinematic deformation histories.In their work, Cates and Sollich (2004) noted that the normal stresses of the tensorial SGR measured in a step strain experiment change their behavior qualitatively at a different noise temperature than the shear stress.This is a subtle historical precursor to the predictions for the MAPS response we have made here.Normal stress differences are notoriously difficult to measure accurately in experiments (especially in soft gels and polymer solutions), whereas intermodulation responses such as the ones measured in a MAPS experiment are straightforward to detect in simple shear with standard instrumentation and offer the additional benefit of spectral noise filtering in the Fourier domain.
The function R(E(t, t )) must be zero with no accumulated strain and can only depend on the scalar invariants of E(t, t ).Therefore, for simple shear and small strain amplitudes, an appropriate functional is as follows: The function Q controlling the stress generated by yielding elements must be traceless and respect the same time-reversal symmetries as the macroscopic stress.At small strain amplitudes, Q, should be proportional to just the accumulated strain in order for the tensorial SGR model to return a linear response.A natural small strain amplitude expansion for the shear component of Q is thus as follows: The coefficients c and d are free parameters and are material specific.Adding a tensorial character to the SGR model changes the MAPS response, G * 3 (ω 1 , ω 2 , ω 3 ) in two key ways by mixing the nonlinearities arising from R and Q: First, the MAPS response of the scalar SGR model is reweighted by the coefficient c, which reflects the fraction of energy stored in the yielding element that is used as work.
Necessarily, c > 0, and ought to be an O(1) quantity.Second, the leading order nonlinearities in the stress produced by a yielding element create an additional contribution to the MAPS response that is time-strain separable and proportional to the stress nonlinearity coefficient, d.If this term arises purely kinematically as in Larson's model of foam rheology (Larson 1997), then the coefficient takes the value d = −1/8.
If, however, there are other modes of nonlinearity in the stress generated by an element, d could potentially deviate from this value in both magnitude and sign.Although the MAPS response of the tensorial SGR model has an additional contribution when compared with the scalar SGR model, this time-strain separable part exhibits the same power law scaling in the low-frequency limit.Consequently, Figs. 2 and 4 apply to the tensorial SGR model as well, while Figs. 1 and 3 are just qualitative descriptors.

Machine learning to accelerate solution of the inverse problem
Recall that the scalar SGR model has three unknown parameters: the attempt frequency f on which frequencies are made dimensionless, the characteristic stiffness of yielding elements k, on which moduli are made dimensions, and the noise temperature x.These appear in both the linear response and the leading order nonlinearities.The tensorial SGR introduces another two parameters: c and d, that appear in the nonlinearities alone.For low oscillatory frequencies, at which the SGR model is expected to be valid, a fit to linear response data cannot accurately determine the three parameters: k, f , and x, independently.From Eq. 20, we see that in the scale-free regime (1 < x < 2), the sensitivity of the linear response in the limit that ω → 0 to the parameter x, given by ∂G * (ω)/∂ x, diverges (note that for the remainder of this work, we return to the convention that ω represents dimensional frequency and G * 3 represents the dimensional third-order complex modulus).The diverging sensitivity to x results in an ill-conditioned problem of inverting the model to find x, k, and f simultaneously.In the region where x > 3, one can determine two unique combinations of k, f , and x, but cannot find all three independently as the model adopts a single mode Maxwell form, which possesses just two parameters: a characteristic modulus and a time scale.Parameter estimation for the SGR model with the linear response alone is thus ill-conditioned.
If we instead consider the weakly nonlinear case, we find that the MAPS response functions are parameterized by k, f , and x as well as the coefficients c and d.Although these functions also possess a scale-free region, that region encompasses the one covering the linear response.Therefore, a combination of the linear response and nonlinear response data can be used to determine each of the five parameters independently.Unlike the case for the linear response alone, parameter estimation when combining linear response and MAPS data is well conditioned across a broad range of parameter values.
One approach to this parameter estimation problem is regression.An estimate of the parameters is formed from minimization of an objective function, which is given by the difference between measured and modeled values of the linear and nonlinear response functions at several different frequencies.A global minimum of this difference could be identified by varying the model parameters systematically through the use of schemes like simulated annealing (Kirkpatrick et al. 1983) and parallel tempering (Earl and Deem 2005).One challenge with this approach is that the model for the nonlinear response functions is computationally expensive to evaluate.The process of searching for a global minimum from a given set of experimental data could become very time consuming.Repeating that process over and over is impractical, especially if one wants to learn this parameterization across a wide variety of samples or in an online fashion during an industrial process.
An alternative approach is to use an artificial neural network (ANN), trained on synthetic data, to learn the inverse functional relationship between the linear response and MAPS data and the model parameters k, f x, c, and d.Artificial neural networks are powerful tools for supervised regression problems, in which the objective is to learn the relationship between input variables and one or more continuous output variables from a data set containing examples with known inputs and known outputs (Ferguson 2018).Neural networks have demonstrated success in learning inverse solutions to other constitutive models when trained on synthetic data (Mahmoudabadbozchelou and Jamali 2021;Mahmoudabadbozchelou et al. 2022).Here, we use a multilayer perceptron to learn this inverse relationship within the SGR model for a data set that combines the linear response and weakly nonlinear response measured in a single three tone MAPS experiment.Schematically, this inverse function gives the relationship (k, f , x, c, d) = f (y), with depending on whether strain-or stress-controlled MAPS experiments are performed.The ellipses describe the 17 other weakly nonlinear response functions measured in a MAPS experiment.The 19 strain-controlled functions are described in Table 1, and the corresponding stress-controlled response functions (third order compliance) are measured at the same sets of frequencies.Artificial neural networks of sufficient size are capable of learning arbitrary functional relationships between their inputs and outputs (Hornik et al. 1989).They do this by composing a set of "layers" in series, with a single layer consisting of a nonlinear (and typically non-parametric) activation function, a(.), and a set of linear equations parameterized by weights W and biases b.The state stored by layer n, denoted y n , is therefore related to the state stored by layer n − 1 by the following: For an ANN with N layers, n = 0 represents the model inputs (the "input layer"), and n = N represents the model outputs (the "output layer"), with 0 < n < N typically called the "hidden layers" of the model.The size of y 0 and y N is set by the shapes of the model inputs and outputs, while the size of the hidden states y n may be selected arbitrarily or by a set of heuristic rules.Denoting the size of state y n as l n , the weight matrix for layer n is l n × l n−1 , and the bias vector is l n × 1.The trainable parameters of the ANN, which are updated to optimize an objective related to the true and predicted outputs, are the combined set of weights and biases from each layer.Figure 5 presents a schematic depiction of the structure of an ANN.Note that like the size of the hidden layers, the total number of layers in the network N is typically chosen arbitrarily or heuristically.
Although we may apply the ANN architecture described above to an arbitrary set of input data, including data from MAPS experiments with arbitrary {n 1 , n 2 , n 3 }, in the present Fig. 5 Schematic depiction of the ANNs used to learn the relationship between data collected in a single three-tone MAPS experiment and the SGR model parameters.The network is composed of an input layer, one or more hidden layers, and an output layer.With the exception of the output, the states stored in successive layers are related by the composition of a (non-parametric) nonlinear activation function and a linear model parameterized by a set of weights and biases work, we apply ANNs to inputs from stress-controlled threetone MAPS experiments with {n 1 , n 2 , n 3 } = {5, 6, 9}.As we will discuss in Sect.5, this set of input tones and the mode of stress control is optimal for a stress-controlled rheometer with inertial limitations at high frequency.However, the same limitations may not apply to other experimental systems; thus, we explore the application of ANNs to other experimental protocols in Appendix 2.
To assess whether the combination of three-tone MAPS experimental data and machine learning can overcome the limitations of health monitoring from linear rheology alone, we construct and train two ANNs: one whose inputs are only the three values of the linear compliance-J * (n 1 ω 0 ), J * (n 2 ω 0 ), and J * (n 3 ω 0 )-obtained from a three-tone experiment and one whose inputs are the entire data set (three linear response measurements and 19 MAPS measurements) obtained from this experiment.Because these response functions are complex and vary over orders of magnitude, the ANN inputs are constructed by concatenating the logarithm and sign of both the real and imaginary components of the raw input data (Eq.30) into a single 88-element vector.In both cases, the network outputs are the five parameters of the tensorial SGR model.The full architectures of these networks are presented in Table 2.These networks are trained using a set of 9000 synthetic experiments generated from the analytical predictions of the tensorial SGR model, with parameters selected uniformly at random from the following intervals: log 10 k ∈ [−1, 5], log 10 f ∈ [1, 5], x ∈ [1, 6], c ∈ [0, 10], and d ∈ [−1, 0].The networks were constructed using the Keras interface for Tensorflow and trained to minimize the mean squared error between the predicted and true SGR parameters.This criterion was optimized using the root mean squared propagation (RMSprop) optimizer.
The trained ANNs were next asked to predict the SGR parameters for a test set of 1000 synthetic MAPS experiments.Figure 6 presents parity plots of the predicted versus true parameters for this test set.The results are in agreement with many of the qualitative statements made throughout this work.Specifically, the ANN trained on linear response data The "Full" ANN is trained on the entire MAPS data set, while the "LR" (linear rheology) ANN is trained only on the linear data from a single MAPS experiment Conversely, the ANN trained on both linear and MAPS data (blue circles) provides more precise predictions of k and f , as the inclusion of MAPS data improves the conditioning of the problem.Moreover, this network is able to accurately track the noise temperature x well beyond x = 3, with parity observed to at least x = 5 as predicted by the qualitative scaling analysis summarized in Figs. 2 and 4. The parameters c and d are still difficult to ascertain, most likely due to the observation that the time-strain separable component of the response gives rise to the same asymptotic behavior as the non-separable component.Therefore, while the inclusion of the tensorial SGR formulation ultimately does not affect our ability to monitor the state of a soft glassy fluid using MAPS, the observation that MAPS data is insufficient to fully parameterize this model is interesting, as the limitations of MAPS rheology are not yet fully understood.However, the conclusion that x may be determined by this network for values of x > 3 is most significant, as x represents the reaction coordinate for gelation or vitrification that is the target of health monitoring.Evidently, it is possible to monitor this coordinate using a single three-tone MAPS experiment well before changes are evident in the linear rheology of the SGR model.
That these accurate predictions of the noise temperature in the SGR model may be obtained from a single MAPS experiment is noteworthy.To obtain the same amount of data from other well-known medium amplitude probes, such as MAOS, would require nearly ten times the number of experiments, because a single three-tone MAPS experiment measures the complex-valued third-order response function at a set of 19 intermodulation and third harmonic frequencies, while a single-tone MAOS experiment excites nonlinearities only on the first and third harmonic (i.e. two measurements of the third order response function).For online health monitoring applications, this amount of experimentation may not be feasible within process time constraints.Such applications are facilitated by the high data density of MAPS experiments, which substantially reduce the time needed to acquire experimental data sufficient for health monitoring.However, the application of machine learning techniques is equally important to these rapid predictions of material health that are necessary for online applications.As previously discussed, the parameters of the SGR model could have been recovered from MAPS data by other means, such as global optimization.In a five-dimensional parameter space, this would likely require thousands of evaluations of the model predictions, which are computationally burdensome.In fact, we find that each evaluation call requires O(1 s) to compute on a single computer (2019 MacBook Pro, 2.4 GHz 8-Core Intel i9 processor, 64 GB RAM); thus, direct optimization could take minutes or hours to converge, which precludes real-time monitoring of the fluid's evolving rheological state.In the machine learning approach, these expensive evaluations are performed only during training data generation, in an offline setting, and only evaluation of the trained model is required in the online setting.Moreover, the trained model needs only be evaluated once to estimate the noise temperature, and inference takes only O(1 ms) on a single computer.Thus, with machine learning, the monitoring timescale is limited only by the data acquisition time.

Experimental demonstration on a Laponite ® dispersion
The previous section demonstrates that, within the framework of the SGR model, MAPS rheology and machine learning may be combined to provide advanced monitoring of the rheological state of a material undergoing a glass transition.Although the SGR model qualitatively describes many real fluids, such fluids may not quantitatively agree with SGR predictions and may be influenced by other effects not captured by the SGR model.Moreover, true experimental MAPS data may not be as precisely determined as the synthetic data used in the previous demonstration.Whether the methodology outlined above extends to physical systems requires independent validation on real experimental data.
In this section, we apply the rheological health monitoring methodology to an aqueous dispersion of Laponite ® clay platelets.Hectorite clays such as Laponite ® and bentonite are employed in a variety of industrial and commercial applications, such as the base component of drilling muds utilized in oilfields, and as cosmetic products (Cummins 2007;Xiong et al. 2019).Numerous studies have investigated the linear rheology of these dispersions as they undergo gelation, showing a transition that is well described by the scaling features of the SGR model outlined in Sect. 3 (Bonn et al. 2002;Suman and Joshi 2020).However, researchers have not investigated whether these transitions in the linear rheology may be anticipated by transitions in the nonlinear rheology of the dispersions.

Materials
The Laponite ® clay powder (Laponite ® RD supplied by BYK) was dissolved in deionized water at 3 wt.%under ambient conditions and mixing at 1000 rotations per minute.The system was mixed for 10 min, at which point a homogeneous, slightly opaque dispersion was produced.This dispersion was immediately transferred to the rheometer for measurements.Fresh samples were prepared for each timeresolved set of gelation measurements.

Experimental methods
All experiments were performed in a DHR-3 Discovery Hybrid Rheometer from TA Instruments using a smooth concentric cylinder geometry (cup diameter of 30 mm and a gap of 1 mm) maintained at 25 • C, with data collected and analyzed using TRIOS v5.2.2.Frequency sweeps were performed in stress control with an amplitude of σ 0 = 0.01 Pa.Three-tone MAPS experiments were performed in stress control using the arbitrary wave feature, with amplitudes of σ 0 = 0.1 Pa and σ 0 = 0.05 Pa, and a fundamental frequency of ω 0 = 0.2 rad/s, with a sampling rate of 30 points/s that is sufficiently above the Nyquist rate for the anticipated thirdorder response.As previously noted, this stress-controlled mode of experimentation is particularly suitable for materials that develop a weak yield stress during aging, such as this clay dispersion, where a fixed strain amplitude may eventually induce microstructural yielding.The MAPS response functions were obtained from the raw time-series data using the data analysis routines outlined in Lennon et al. (2020a).

Results
We first perform a series of small amplitude oscillatory frequency sweeps to monitor the linear response of the Laponite ® dispersion as it approaches the gel transition.We sweep frequencies from 10 to 0.5 rad/s, with five frequencies per decade.At each frequency, the linear complex modulus is computed from two oscillation cycles.A single frequency sweep takes slightly less than 5 min to complete; thus, we run frequency sweeps every 5 min to monitor the material.We find that the material enters the scale-free regime at approximately 45 min after the cessation of mixing; thus, we conduct frequency sweeps until 50 min after mixing to monitor the linear rheology up to the point where noticeable changes are observed.Figure 7 depicts the measured loss modulus (G (ω)) for each frequency sweep.The storage modulus (G (ω)) for many experiments is below the noise floor, resulting in negative values.Thus, we exclude the storage modulus from the present analysis.
We similarly conduct three-tone MAPS experiments with {n 1 , n 2 , n 3 } = {5, 6, 9} for a freshly prepared sample.The fundamental frequency of ω 0 = 0.2 rad/s was selected to remain in the low-frequency regime, while not requiring unduly long experiments.Data was collected for eight periods with respect to this fundamental frequency, resulting in experiments just under 5 min in duration.These MAPS tests were run every 5 min, with a stress amplitude of σ 0 = 0.1 Pa.This amplitude was selected via a separate stress-controlled Fig. 7 The loss modulus G (ω) measured for the Laponite ® dispersion using small amplitude oscillatory shear.Experiments are conducted every 5 min beginning at 5 min after mixing, and continuing until 50 min after mixing.The lighter points with higher transparency represent data taken at earlier times, with data points connected by dashed lines to aid in readability amplitude sweep at ω = 0.2 rad/s, in order to mitigate the effects of variance due to noise and bias due to the growth of higher-order effects.The tests were repeated with a fresh sample for an amplitude of σ 0 = 0.05 Pa.For each MAPS experiment, we truncate the first four cycles of the time-series data to eliminate transient features in the response, and the resulting time series are analyzed to obtain the third-order complex compliance at the coordinates listed in Table 1.
We do not visualize the response functions computed using the MAPS experiments here, as previous work has noted the difficulty in visualizing such high dimensional data, and our previous visualization schemes rely on a sweep of the fundamental frequency ω 0 .However, the raw and processed MAPS data are included in an online repository associated with this work (https://github.com/krlennon/rheologyhealth-monitoring),and an example waveform and frequency spectrum of the MAPS stress-controlled protocol and shear strain response are shown in Fig. 8.We see from this figure that the input tones, first harmonic response, and multiple intermodulation effects are visible above the floor of the data.However, many features of the frequency spectra in particular make direct interpretation of the high frequency content challenging, such as peaks and minima not expected from the input or output waveforms, and the nonzero floor of the signals (which can be affected by the signal duration or the windowing applied during calculation of the discrete Fourier transform (Rathinaraj and McKinley 2023)).We also emphasize that the methods used to obtain the MAPS response Fig. 8 An example of the three-tone MAPS protocol used to obtain weakly nonlinear data characterizing the Laponite ® dispersion.The data depicted here is for the MAPS experiment run at t = 15 min after mixing, with a stress amplitude of σ 0 = 0.1 Pa. a The input three-tone stress waveform (with an input tone set {5, 6, 9} and ω 0 = 0.2 rad/s) and b the resulting absolute shear strain response (dimensionless), windowed to exclude the transient response associated with starting the deforma-tion protocol.c The real and d imaginary components of the discrete Fourier transform of both the stress (red) and strain (blue) waveforms.Solid lines depict the absolute value of the frequency spectrum.Peaks associated with the input tones and first harmonic response are marked with squares, and peaks associated with third-order intermodulation are marked with circles, with filled and unfilled symbols representing positive and negative values, respectively function uses the observed stress signal, rather than simply assuming that the commanded input is achieved exactly, and therefore, any deviations from the expected input waveform are accounted for during data processing.Figure 8, along with a figure depicting the corresponding MAPS test with σ 0 = 0.05 Pa, is also included in the GitHub repository associated with this work.
The set of input tones used for the MAPS experiments was {n 1 , n 2 , n 3 } = {5, 6, 9}.These tones result in distinct sums and differences for all triplet sets-a requirement for maximizing the data density of the experiment.This specific tone set was chosen to ensure that the highest frequency component of the weakly nonlinear response, occurring at three times the highest input frequency (3n 3 ω 0 = 3 × 9 × 0.2 rad/s = 5.4 rad/s), remained below the high-frequency inertial limit of the controlled stress rheometer fixture, which was observed to be 10 rad/s.Thus, this triplet set and the selected fundamental frequency ensure that experiments are rapid enough to enable real-time monitoring, while still operating within the capabilities of the instrumentation.
During data acquisition, we observed that the amplitudes required to distinguish the nonlinear features of the in-phase component of the strain response were large enough to disrupt the gelation process and fluidize the dispersion.This difficulty likely stems from the relatively weak fluid response, combined with spectral leakage from the linear response.Thus, only the out-of-phase MAPS response (corresponding to the imaginary component of J * 3 ) was retained in subsequent analysis.Similarly, because the measured storage modulus (and compliance) was below the noise floor, resulting in unphysical negative values, for many experiments, only the loss compliance J was retained.We note that the storage modulus does eventually become significant, and in the glassy state eventually dominates the linear response.This regime has been thoroughly characterized by other work (Rathinaraj et al. 2023); however, it is outside of the scope of the present study.
Following the methodology outlined in Sect.4, we trained ANNs to predict the parameters of the SGR model from both linear response and MAPS data.For the ANN trained on only the linear response, the network was trained using 10,000 synthetic frequency sweeps, with inputs consisting of the logarithm of the loss modulus at the eight measured frequencies.This ANN contained 3 hidden layers with 16 neurons each.For the ANN trained on MAPS data, the network was trained using 10,000 synthetic MAPS experiments, with inputs consisting of the logarithm of the absolute value of the linear response and MAPS response measured on each of the 22 output channels (three linear and 19 intermodulation).This ANN contained only a single hidden layer with 8 neurons.These network sizes were selected by conducted a brief hyperparameter optimization survey to determine the minimal network size to sufficiently resolve the model parameters from synthetic data (similar to the L-curve method).The networks were trained using the same optimizer and criterion described in Sect. 4.
In the previous sections, we have noted that models trained to invert the SGR model from linear response data provide no (or very limited) sensitivity to x for x > 3. The same lack of sensitivity occurs for MAPS-trained models when x > 5. Predictions in these regimes should be taken as highly uncertain, but the machine learning approach described in Sect. 4 provides no quantification of this uncertainty.Thus, we train additional neural networks to parameterize a distribution over x for each input data set.Because we expect networks to be biased towards lower values of x, we choose a Gumbel distribution (Gumbel 1941), which has a natural skew to account for a longer tail at higher x.The probability density function for this distribution is given by the following: We take the output value of x from the former neural networks as the mode μ of this distribution and the outputs of the second set of neural networks as the scale parameter β, trained to minimize the negative log-likelihood of the Gumbel distribution.Figure 9a depicts the SGR noise temperature inferred by the ANNs for SAOS and MAPS experiments as a function of the elapsed time between initial mixing of the dispersion and the beginning of the experiment, with symbols representing the mode μ of p(x; μ, β) and bars spanning the lower and upper quartiles: 0.25 < p(x; μ, β) < 0.75.At 45 and 50 min, when the dispersion transitions to a critical gel, the fixed amplitude of the MAPS experiments begins to drive a strongly, rather than weakly, nonlinear response, resulting in the prediction of unphysical values of the noise temperature.Thus, these experiments are discarded, as measurements of the linear rheology alone (as shown in Fig. 7) are sufficient to monitor the material beyond this point.
Before 45 min, we observe that the linear viscoelastic data alone (red) are only sufficient to monitor values of x below approximately 3.5.This is consistent with the qualitative predictions of the SGR model.Before 35 min, the predicted values of x from linear response data are highly uncertain, with a 50% confidence interval typically spanning values of 3 < x < 5.No substantial change in the linear viscoelastic response is detected until between 35 and 40 min, at which point the noise temperature transitions into the x < 3 regime, wherein the linear rheology becomes qualitatively consistent with that predicted by the fractional Maxwell model.Using the measured MAPS signature, however, the noise temperature can be monitored to values exceeding x = 5.In fact, we observe that substantial changes in the noise temperature can be detected between 25 and 30 min, a full 10 min (approximately 0.2t g , with t g ≈ 50 min) earlier than transitions in the linear rheology.That the behavior of the noise temperature inferred from MAPS experiments appears to follow our expectations from physical considerations, and eventually converges onto the well-characterized behavior inferred from linear rheology, also provides confidence that the MAPS tests have indeed measured intrinsic nonlinearities from the clay dispersion and that the tests do not disrupt the delicate aging microstructure.
In addition to the earlier detection threshold, the values of x predicted using MAPS data have smaller uncertainty than those predicted using linear response data alone, even for x < 3. Thus, monitoring the state of the fluid using the MAPS protocol is apparently more robust, even in the regime where linear data does exhibit some sensitivity to the noise temperature of the system.This is particularly clear in Fig. 9b), which depicts the Gumbel distribution inferred by each ANN for two times: t = 30 min and t = 40 min.In the case of the distributions inferred from the linear response alone (red curves), although at t = 40 min the mode falls below the detection threshold of x = 3, there remains substantial overlap at that time with the distribution at t = 30 min.It is therefore difficult to state with certainty that the linear response of the system has undergone a transition.The MAPS-inferred distributions, on the other hand, are much narrower and exhibit very little overlap between t = 30 min and t = 40 min.This provides a high level of confidence in the detection of the rheological transition, which evidently cannot be achieved with such certainty using linear measurements alone.
The results presented in Fig. 9 represent a significant advance for potential health monitoring applications in rheology.Not only does it indicate that weak nonlinearities can be used in practice to anticipate transitions in the linear response of complex fluids, but this realization is also based on first principles models rather than on empirical evidence alone.Thus, one might expect the same conclusions to hold true for other materials undergoing phase transitions that are qualitatively described by the SGR model, and the combination of MAPS rheology and machine learning may be similarly applied to health monitoring or fault detection in those systems.These results also represent the development of machine learning tools that may be employed in future online monitoring of gelling dispersions such as Laponite ® or other discotic clays.Indeed, all that is required to implement this method in practice is the online experimental capabilities to conduct two MAPS (at distinct stress amplitudes), along with the analysis software to obtain the MAPS response functions and to infer x from a trained ANN.We provide an example of this software, which was used to obtain the results in this work, on GitHub (https://github.com/krlennon/rheology-health-monitoring).Although MAPS experiments have only been conducted in a laboratory setting at present, the recent development of devices such as immersed microelectromechanical oscillators may enable online measurements of intermodulation effects within complex fluids (Gonzalez et al. 2018).With these tools, it may soon be possible to continually monitor the evolution of gelling or vitrifying materials within industrial processes.

Conclusions
We have demonstrated that the weak intermodulation nonlinearities measured by MAPS rheology are sensitive predictors of a qualitative change in the linear response of soft glassy materials.The sensitivity to model parameters in the weakly nonlinear response is strong enough that an artificial neural network is able to robustly learn an inverse function mapping the results of a three-tone MAPS experiment onto the model parameters x, f , and k with high confidence.In the weakly nonlinear limit, where the effects of the parameters c and d are nearly indistinguishable, the three parameters x, f , and k fully specify the behavior of the tensorial, frame-invariant SGR framework.This suggests that a minimal data set collected non-destructively which probes both the linear and weakly nonlinear response of soft glassy materials could be used to fully characterize the rheological response, not just to simple shear, but to any medium amplitude deformation protocol.
A more modest application of using ANNs to learn the inverse mapping between MAPS data and the SGR model parameters is the identification of incipient transitions in the qualitative features of the linear viscoelastic response of the gelling material by classification of the rheological state of the fluid in terms of the noise temperature, x.We have shown how robust inference of x allows one to identify clearly the boundaries for the transition to a gel-like rheological regime-from Maxwell-like to fractional Maxwell-like-as well as the transition to a vitrified rheological regime-from fractional Maxwell-like to fractional Kelvin-Voigt-like.We finally applied this inference to a real fluid system-an aqueous dispersion of Laponite ® clay platelets-and successfully identified transitions in the nonlinear rheology 10 min (≈ 0.2t g ) in advance of measurable changes in the linear rheology.This demonstrated ability to correctly identify the rheological state of a fluid non-destructively may prove useful for deciding whether a soft glassy material meets the specifications for a particular field application.
While the perspective taken in this work is oriented around detailed analytical calculations with a well-known, if specific, constitutive model, it suggests broader opportunities for the characterization of viscoelastic fluids in rheology.As indicated by the present study, there likely exists a similar sensitivity to changes in rheological state in many distinct sub-classes of soft materials.For example, careful study of the leading order of a power-law expansion in weakly nonlinear tests on capillary suspensions provides insight into the force laws connecting the bridging particles (Natalia et al. 2020(Natalia et al. , 2022)).Although the underlying constitutive behavior may not always be known precisely for many materials, monitoring of the weakly nonlinear response coupled to purely data-driven methods of regression or classification may be able to accomplish the same tasks performed in this work with ANNs trained purely on synthetically generated model data.

Outlook on machine learning in rheology
The role of machine learning methods-here, artificial neural networks-in this work has been to learn an approximation of the inverse of a computationally expensive rheological constitutive equation from synthetic data, thereby allowing one to rapidly infer model parameters from experimental data.Model inversion and acceleration are well-studied applications of machine learning in many fields, including recent endeavors within the field of rheology (Mahmoudabadbozchelou and Jamali 2021; Mahmoudabadbozchelou et al. 2022).A number of previous approaches have employed physicsinformed neural networks (PINNs), which seek a machine learning surrogate for a constitutive model by augmenting the traditional data residual error with the residual between the neural network predictions and the predictions of the constitutive equations.This technique essentially produces an accelerated solver for the constitutive equation.Recent work has adapted these PINN approaches for parameter identification (Saadat et al. 2022); however, this requires repeated evaluation of the model equations during inference.As we have discussed, off-loading the computational burden to the training of a machine learning model, rather than the repeated evaluation of an expensive model during parameter inference, may enable the online and real-time rheological monitoring of soft materials in industrial environments.Similarly, machine-learning (ML)-accelerated scientific simulations may expedite the design and engineering of materials and industrial processes.
In such applications of machine learning, however, we must recognize that the accuracy of either the inferred material properties (in this work, presented in the form of learned model parameters) or the scientific simulations is limited by that of the underlying constitutive equations describing the material.More robust engineering predictions of soft material systems often require the models themselves to be more accurate.Here, rheology stands to make substantial gains from the adoption of machine learning-in particular because rheological constitutive equations that are simple enough to be embedded within engineering pipelines tend not to be highly accurate, while the most accurate models are often too complex to evaluate numerically within simulation or inference workflows.
However, simply replacing rheological constitutive equations with "end-to-end" machine learning models, eschewing well-established rules and domain knowledge regarding the structure of admissible constitutive equations, is likely insufficient to generate accurate models in the face of limited rheological data.While imposing such knowledge at training time, for instance, by embedding physical laws as terms within the loss function (Raissi et al. 2019), can lessen the data burden, it does so at the expense of an increase in the training required for convergence and does not guarantee that the model will obey these laws when presented with test data.Instead, it is the current authors' view that the synthesis of machine learning methods and existing frame-invariant constitutive equations will facilitate the most rapid and robust advancements in rheological constitutive modeling.For instance, by embedding neural networks within differential constitutive equations-a so-called "universal differential equation"-rheologists can develop extensible surrogates that accelerate known constitutive equations (Rackauckas et al. 2020).Such frameworks can be automatically constructed to follow well-known rheological principles, such as the concepts of rheological invariance established by Oldroyd (Oldroyd 1984;Lennon et al. 2022), and computational effort may therefore be dedicated to learning the nonlinear features particular to a certain material.This paradigm in machinelearning-augmented rheological modeling produces models that are at once computationally efficient and highly accurate, while remaining general enough to suit the many goals of soft materials science and engineering.

Appendix 2. ANNs for other MAPS protocols
Previous work has explored experimental design within the MAPS rheology framework, with investigations of the relative merits of stress-versus strain-controlled experiments and the effect of the input tone set {n 1 , n 2 , n 3 }.The stress-controlled mode with tone set {5, 6, 9} was selected in this work due to practical considerations, such as the high-frequency inertial limitations of the stress-controlled rheometer; however, other experimental setups may admit other protocols.Here, we demonstrate that the same MAPS In each case, the data generation, processing, augmentation, and training procedures, as well as the architecture of the ANNs, are identical to those described in Sect.4(Table 2).

Fig. 2
Fig.2Logarithmic slopes of the linear and weakly nonlinear response functions of the SGR model in a strain controlled experiment-moduli:G * (ω), G * 3 (ω 1 , ω 2 , ω 3 )-in the low-frequency limit: ω, |ω 1 | + |ω 2 | + |ω 3 |1, as a function of noise temperature x.There exists a distinct regime (shaded) in which the logarithmic slopes of G 3 and G 3 serve as antecedents of a transition in the character of the linear response

Fig. 6
Fig. 6 Parity plots for the predicted versus true values of the SGR model parameters (ln k, ln f , x, c, and d) for both the ANN trained on linear response data alone (red triangles) and the ANN trained on MAPS

Fig. 9 a
Fig. 9 a The noise temperature inferred by ANNs from experimental data characterizing an aqueous dispersion of Laponite ® clay platelets.Values inferred from linear response (LR) data alone are shown in red, and values inferred from both MAPS and linear response data are shown in blue.A transition to a fractional linear response (x < 3, c.f., Fig. 4) occurs between 35 and 40 min, while a transition to fractional MAPS behavior (x < 5) occurs between 25 and 30 min, providing the early Pochhammer symbol defined by the following:

Fig. 10
Fig. 10 Parity plots for the predicted versus true values of the SGR model parameters (ln k, ln f , x, c, and d) for both the ANN trained on linear response data alone (red triangles), and the ANN trained on

Table 2
The activation function a n (.) and size of each layer l n for the two ANNs explored in this section, for MAPS experiments with a tone