Large amplitude oscillatory shear study of a colloidal gel near the critical state

The linear viscoelastic (LVE) properties of soft materials give crucial insights into the microstructure of the same.^1–5 Understanding how the latter affects the former is vital not only just from an academic point of view but also from an industrial processing perspective.^6–8 Rheological behavior under nonlinear flow field associated with industrial processes cannot be described by merely the knowledge of linear viscoelastic response. The strains that a material is subjected to during a practical application can be quite large. Consequently, understanding microstructure specific linear—nonlinear viscoelastic transition in soft materials is of critical importance. In this work, we analyze linear—nonlinear viscoelastic transition of the space spanning percolated fractal network of colloidal particles (colloidal gel at the critical state) using a large amplitude oscillatory shear (LAOS) approach.

The investigation of soft materials by applying large amplitude oscillatory shear dates back to the early 1960s.⁹ The measurement of a response (the stress or the strain) was conducted in an analog fashion, where the nonlinearity was represented by a deviation in the closed-loop plots (Lissajous–Bowditch curves) of stress vs strain or stress vs strain rate.^10,11 The graphical way of representation of the measured data by Lissajous–Bowditch curves has been a major characterization feature of LAOS.^10,12 Matsumoto et al.¹³ were the first to quantify the nonlinearity in terms of the third harmonic of moduli for a particulate suspension in polystyrene solution. Giacomin and Oakley¹⁴ further analyzed the Lissajous-Bowditch curve to obtain the Fourier series from it. The Fourier transform (FT) analysis of the stress data had become one of the prevailing methods of quantifying the nonlinear viscoelasticity.^{11,13,15–17} This led to the birth of the term “Fourier transform (FT) rheology.”¹⁸ Since then, it has been routinely used to study the nonlinear rheological behavior of soft materials.^19–22 The FT rheology framework is mathematically robust. Especially, with the knowledge of higher harmonics, the output wave can be reconstructed by filtering out the noise. Analysis through FT rheology can lead to profound insights into the material behavior, especially in distinguishing one material from the other.²³ Nevertheless, the FT rheology leads to the limited physical interpretation of the higher harmonic coefficients.²⁴ In order to overcome this shortcoming, Cho et al.²⁵ developed the stress decomposition (SD) method, wherein the total stress is decomposed into elastic and viscous components. The SD methodology has been proven to be mathematically equivalent to FT rheology.²⁶ Ewoldt et al.²⁷ replaced the unknown material functions in the SD method with orthogonal Chebyshev polynomials and coefficients. Later, Hyun and Wilhelm²⁸ introduced a nonlinear parameter known as the Q parameter, which is expressed in terms of the ratio of intensities of the first and third harmonic. The shape of the Q curve has been analyzed to distinguish between linear and branched polymers.⁵ Rogers and co-workers^24,29 suggested a new qualitative analysis and described the response of yield stress materials as a sequence of physical processes (SPP) by introducing the locally defined measures including time dependent dynamic moduli. However, detailed information at every time step throughout the cycle is required for SPP to describe the material response.^24,30 Furthermore, SPP approach is a high-dimensional representation even in the linear viscoelastic regime and is observed to inaccurately report elasticity for a generalized Newtonian fluid model, which makes it a challenging approach.^27,30 Recently Erturk et al.³¹ compared SPP approach and Chebyshev polynomials approach to analyze nonlinear viscoelastic behavior of wheat flour dough and pectin solution. Over the past two decades, various approaches have been employed to study LAOS response to various kinds of materials. How LAOS parameters can be analyzed to obtain meaningful information about microstructure, and how it gets affected by nonlinear deformation has been reviewed by Hyun et al.³² and more recently by Cho.¹⁸ Lately, there has been a greater thrust on analyzing intrinsic nonlinear complex moduli, which depend only on frequency.^33–36 The nonlinear complex moduli can be obtained by expressing the stress waveform in response to sinusoidal strain input as a power series.^37,38 For medium amplitude oscillatory shear, Ewoldt and co-workers²⁷ analyzed nonlinear complex moduli and proposed that their signs give physically meaningful interpretations. However, experimental measurement of these intrinsic parameters requires fine resolution and high precision at low torque limits.

Apart from experimental characterization, there have been theoretical efforts using mode-coupling theory (MCT) to understand the behavior of material upon large deformations. In order to extend the mode-coupling theory to model the system behavior under flow, few additional approximations beyond those of quiescent MCT are required.³⁹ One of the approaches involves an integration through transients (ITT) formalism, originally developed by Fuchs and Cates.⁶ ITT approach has been subsequently developed and has been comprehensively studied in case of steady shear,⁶ time-dependent shear⁴⁰ and time dependent flows.⁴¹ Furthermore, MCT-ITT approach has been utilized to study the nonlinear response of dense colloidal suspensions upon application of a large amplitude of oscillatory shear.⁴² The MCT-ITT theory predicts the final power-law decay of moduli with increasing strain amplitude. However, despite the overall qualitative agreement of the MCT-ITT theory with the experiments, there is scope for improvements in the theory, which would lead to a better quantitative agreement with experiments (see Ref. 42 for details).

In this work, we study LAOS behavior of the unique state called the critical gel state in a colloidal gel. Sol–gel transition, wherein material transforms from a free-flowing liquid state to a non-flowing state with network-like microstructure, has been observed for polymeric^43,44 as well as colloidal systems.^45,46 Such sol–gel transition could be spontaneous,⁴⁷ because of change in temperature⁴³ or concentration of one of the species.⁴⁸ It has been observed that in a limit of free flowing liquid state, elastic (G′) and viscous (G″) moduli show a terminal behavior given by G′ ∼ ω² and G″ ∼ ω, where ω is the probed frequency associated with the linear response regime.⁴⁹ As the extent of crosslinking increases, ω dependence of both the moduli decreases, with that associated with G″ decreasing at a faster rate. Eventually growing clusters connect with each other such that they just form a space spanning percolated network. Such weakest network has been represented as a critical gel state in the literature.⁴⁹ The critical gel has a unique signature wherein G′ and G″ are observed to show an identical power law dependence on ω given by⁴⁹ 

where n_c is the power law index (0 < n_c < 1), S is the gel strength, while Γ(⋯) represents the Euler Gamma function. Knowledge of the power law index n_c leads to fractal dimension (f_d) of the percolated network $fd=52nc−3/2nc−3$⁠.⁵⁰ Remarkably, this relationship between n_c and f_d has been validated for fractal gels that were subjected to simultaneous rheology and scattering studies.^51,52 Nonlinear viscoelastic behavior of the critical gels has been investigated by performing step strain experiments at different magnitudes of strains in the linear and nonlinear regime.^53,54 The resultant self-similar nature of the relaxation modulus leads to damping function, which has been used to analyze linear—nonlinear transition in the critical gels. Suman and Joshi⁵⁴ reported strain softening for colloidal clay gel while Keshavarz et al.⁵³ reported the strain hardening followed by strain softening behavior for casein network formed through physical interactions. Nonlinear rheological behavior at the critical state has also been studied for a food gel⁵⁵ as well as a polymer gel.⁵⁶ Gisler et al.⁵⁷ investigated diffusion limited aggregate of fractal colloidal gel and reported strain hardening originating from the rigid backbone of the same. However, a comprehensive experimental characterization and theoretical modeling of the nonlinear viscoelastic behavior of the uniquely defined critical gel state is still in its infancy.

In this work, we present a detailed analysis of nonlinear measurements of a colloidal dispersion of synthetic hectorite clay at the critical gel state. Based on our observations, we propose a picture of the microstructural changes occurring during the transition from linear to nonlinear regime. We also model the nonlinear behavior using an integral constitutive equation that uses information of linear viscoelasticity and damping function obtained from nonlinear stress relaxation modulus. The proposed model with further modification of the damping function predicts the intracycle stress response in the linear and nonlinear region quite well and provides physical insights during the nonlinear transition.

In order to characterize the linear and nonlinear response, the material is subjected to a sinusoidal strain input with amplitude γ₀ at a constant frequency (ω) given by

The linear stress response on the application of a sinusoidal strain input can be represented as⁵⁸ 

where G′ and G″ are elastic and viscous moduli, respectively. Furthermore, on eliminating time variable in the above equations, we get the direct relationship between stress and strain given by

The above equation suggests the Lissajous-Bowditch curve (stress–strain curve) corresponds to a straight line for a linear elastic material where stress is in phase with strain. On the other hand, stress is out of phase with strain for purely viscous material and the Lissajous-Bowditch curve is represented by a circle. Importantly, for a linear viscoelastic material, the Lissajous-Bowditch curve is elliptical.

In the nonlinear region, G′ and G″ are not uniquely defined since stress response is not a single harmonic sinusoid.⁵⁹ The most common method to access the nonlinearity is Fourier transform rheology.⁶⁰ The stress response on the application of a sinusoidal strain input can be represented by a Fourier series as⁵⁸ 

where $Gn′$ and $Gn″$ are the nth order elastic and viscous moduli, respectively. The Fourier series contains only odd harmonics in order to obey the directionality of strain, which suggests that on reversing the coordinate system, the material response remains unchanged. The stress response consists of first harmonic only in the limit of linear viscoelasticity and the material properties are characterized by G′ and G″, while the higher harmonic appears with an increase in the strain amplitude in the nonlinear region. Even though the FT framework is mathematically complete, it fails to provide a physical interpretation of the higher-order coefficients.

In order to interpret the LAOS data meaningfully, Ewoldt et al.²⁷ adopted the method of orthogonal stress decomposition proposed by Cho et al.²⁵ The idea was to decompose the stress response into a summation of elastic stress $σ′$ and viscous stress $σ″$ using symmetry arguments as

and

The analysis of the nonlinear parameters underscores the microstructural changes occurring in the material on application of the deformation field in the nonlinear domain.

In this work, we prepare a dispersion of 3 wt. % synthetic hectorite clay LAPONITE XLG^® having 3 mM NaCl. Hectorite clay (LAPONITE XLG^®) consists of disk shaped particles of around 25–30 nm in diameter and 1 nm in thickness. On dispersing the LAPONITE particles in water, the negatively charged faces forms attractive bonds with the positively charged edges and leads to the formation of a network structure. The comprehensive characterization of LAPONITE XLG^® as a particle and detailed protocol to prepare the dispersion is discussed in our previous publication.⁵⁴ Briefly, oven dried (120 °C for 4 h) white powder of hectorite was added to ultrapure millipore water with added salt and stirred using ultra Turrax drive for 30 min. Hereafter, we refer to this dispersion as hectorite dispersion. The freshly prepared dispersion was loaded in the sandblasted concentric cylinder geometry having a cup diameter of 30.38 mm and a gap of 1.17 mm on TA Instruments DHR 3 rheometer for rheological measurements. Subsequent to the loading of the sample, it was subjected to cyclic frequency sweep within the linear viscoelastic regime with a stress magnitude of 0.1 Pa and frequency range of 0.5–25 rad/s. It takes 100 s to complete one frequency sweep. The amplitude sweep and all the nonlinear measurements were performed at the critical gel state of the colloidal dispersion. In a large amplitude oscillatory shear (LAOS) test, the critical gel is subjected to a sinusoidal strain of the form: γ = γ₀ sin ωt where ω is a constant frequency and γ₀ is the amplitude of strain varying logarithmically from 0.001 to 10. The LAOS tests were performed at various angular frequencies of 0.63, 3.14, and 31.4 rad/s, which gets completed within 1500 s. We did not explore higher strain amplitudes and frequencies to avoid fracture of sample at the free surface. The transient data were collected using the rheometer software Trios. For every strain input, the raw values of 900 data points were recorded in each cycle. The smoothing of the data and analysis of the stress and strain waveform data were performed using MITlaos program (Version 2.2 beta).⁶¹ The elastic and viscous stress components, Chebyshev coefficients, higher harmonic viscoelastic moduli, and Fourier spectrum were determined using MITlaos. The Fourier coefficient obtained from MITlaos program was also verified independently by performing Fourier transform using Matlab. All the experiments were conducted at a constant temperature of 30 °C maintained by a Peltier concentric cylinder temperature system. We employ a solvent trap and also add a thin layer of low viscosity silicone oil (320 mPas) to the free surface of the sample to prevent the evaporation of water during the measurements.

The freshly prepared colloidal dispersion is subjected to cyclic frequency sweep, which enables the estimation of viscoelastic properties at different extents of gelation. In Fig. 1(a), we plot the evolution of elastic $G′$ and viscous $G″$ modulus with respect to angular frequency $ω$ during the process of sol–gel transition. The modulus values at the respective times have been shifted vertically for clarity and the corresponding shift factors are mentioned in the legend. At very small times, G″ dominates G′ and the dispersion exhibits characteristic terminal behavior with G′ ∼ ω² and G″ ∼ ω dependence, which is suggestive of the sol state. In Fig. 1(b) we plot the evolution of loss tangent $tanδ$ as a function of time at different frequencies, which also shows a prominent decrease with ω. However, as time progress, both the moduli grow with the rate of growth of G′ being higher than G″. Eventually, G′ surpasses G″ and the dependence of G′, G″ and tan δ on ω weakens with an increase in time. At a particular time, both the moduli exhibit identical dependence on ω $G′∼G″∼ω0.27$ while tan δ becomes completely independent of the applied frequency. Such dependence is the characteristic rheological signature of the critical gel transition. The time to achieve the critical gel transition is indicated by a solid arrow in Fig. 1(b). Very interestingly, the critical relaxation exponent $nc$ computed at the critical gel state using the relation n_c = 2δ/π is identical to the power-law exponent of viscoelastic moduli on ω [As suggested by Eq. (1)]. At further higher times, G′ and G″ are only weakly frequency-dependent while the dependence of tan δ on ω inverts and eventually tan δ increases with an increase in ω. This is a distinctive feature associated with the post-gel state. Therefore, Fig. 1 clearly suggests that the results exhibited by a freshly prepared dispersion of synthetic hectorite clay are consistent with a system undergoing sol–gel transition.

We now aim to study the nonlinear rheological behavior of the material at the critical gel state. In order to ensure the colloidal dispersion is at the critical gel state, we first subject the freshly prepared dispersion to cyclic frequency sweep experiment. At the point when tan δ becomes frequency independent, we subject the critical gel to an oscillatory strain sweep. The time window in which the oscillatory strain sweep is performed is indicated by the shaded region in Fig. 1(b). The critical gel state is the weakest space spanning fractal network as it is on the verge of a pre-gel state. In Fig. 2, we plot the variation in G′ and G″ as a function of strain amplitude (γ₀) at three different values of explored frequencies. The frequency values are chosen in such a way that the system remains in the critical gel state over the complete duration of strain sweep. At low strain amplitude, as expected, the critical gel exhibits constant moduli. The value of G′ is greater than G″ in the linear viscoelastic region, by virtue of n < 0.5, thus indicating the critical gel to be dominated by elastic effects. In this regime, G′ and G″ increases with an increase in the magnitude of the applied frequency. The linear viscoelastic regime is observed until γ₀ ≈ 0.1 for all the frequencies after that the value of G′ starts to decrease as the critical gel enters the nonlinear viscoelastic region. It can be seen in Fig. 2 that both G′ and G″ decrease with an increase in the strain amplitude in the nonlinear regime. Since critical gel state is the weakest percolated network, the energy dissipated during the linear—nonlinear transition is less, which does not lead to any peak in G″.⁶² In the nonlinear regime, G′ and G″ lose their original meaning, and represent $G1′$ and $G1″$⁠. In the limit of large γ, the decay in $G1′$ and $G1″$ follows a power law given by $G1′∼γ0−α$ and $G1″∼γ0−β$⁠. The power-law exponents are observed to be independent of the applied frequency and take the value as α = 1.62 ± 0.015 and β = 0.74 ± 0.035 for all the explored values of frequency. Similar value of power-law exponents has been reported for Hectorite-PEG dispersion⁶³ and pluronic-hyaluronic acid gels.²² The critical gel shows a similar qualitative response at all three values of the explored frequencies.

The transient data collected during an oscillatory strain sweep is analyzed to obtain various linear and nonlinear viscoelastic parameters. We first analyze the transient data using Fourier transform (FT) Rheology. The viscoelastic nonlinearity can be identified by the appearance of high-order harmonics on the Fourier transform spectra of the stress in the frequency domain. In Fig. 3, we plot the intensity of the nth harmonic related to the primary one (I_n/I₁) at various strain amplitudes ranging from the linear to nonlinear region with a primary frequency of ω = 3.14 rad/s. The noise floor appears near a normalized intensity of I_n/I₁ ≈ 2 × 10⁻⁴. The signal-to-noise (S/N) ratio was computed following the definition by Merger and Wilhelm,⁶⁴ with S = 1 and N being the average of the normalized spectra data in the harmonic window of 15.2 < n < 16.2, where ideally no peak should be present. For all the values of strain amplitudes, S/N ratio was always above 10³, which substantiates the reliable detection of the harmonic intensities. For smaller strain amplitudes belonging to the linear viscoelastic region (γ₀ = 0.001 and 0.01), there is a primary peak at n = 1 while all higher harmonics are negligible (I_n/I₁ ≈ 2 × 10⁻³). For the strain of γ₀ = 10⁻² the peak associated with the second harmonic can be seen to be as strong as that of third harmonic. However, both the peaks are around two orders of magnitude smaller than that associated with the primary harmonic indicating the nonlinear effects are sufficiently weak to be negligible. Such low value of even harmonics suggests the effect of wall slip^65,66 or secondary flows⁶⁷ to be small and not very important as compared with the nonlinear effect demonstrated by presence of odd harmonics. With increase in the strain amplitude to γ₀ = 1.1, the intensity peaks at higher odd harmonics become dominant, thus signifying the transition into the nonlinear region. At the highest strain amplitude of γ₀ = 5.8, the nonlinearity is more pronounced by the peaks at high odd harmonics.

In Fig. 4(a) we plot the variation of G′ and G″ as a function of strain amplitude at ω = 3.14 rad/s. The linear to nonlinear transition can be seen in Fig. 4(a) with an increase in the strain magnitude. In Fig. 4(b), we show the intracycle response by plotting the total stress as a function of strain for the strain magnitudes highlighted in Fig. 4(a) with identical symbols. It can be seen in Fig. 4(b), that the stress induced in the critical gel at low amplitudes of strain is small, which as expected increases with an increase in the strain amplitude. The response of the critical gel to the strain loading in the linear viscoelastic region is primarily elliptical as suggested by Eq. (4). Furthermore, the area enclosed by the stress–strain closed curves increases with increasing strain amplitude. At large strain amplitude, as the system gets into the nonlinear regime, the shape of the stress–strain curve becomes qualitatively different. Using the framework discussed by Cho et al.²⁵ and Ewoldt et al.,²⁷ the total stress can be decomposed into the elastic and viscous contributions, which are shown in Figs. 4(c) and 4(d), respectively. While the total stress response is a hysteresis loop in the stress–strain plane, the elastic and viscous stresses are single-valued functions, when plotted, respectively, as a function of γ and $γ̇$⁠. It can be seen in Figs. 4(c) and 4(d) that with an increase in the strain amplitude, the nonlinearity is evident through the deformation of the stress–strain curves. While the elastic stress curves plotted against γ turn upward, the viscous stress plotted as a function of $γ̇$ can be seen to curve downward at the highest amplitude of strain.

The nonlinearities can also be examined graphically from the 3D Lissajous-Bowditch plots in Fig. 5. The projections of the Lissajous-Bowditch curves on stress–strain, stress-shear rate, and strain-shear rate planes are shown in Fig. 5. In an oscillatory cycle, the variation of strain is orthogonal to the strain rate. Consequently, the strain achieves the maximum value when the strain rate is zero and the other way around. While the stress–strain curves are elliptical in the linear viscoelastic region, the shape deviates from a typical elliptical response as the flow field enters the nonlinear regime. As the stress–strain Lissajous-Bowditch curve approaches a rectangular shape at high amplitude of strain field, the stress–strain rate Lissajous-Bowditch curve bend over at maximum strain rate. It is interesting to note that even though the oscillatory strain sweep does not show any variation in the viscoelastic properties in the linear domain, the Lissajous-Bowditch curve associated with each strain amplitude in the linear domain is different. The unique shape associated with all the strain amplitudes sets an unambiguous measure of quantifying the linear as well as the nonlinear response of the colloidal dispersion. Therefore, the Lissajous-Bowditch curves act as a rheological fingerprint test of a material. However, it is important to note that Lissajous-Bowditch curves by itself, gives a qualitative estimation of the material behavior at different strain amplitudes and further quantitative assessment is required for an exact estimation of nonlinearity.

In order to analyze the shape of the Lissajous-Bowditch curve and account for the nonlinearities in the stress response, the following new parameters are defined. The minimum-strain $(GM′)$ and large-strain moduli $(GL′)$ are defined as²⁷ 

These moduli are defined in such a way that each of the moduli reduces to G′ in the limit of linear viscoelasticity, which, in principle, represents the first harmonic of elastic modulus in the linear viscoelastic domain $(G1′)$⁠. The graphical interpretation of nonlinear moduli at a strain amplitude in the nonlinear region is shown in Fig. 6(a). It can be seen that $GM′$ is the tangent modulus at minimum strain amplitude while $GL′$ is the secant modulus at maximum strain amplitude. The nonlinear moduli quantify the extent of distortion within a cycle. The first harmonic of elastic modulus $(G1′)$⁠, minimum strain modulus $(GM′)$ and large strain modulus $(GL′)$ are plotted in Fig. 6(c) as a function of strain amplitude. At low amplitudes of strain (linear viscoelastic region), both the moduli $GM′$ and $GL′$⁠, reduce to $G1′$⁠, where the intensity ratio I₃/I₁ is below 0.3% and Lissajous-Bowditch curve has a perfectly ellipsoidal shape. As expected, the increase in the magnitude of strain, at the onset of nonlinearity, results in a decrease in $GM′$ and $GL′$⁠, such that in the nonlinear region, $GM′$ and $GL′$ are found to deviate dramatically from $G1′$⁠. This deviation of $GM′$ and $GL′$ from $G1′$ signifies the divergence of a Lissajous-Bowditch curve from a linear viscoelastic response. The increasing extent of deviation of nonlinear moduli from $G1′$ with increase in strain amplitude suggests a greater extent of distortion in the Lissajous-Bowditch curves. Furthermore, $GM′$ decreases at a faster rate than $GL′$ as evident at higher values of strain amplitude. This suggests that within a particular strain cycle at constant strain amplitude, $GL′$ is greater than $GM′$⁠. The dominance of $GL′$ over $GM′$ suggests that the system is able to store more energy at the highest amplitude of strain than at the minimum value of strain. This advocates toward intracycle strain stiffening behavior in the critical gel state.²⁷ However, it is important to note that, overall, the value of moduli is decreasing with an increase in the strain amplitude, which suggests overall strain thinning behavior in the critical gel state.

In order to quantify the dissipative response of the material in the nonlinear domain, a set of dynamic viscosities have also been defined as²⁷ 

where $ηM′$ is the minimum-rate dynamic viscosity and $ηL′$ is the large-rate dynamic viscosity. In the limit of linear viscoelasticity, both the viscosities reduce to $ηM′=ηL′=η1′=G″1/ω$⁠. The nonlinear viscosities are graphically represented in Fig. 6(b). The viscous nonlinearities are plotted as a function of strain amplitude in Fig. 6(d). As expected, $ηM′$ and $ηL′$ remain identical to $η1′$ in the linear viscoelastic region. In the nonlinear region, $ηM′$ and $ηL′$ decreases with the rate of decay of $ηL′$ being higher than $ηM′$⁠. This suggests that the value of $ηM′$ is always higher than $ηL′$ in the nonlinear region. Such a dependence of higher viscosity at the minimum strain rate than at highest strain rate is an indicator of intracycle shear thinning behavior of the critical gel state.²⁷ We also analyze behavior of Chebyshev coefficients as a function of amplitude of strain amplitude at various exploited frequencies. The corresponding behavior is discussed in the supplementary material. The inference of the Chebyshev coefficients analysis matches well with that reported in Fig. 6.

The elastic Lissajous-Bowditch curves at eight different strain amplitudes and varying frequencies represented by solid black lines are shown as a Pipkin diagram in Fig. 7. This diagram is defined using two independent parameters: the strain amplitude on the vertical axis and the frequency on the horizontal axis. Each Lissajous-Bowditch curve within the Pipkin diagram represents a plot of normalized stress on the vertical axis with normalized strain on the horizontal axis. The normalization of strain is carried out using the amplitude of imposed strain (γ₀ indicated on the vertical axis of Pipkin space) and stress is normalized using the maximum value of stress achieved in each cycle (indicated at bottom right corner of or inside each Lissajous-Bowditch curve). The local shape of the Lissajous-Bowditch curves depends on all the harmonics present in the system. For the linear viscoelastic region in Fig. 7 (γ₀ < 0.1), the stress–strain curves are elliptical at all the explored frequencies. As strain amplitude increases (γ₀ > 0.1), Lissajous-Bowditch curves deviate from an elliptical shape and get distorted as higher-order harmonics become non-zero. The elastic Lissajous-Bowditch curves broaden in the center, and the area enclosed by a Lissajous-Bowditch curve increases with an increase in strain amplitude. The qualitative effects of increasing strain amplitude are similar at all the explored frequencies. Furthermore, the maximum stress attained in each cycle increases as the oscillation frequency increases. The weak change in Lissajous-Bowditch curves upon increase in the frequency is a consequence of the weak increase of the moduli with frequency with a power-law exponent of 0.27 for the critical gel state. Such weak change in Lissajous-Bowditch curves across the frequencies is also reported in critical gel of gluten,⁶⁸ and colloidal glass of star-like micelles system.⁶⁹ 

To model the behavior of the critical gel state, we conducted stress relaxation experiments over a strain range of γ₀ between 0.002 and 1. A freshly prepared sample is used in each stress relaxation experiment, and the step strain is applied at the point of critical gel transition. The resulting decay of the stress relaxation modulus (G) is shown in Fig. 8(a). In the linear viscoelastic limit, G for a critical gel exhibits a power-law dependence on time given by⁷⁰ 

where S is the gel strength. Consequently, the experimental relaxation modulus data are independent of the amplitude of strain within the linear viscoelastic region (γ₀ ≤ 0.1), and overlap on top of each other. The corresponding power-law index is identical to the value of n_c obtained from the oscillatory studies as plotted in Fig. 8(a). Interestingly, as shown in Fig. 8(a), for the applied step strains beyond the linear regime, the stress relaxation curve obeys the same power-law dependence given by Eq. (12). However, the magnitude of the relaxation modulus decreases with an increase in γ₀. The strain-dependent stress relaxation modulus is represented as G(t, γ₀). The self-similar behavior in G(t, γ₀) over the explored duration of γ allows vertical shifting using a damping function h(γ₀). This leads to the superposition in G(t, γ₀) as shown in Fig. 8(b). The existence of the superposition in G(t, γ₀) suggests that the strain-dependent relaxation modulus can be factorized into time-dependent $Gt$ and strain-dependent h(γ₀) as⁵⁶ 

The vertical shift factor h(γ₀) is plotted as a function of strain in Fig. 8(c). As expected, the value of h(γ₀) is seen to remain close to 1 in the linear viscoelastic region. With an increase in the amplitude of strain, the value of h(γ₀) decreases below 1. In the literature, a number of analytical expressions of exponential, sigmoidal, and other types have been proposed for h(γ₀).⁷¹ The data shown in Fig. 8(c) can be described by a popular functional form of damping function given by^72–76 

where a is the model parameter. This functional form closely resembles the predicted damping function from the Doi-Edwards theory formulated with and without the independent alignment approximation.⁷⁷ It captures the deviation of damping function from unity at increasing strain beyond linear viscoelastic limit and approximates the predicted behavior of damping function given by the Doi-Edwards reptation theory.^76,77 We rigorously fit Eq. (14) to the data in Fig. 8(c) and determined the best-fit parameter a by minimizing the least squares difference. The best-fit a = 2.9 ± 0.6 that describes the experimental data well is shown by the dashed blue line. This damping function is valid up to γ₀ ≈ 1, based on the range of experimental data in Fig. 8(b).

To model the nonlinearity, we employ the time-strain separable (TSS) K-BKZ (Kaye–Bernstein–Kearsley–Zapas) model.⁷⁶ The integrand in the K-BKZ model can be factored into two components: the memory function m(t − t′) = ∂G(t − t′)/∂t′, and a purely periodic part $hγt−γt′γt−γt′$⁠, where γ(t) − γ(t′) = γ₀(sin ωt − sin ωt′) is the strain accumulated between times t and t′. The shear stress is given by¹⁸ 

The nonlinear response stems from the damping function. In the limit of small values of accumulated strain, the damping function is close to unity, and the time-strain separable K-BKZ equation becomes equivalent to the Boltzmann superposition equation. We incorporate the damping function given by Eq. (14) into the K-BKZ equation for an oscillatory strain input γ(t′) = γ₀ sin ωt′ to give the quasilinear constitutive equation:

Since all the material constants (S, n_c, and a) are known from LVE and step-strain experiments, the integral in Eq. (16) can be numerically solved for $σt$ once the strain amplitude and frequency are specified. However, on the application of oscillatory shear, numerical integration of TSS K-BKZ equation given by Eq. (16) becomes complicated due to the infinite domain of the integration and highly oscillatory nature of the strain at large frequencies.⁷⁸ 

Recently, we developed a spectrally accurate method for solving the time-strain separable K-BKZ model in LAOS,⁷⁸ which has the same form given by Eq. (15). The efficiency and accuracy of the spectral method over other methods of numerical integration have already been established.⁷⁸ The spectral method takes care of the infinite limit of integration and oscillatory integrand without any uncontrolled assumptions. By changing variable s = t − t′, Eq. (15) can be rewritten as⁷⁸ 

where $p(ωs;ωt,γ0)=hγt−γt′γt−γt′$ is a smooth periodic function of ωs with a period 2π/ω. We can accurately approximate p(ωs) with a discrete Fourier series (DFS), $Snf(ωs)$⁠, where the subscript n_f represents the highest harmonic resolved in the Fourier series,⁷⁸ 

In practice, n_f can be increased and automatically selected to ensure that the maximum deviation $p(ωs)−Snf(ωs)$ is below a pre-specified error tolerance, such as 10⁻⁸ used here. The DFS coefficients can be determined using fast Fourier transform with an asymptotic computational cost of $O(nflognf)$⁠.⁷⁸ When s = t − t′ = 0, $γt−γt′=0$⁠, which implies $p(ωs=0)=hγt−γt′γt−γt′=0$⁠, and hence $a0/2=−∑k=1nfak$ from Eq. (18). Using this relation, and substituting Eq. (18) into Eq. (17), we obtain an approximation for the shear stress $σ(t)≈σnf(t)$ as⁷⁸ 

Interestingly, the dynamic moduli in the LVE regime are related to the memory function via,⁷⁹ 

and

where the second equality in both expressions above is the dynamic moduli for critical gels. This allows us to analytically evaluate the shear stress as

The only approximation invoked in this method is the representation of p(ωs) with a DFS given by Eq. (18). Consequently, the convergence of this method is spectral, and is insensitive to large values of ω, which result in highly oscillatory integrals that degrade the performance of standard quadrature methods. Predictions using this spectral method, which avoids numerical quadrature, and conveniently side-steps the aforementioned singularity in the integrand, are shown by dashed blue lines in Fig. 7. At the smallest strain amplitudes, the agreement between these predictions and experimental data is excellent. The model prediction of the stress response is perfectly elliptical in the linear domain and fits the experimental data well at all explored frequencies as seen in Fig. 7. As strain amplitude increases beyond γ₀ = 0.1, predictions start to deviate systematically from the experimental data. Although the predicted stress response also exhibits broadening at the center, it underestimates its magnitude. Thus, we attribute most of the discrepancy between model predictions and experiments to the uncertainty in estimating the damping function, and the breakdown of the assumption of time-strain separability. As mentioned earlier, beyond a certain strain amplitude, the power-law behavior of the experimentally obtained $Gt,γ$ breaks down. Therefore, like most real materials, the critical gel empirically exhibits TSS over a finite strain amplitude window. However, the TSS K-BKZ equation imagines an ideal material that obeys time-strain separability at all strain amplitudes.

Nevertheless, it is perhaps still meaningful to ask the question, “if we assume the material to be time-strain separable, what damping function is implied by the Pipkin diagram shown in Fig. 7?” In other words, can we find parameters for the damping function that yield a more satisfactory fit of the experimentally observed Lissajous-Bowditch curves? Due to the efficiency of the spectral method, we can pose this as a least squares optimization problem by considering only a subset of the Lissajous-Bowditch curves in Fig. 7. Here, we consider the data corresponding to the lowest frequency (0.63 rad/s), i.e., the left-most column of the Pipkin diagram, and seek the model parameter (a) that simultaneously provides the best fit to all the Lissajous-Bowditch curves in this subset. This exercise yields the optimized value for the parameter: a = 2.1, which is less than the value (a = 2.9) obtained from the step-strain experiment. The value of a has been analyzed in the literature to relate to branching in polymers, and a smaller value is reported for branched topologies compared with linear polymers.^71,80 The modified damping function with smaller value of model parameter a = 2.1 is, therefore, suggestive of a denser percolated structure than reflected by the step strain measurements. A careful study combining rheology and scattering at a high strain field might be helpful in rendering insights about the microstructure.

The damping function with the modified parameter is shown as a solid red line in Fig. 8(c). Fits (0.63 rad/s) and predictions (3.14 and 31.4 rad/s) obtained using these optimized parameters and the spectral method are shown by solid red lines in Fig. 7. Interestingly, the damping function inferred from a single frequency works remarkably well at other frequencies as indicated by solid red lines in Fig. 7. It is worthwhile to mention that the prediction from the TSS K-BKZ equation with modified damping function matches the distortion in the experimental Lissajous-Bowditch curves exceedingly well in the nonlinear region. Similar approach was previously reported for polymeric liquids by Giacomin and Oakley.¹⁴ They also found that TSS K-BKZ model using the damping function determined from step-strain experiments produced poor predictions of the LAOS response. However, fitting the damping function to a subset of the LAOS data produced a satisfactory description of the response at other frequencies and strain amplitudes.

To investigate whether the difference in the values of the damping function parameter is significant, we perform Bayesian analysis to estimate the plausible range of values for a at each strain amplitude and frequency (see  Appendix A). We find that both a = 2.1 and a = 2.9 are within this plausible range for most of the datasets, which explains their effectiveness at describing the data. Nevertheless, we note that the predicted loops in Figs. 7 and 11 ( Appendix A) bulge beyond the experimental curves approximately along the major and minor axes, which are particularly visible at large values of strain amplitude. The nature of this discrepancy, and the difference in the damping function parameter obtained from step-strain experiments and LAOS measurements can be understood in terms of the difference in which the relaxation spectrum is modified under these flow conditions (see  Appendix B).^81–83 

The analysis of the large amplitude oscillatory shear and the resultant nonlinear parameters obtained from the transient data suggest that the studied material undergoes intracycle strain hardening, intracycle shear thinning, and overall strain softening. In this section, we propose the possible microstructural changes occurring in the critical gel state leading to above-mentioned phenomena. The schematic represented in Fig. 9 corresponds to a percolated network formed by the edge–face association between the colloidal particles. The percolated network remains unaltered at low amplitudes of strain and the response is viscoelastic without the presence of any nonlinear effects. The dashed rectangle highlights a particular section of the network, which we shall analyze when subjected to strain in the nonlinear region. The estimated elastic nonlinear parameters $(GM′andGL′)$ suggest intracycle strain hardening, while the viscous parameters $(ηM′andηL′)$ indicate intracycle shear thinning in the system. We believe that the intracycle strain hardening is likely due to deformation of the network of clay particles, which offers greater resistance to the elastic deformation due to the stretching of the same. The dashed rectangle in Fig. 9 representing the interparticle associations can be seen to get stretched upon the application of nonlinear strain deformation. Furthermore, the particles also get aligned in the direction of flow as indicated by the horizontal shift of the highlighted rectangle in Fig. 9(b), which results in intracycle shear thinning. At high amplitudes in the nonlinear region, the network junctions might break, and the broken clay clusters have little opportunity to rejoin the network structure under continuous deformation. This results in overall strain thinning behavior in the critical gel as evident by the decrease of moduli in a regular strain sweep test. Although it may seem paradoxical to have intracycle strain hardening and overall shear thinning simultaneously in a system, it is important to note that LAOS analysis of a single oscillation cycle suggesting strain hardening is a local effect while the regular strain sweep suggesting strain thinning is an overall behavior.⁸⁴ 

This work, therefore, presents a comprehensive characterization of a critical gel in the linear as well as in nonlinear domain. The TSS K-BKZ model describes the rheological behavior of the critical gel ranging from the linear to nonlinear domain. The damping function drives the nonlinear response of the model. Remarkably, the TSS K-BKZ model with optimized damping function parameters solved using an efficient numerical method describes the experimental data well. This work is a complete report on the experimental characterization and theoretical modeling of LAOS behavior of a critical gel. The proposed model describes the complete rheological behavior of the critical gel from linear to nonlinear domain. We believe that the simple model proposed in this work can be extended to a wide range of soft materials at the critical gel state to study the linear, weakly nonlinear, and nonlinear viscoelastic behavior.

We investigate the nonlinear behavior of synthetic hectorite clay at the critical gel state by estimating higher-order elastic and viscous moduli. Inspection of the intracycle data using Lissajous-Bowditch curves renders physical insights about the nonlinearity. The dominance of large-strain modulus over minimum-strain modulus suggests intracycle strain hardening in the critical gel. On the other hand, a higher value of viscosity at the minimum strain rate than at the highest strain rate indicates intracycle shear thinning in the system. We believe that the extension of the network of clay particles offers more resistance to elastic deformation that manifests in intracycle strain hardening. Furthermore, as the network gets strained, the particles get aligned in the direction of shear, which results in intracycle shear thinning. In order to model the linear–nonlinear transition, we perform a stress relaxation experiment on the critical gel state at different magnitudes of strains. Up to a certain magnitude of strain (γ ≤ 1), the vertical shifting of the independent stress relaxation curves leads to a superposed master curve. The corresponding strain-dependent vertical shift factor, also known as a damping function, is observed to decay in a quadratic fashion with strain beyond the linear region. Incorporating the damping function into the time strain separable K-BKZ integral equation leads to a quasilinear constitutive equation. We solve the TSS K-BKZ model on the application of an oscillatory strain using the spectral method. Although the proposed model predicts the intracycle stress response quite well in the linear as well as the nonlinear domain with model parameter a = 2.9, a slight modification in the damping function parameter (a = 2.1) leads to enhanced accuracy in the prediction of the experimental data. This leads to a simple formulation to analyze the nonlinear response of a critical gel.

See the supplementary material for the analysis of nonlinear behavior in terms of Chebyshev coefficients.

K.S. and Y.M.J. acknowledge the financial support provided by the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India. S.S. acknowledges the support by National Science Foundation under Grant No. NSF DMR-1727870. The authors were grateful to Professor Randy H. Ewoldt for providing us with the MITlaos software package. The authors were thankful to the anonymous reviewers for their insightful suggestions that have helped us improve the manuscript.

The authors have no conflicts to disclose.

Khushboo Suman: Data curation (lead); Formal analysis (equal); Investigation (lead); Methodology (lead); Writing – original draft (lead); Writing – review & editing (equal). Sachin Shanbhag: Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Yogesh M. Joshi: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

We used a Bayesian approach to infer the plausible range of the parameter a for each experimental setting of frequency and strain amplitude in the oscillatory shear experiments. We assumed a uniform prior distribution for the parameter, π(a) = U[0, 8], and defined a least-squares loss function,

where $σ̂i$ are the Fourier coefficients of the periodic shear stress response, N_f is the number of significant Fourier modes in the experimental data, and σ_d = 0.05 max(σ^exp(ω, γ₀)) assumes a 5% relative noise level. The superscripts “exp” and “model” stand for experimental data D, and the results of TSS-KBKZ model for a particular choice of a, respectively. The same loss function was used to fit a to the low-frequency oscillatory shear data in Fig. 7.

We assume that the likelihood function is a normal distribution,

According to Bayes theorem, the posterior distribution is given by

From the posterior distribution, we can find the 90% credible interval for a, which is shown in Fig. 10. We make the following observations:

• At small strain amplitudes, the posterior and prior distributions are nearly identical to the prior. This is expected since strains experienced by the sample in this (approximately linear) regime are too small to reveal any information about the nonlinear response.

• Figure 11 shows the corresponding Lissajous plots. The agreement is quite satisfactory in general. Except perhaps at large strain amplitudes, experimental data lie within the family of curves corresponding to the 90% credible interval for a. At large strain amplitudes, the predicted curves bulge beyond the experimental curves along the major and minor axes.

• From Fig. 10, a ≈ 2.1, shown by the dashed horizontal line, is found to be reasonable choice for most datasets explored. This observation explains the effectiveness of a = 2.1 ± 0.3 inferred from fitting the low-frequency data in Fig. 7. Interestingly, the credible interval of the damping function parameter is also consistent with a ≈ 2.9 obtained from step-strain experiments, which is shown by a dotted line on the same figure.

• At the highest frequency, which also corresponds to the highest shear rates at a given strain amplitude, there is a small but discernible preference for larger values of a as demonstrated by data corresponding to 31.4 rad/s in Fig. 10. This corresponds to stronger damping with higher shear rates.

We explore the residual error in the model and experimental Lissajous-Bowditch curves at large strain amplitudes in Fig. 11, and the preference for stronger damping at higher strain-rates in terms of changes in the relaxation spectrum within the framework of Yamamoto.⁸⁵ Rathinaraj et al.⁸³ used a similar system and found that the relaxation spectrum gets truncated as shear rate increases. Do we observe a similar trend in our analysis?

The fast spectral method used in this work exploits time-strain separability. In particular, it assumes that the memory function (or relaxation spectrum) in the K-BKZ model is independent of strain. Thus, we cannot use it to directly infer a strain-rate dependent memory function or spectrum. However, we can still ask whether our observations are consistent with a truncated relaxation spectrum at higher strain rates.

We demonstrate this with an illustrative example. Recall that the memory function of the critical gel m(s) = Sn/s¹⁺ⁿ. This is shown by the dashed line in Fig. 12(a). We synthetically truncated m(s) at s* = 10 s to obtain a truncated memory function shown by the solid red line,

The Lissajous curves corresponding to the full and truncated memory functions with a = 2.1 are shown in Fig. 12(b) at γ₀ = 2 and ω = 1 rad/s. Qualitatively, the shift in the curves due to truncation can be partially captured by using a stronger damping function with a = 4, as shown by the dotted line. The fit between the blue dotted line and the solid red line is not perfect; nevertheless, we note that the pattern of discrepancy is similar to that observed in Fig. 10. The approximately ellipsoidal curve, predicted assuming that the relaxation spectrum is unchanged, bulges beyond the target curve along the major and minor axes. Thus, the preference for larger a at higher strain amplitudes is consistent with a strain-rate dependent relaxation spectrum wherein slower modes are suppressed.