Viscoelastic phenomena in methylcellulose aqueous systems : Application of fractional calculus

Fractional calculus models can potentially describe the viscoelastic phenomena in soft solids. Nevertheless, their successful application is limited. This paper explored the potential of using fractional calculus models to describe the viscoelastic properties of soft solids, focusing on methylcellulose aqueous systems. Methylcellulose is an important food additive


Introduction
Rheology is a critical field in material science, which allows for describing the time-dependent properties of polymers.Viscoelastic phenomenological models constitute a key approach to achieve this goal (Fleischhauer et al., 2012;Peterson & Cates, 2021;Reyes-Melo et al., 2008).In the theory of linear viscoelasticity, such models usually rely on the separation of elastic and viscous properties of materials into two rheological elements, namely, springs and dashpots (Morrison, 2001).A spring is equivalent to an elastic solid ruled by Hooke's law, and a dashpot is comparable to a viscous liquid described by Newton's law of viscosity (Mezger, 2020).However, representing the rheology of polymers with only springs and dashpots neglects the existence of intermediate viscoelastic behaviors (Heymans & Bauwens, 1994;Jaishankar & Mckinley, 2013;Reyes-Melo et al., 2008).This issue is particularly common in materials with soft solid consistency, such as gels, since their rheological behavior lies in between a spring and a dashpot (Faber et al., 2017b;Warlus & Ponton, 2009).To address this problem, constitutive fractional calculus models can be utilized to define the stress-strain relationship of such complex materials.Previous studies have shown that models based on fractional calculus provide a viable solution to understand the rheology of complex materials (Bagley & Torvik, 1983;Faber et al., 2017b;Puente-Córdova et al., 2018;Reyes-Melo et al., 2008;Wagner et al., 2017).
Fractional calculus is a mathematical technique that can describe phenomena through derivatives and integrals of fractional or nonarbitrary order (e.g., calculating the 0.5th derivative of an x variable).In rheology of polymers, fractional calculus has practical applications in the viscoelastic phenomenological models of Maxwell, Voigt-Kelvin, and Zener.In these modified models, a "springpot" (also known as the Scott-Blair element) replaces the dashpot constitutive equation.The springpot, as depicted in Fig. 1, couples a Hookean spring and Newtonian dashpot into a single rheological element, interpolating from one to another to describe the viscoelasticity of a polymer (Koeller, 1984).Using Euler's notation, Eq. (1) defines a springpot as an α (0 ≤ α ≤ 1) fractional order derivative (D) of strain (γ) respective to time (t), σ = Gτ α D α t γ. (1) In Eq. (1), σ is the stress, G is the modulus, and τ is the relaxation time required for the polymer chains to reorganize into a different equilibrium state.As Fig. 1 shows, Eq. ( 1) generalizes Hooke's law (σ = Gγ) when α = 0 (i.e., spring element) and Newton's law (σ = ηD 1 t γ) when α = 1 (i.e., dashpot) (Faber et al., 2017a;Puente-Córdova et al., 2018;Reyes-Melo et al., 2008).This fractional order captures intermediate responses to the applied strain and temperature that are neither fully elastic nor viscous.
Given the power-law responses in viscoelastic materials, fractional frameworks, such as the one in Fig. 1, have the potential to describe intermediate behaviors between solid and liquid.For example, with a few meaningful parameters, fractional models have quantified material properties of gels, glassy polymers, and food systems (Faber et al., 2017b;Fang et al., 2015;Gao et al., 2022;Puente-Córdova et al., 2018;Wagner et al., 2017).Nevertheless, only a few studies have successfully translated fractional models into the rheological study of polysaccharides (Chen et al., 2013;Jóźwiak et al., 2015;Orczykowska & Dziubiński, 2014).Particularly, fractional phenomenological models can integrate understanding the impact of different non-covalent interactions in systems (e.g., food systems) where polysaccharides are rheology modifiers (Alamprese & Mariotti, 2015;Ginzburg et al., 2016;Moreira et al., 2017).Therefore, the present study models the rheological data of methylcellulose, a polymer widely used as a food additive (Hedayati et al., 2022;Sanz et al., 2005;Tanti et al., 2016).It further shows the potential of fractional calculus in describing the viscoelastic phenomena in methylcellulose.
Methylcellulose (MC) is a cellulose derivative with complex rheological behavior characterized by a thermogelation process (Coughlin et al., 2021;Miranda-Valdez et al., 2023;Reichler et al., 2021).When diluted in water and further heated, MC experiences a gel transition from a viscoelastic liquid to a viscoelastic solid (Coughlin et al., 2021;Mcallister et al., 2015aMcallister et al., , 2015b;;Schmidt et al., 2020).The relevant part for the present manuscript about the thermogelation of MC is its similarity to the glass transition in polymers (Rentería-Baltiérrez et al., 2020;Reyes-Melo et al., 2021).The thermogelation of MC observed during dynamic mechanical thermal analysis exhibits an antisymmetric shape in comparison to that of glass transition.Given the success of fractional frameworks in modeling primary relaxation phenomena in synthetic viscoelastic materials, it prompts the question of whether a fractional model can similarly describe the thermogelation process of MC.
The rheological behavior of MC systems has been extensively characterized in the literature using various approaches, including empirical equations, statistical design of experiments, molecular simulations, and machine learning (Alamprese & Mariotti, 2015;Desbrières et al., 2000;Ginzburg et al., 2016;Miranda-Valdez et al., 2022).Despite these efforts, a phenomenological model based on fractional calculus has not been applied to model the rheological properties of MC to date.This Greek letters α fractional derivative order (0 ≤ α ≤ 1) η* I.Y.Miranda-Valdez et al. in the linear viscoelastic regime.The FZM in this work is a fractional version of the classic Zener model extended to include a fractional derivative term that captures power-law responses in viscoelastic materials.This manuscript further shows that traditional frameworks fail to model the rheology of MC, whereas fractional rheology successfully captures responses to the applied strain and temperature.
Here, our research objective is to describe the mechanical and thermomechanical spectra of MC using fractional calculus.Furthermore, we aim to provide physical insights into the viscoelastic gel transition of MC and its complex behavior near the liquid-solid transition.We use a springpot (Scott-Blair model) and FZM with two springpots to describe the response of MC in frequency and temperature, respectively.The Supplementary Information provides the reader with a brief background about fractional rheology.Our hypothesis is that the fractional models will allow us to characterize the linear viscoelastic behavior of MC aqueous systems.To prove our hypothesis, we prepared MC solutions with polymer concentrations (c) ranging from 0.8 to 4.0 wt.%.Then, we assessed the rheological properties through oscillatory and rotational experiments; subsequently, we discussed the results using viscoelasticity concepts.The experimental results were finally modeled using fractional rheology and compared with traditional phenomenological models to show the potential of fractional calculus in viscoelastic materials.

Materials and methods
We used a food grade methylcellulose, Benecel™ MX MC-50000 (Ashland Specialties Belgium).According to the provider, the name MC-50000 implies a η of ca.50 Pa s (2 wt.% solution); more details are unavailable from the provider.This MC was used to prepare suspensions with polymer concentrations ranging from 0.8 to 4.0 wt.%.The MC has a density of 1322.7 kg m − 3 , which was determined using a gas pycnometer (Ultrapyc 1200e, Quantachrome, USA) after averaging 15 measurements (standard deviation of 0.3 kg m − 3 ).With the measured density of MC, we calculated its intrinsic viscosity [η] at 25 • C from a single-point viscosity measurement using a 0.1 wt.% aqueous solution; this was the minimum measurable concentration with a Couette geometry CC27.[η] was approximated using the Solomon-Ciuta equation given next in Eq. (2) (Solomon & Ciuta, 1962;Solomon & Gotesman, 1967), where η sp is the specific viscosity, η r is the relative viscosity, c is the polymer concentration in g ml − 1 , η is the dynamic viscosity of the solution taken from the zero shear region, and η 0 is the viscosity of the solvent (i.e., water).Following this approach, [η] was approximated to 1351.0 ml g − 1 .From the reciprocal of [η], we estimated the critical overlap concentration (c*), which indicates the transition from a dilute to a semidilute solution.For the studied MC, c* ≈ 0.0007 g ml − 1 or 0.0740 wt.%.Other authors have reported similar values of [η] and c* for an MC with comparable high molecular weight, degree of substitution, and nominal η (Arvidson et al., 2013).
We obtained the suspensions by dispersing MC in deionized ultrapure Milli-Q water at 50 • C with (resistivity of 18.2 MΩ cm) and then tempering the suspension at 4 • C overnight.According to the manufacturer, their MC satisfies the food additive requirements in standard E461 (JECFA, 2006).However, to achieve a thorough characterization of the polymer, we used Size Exclusion Chromatography (SEC), Raman spectroscopy, and Nuclear Magnetic Resonance (NMR) spectroscopy to estimate its weight average molecular weight (M w ) and chemical features.

Size exclusion chromatography
SEC experiments were carried out in 0.1 M NaNO 3 with an Agilent 1260 Infinity II Multi-Detector GPC/SEC System, including lightscattering (two measurement angles: 15 • and 90 • ), refractive index, and differential viscometer (VISC) detectors.Three Waters 7.8 mm × 300 mm Ultrahydrogel columns (500 Å, 250 Å, and 120 Å) with a 6 mm × 40 mm Ultrahydrogel guard column were used for separation with a flow rate of 0.5 ml min − 1 .The injection volume was 100 μl.The detectors were calibrated using a narrow dispersity pullulan standard with a nominal M w of 110 kg mol − 1 .Agilent OpenLAB CDS ChemStation Edition software was used for instrument control, and Agilent GPC/SEC software for data collection and handling.A refractive index increment of 0.140 ml g − 1 was used for molar mass calculation.

Raman spectroscopy
Raman spectroscopy was used to identify the functional groups of MC.For comparison, we measured the spectrum of commercial microcrystalline cellulose, Avicel® PH-101 (particle size ~50 μm), which was purchased from Sigma-Aldrich (Japan).On a Raman microscope inVia™ Renishaw® Qontor (United Kingdom), a green laser (λ of 532 nm and 50% intensity) scanned three times the polymer samples using a 20 × objective with high confocality.The spectra were acquired in a range from 4000 to 100 cm − 1 (exposure time of 20 s per scan at room temperature).A baseline correction was then made using the intelligent fitting tool of the device software, WiRE™ 5.3.

Nuclear magnetic resonance spectroscopy
NMR characterization of the methylcellulose sample was performed on a Bruker NMR AV III 400 spectrometer ( 1 H, diffusion edited 1 H and HSQC experiments, refer to King et al. (2018)) and on a Bruker NMR AV NEO 600 spectrometer (quantitative 13 C experiment).Peak assignments and degree of substitution (DS) determination were performed according to Kono et al. (2017) and Kono (2018).The sample was prepared by dissolving 40 mg of the methylcellulose in 1960 mg of DMSO− d 6 at elevated temperatures (65-90 • C) under stirring and transferring an aliquot of the resulting 2 wt.% solution into a standard 5 mm NMR tube.Due to the high molecular weight of MC, see Fig. 2a, and the thus resulting high viscosity of the solution, a lower measurement concentration was required than reported in the original protocol.Furthermore, the experiments were carried out at an acquisition temperature of 65 • C (compared to 90 • C), which was the maximum recommended temperature of the probe head used for overnight measurements.

Rheometry
The rheological characterization was performed using a Modular Compact Rheometer 302 (Anton Paar, Austria).We carried out three types of rheological experiments: dynamic mechanical analysis (DMA), dynamic mechanical thermal analysis (DMTA), and rotational tests.DMA and DMTA are techniques (partially inspired by impedance spectroscopy) suitable to measure systems that gradually change from viscoelastic liquids to viscoelastic solids or vice versa (Scott Blair, 1969).All the experiments were triplicated to confirm the shape of the rheograms.For DMA, we chose a parallel-plate geometry composed of a serrated plate (PP25/P2) and a Peltier plate temperature device (P-PTD200).During the experiment, the serrated plate deformed the samples within the linear viscoelastic region (LVR) applying a strain amplitude (γ) of 0.5% (1 mm gap between plates).For DMTA, we used a serrated Couette geometry.The bob tool (CP17/P6) applied a dynamic γ within the LVR of 1%, and the mechanical stimulus oscillated with an angular frequency (ω) of 10 rad s − 1 .
During DMTA, the rheometer controlled the temperature using a C-PTD 180/A accessory.The heating of the sample went from 15 to 80 • C, and the cooling from 80 to 15 • C (1 • C min − 1 ).On the other hand, DMA consisted of deforming the sample within a ω range from 0.01 to 100 rad s − 1 at 25 • C. The parallel-plate geometry was chosen for DMA because of inertial restrictions at high angular frequencies for samples with MC concentrations of 2.3-4.0 wt.%.Alternatively, the Couette geometry was ideal for DMTA, as this assembly allows more mass to be sampled and minimizes sample dehydration problems pointed out by Niemczyk-Soczynska et al. ( 2022).Measuring the thermal behavior of systems with concentrations of 2.3-4.0 wt.% at higher temperatures than 25 • C with parallel-plates ends up drying the MC sample, whereas the Coutte geometry has a maximum torque that is exceeded while measuring MC samples with concentrations higher than 2.3 wt.%.
In both DMA and DMTA, the theory of linear viscoelasticity assumes that γ generates in the material a sinusoidal σ, which has an in-phase and out-of-phase component.These components, properly named the storage modulus (G′) and the loss modulus (G′′), together with the phase shift angle (δ), are the main experimental results of the oscillatory tests (Alcoutlabi & Martinez-Vega, 1999).Eq. ( 3) shows a definition of δ as a function of ω and temperature (T), where the inverse function, tan δ, represents the energy dissipation associated with molecular frictions.
In this work, we used the complex modulus (G*) to assess the viscoelastic behavior in MC systems.G* is a complex number representing the sum of elastic and viscous contributions.From Pythagoras' theorem, G* is defined as in Eq. ( 4), The ratio of G* to ω is known as the complex viscosity (η*) (Morrison, 2001).When the Cox− Merz rule is valid, η* is equivalent to η, as Eq. ( 5) shows (Cox & Merz, 1958): (5) To evaluate the validity of Eq. ( 5), η and σ were recorded from rotational tests as a function of the shear rate (γ) in a logarithmic range from 0.01 to 1000 s − 1 at 25 • C, as in Eq. ( 6): The duration of each experimental point obtained for η was adjusted according to γ− 1 .Initially, we set the acquisition time as the reciprocal of γ (that is, when γ = 0.01 s − 1 , the measurement time was 100 s).Then, we logarithmically reduced the time to 10 s once γ = 1.0 s − 1 .Similarly, for the next interval, the measurement time was reduced from 10 to 1 s when it reached γ = 10.0 s − 1 .For the final interval, the measuring time was constant (1 s).We replicated all the rheological experiments to confirm the shape of the curves.Furthermore, before starting any experiment, the samples rested for 50 (DMTA) or 15 min (DMA and rotational tests) in the measuring position to ensure the thermal equilibrium of the sample.The difference in time is due to the sample mass required for each experiment, DTMA being the one sampling more mass.
It is worth mentioning that the samples were enclosed with a hood during every experiment.

Modeling
In this work, we used the springpot constitutive equation and fractional Zener model (FZM) in Fig. 1 and Fig. S3 (Puente- Córdova et al., 2018;Reyes-Melo et al., 2004).To describe the rheological mechanical spectra of the MC systems, we used the springpot in Fig. 1, which describes the complex rheological response of a viscoelastic material based on Eq. ( 7) (West et al., 2003, pp. 235-270), Alternatively, for the thermomechanical spectra, an FZM with two springpots is preferred because the version with one springpot fails to describe the complex viscoelastic response.Its solution in terms of G* is given by Eq. ( 8), The solutions to both models result from transforming the constitutive equations of the fractional frameworks to the frequency domain ω; this is best described in the Supplementary Information, and the literature cited next (Heymans & Bauwens, 1994;Puente-Córdova et al., 2018;Reyes-Melo et al., 2004;Schiessel & Blumen, 1993).In Eq. ( 8), the terms τ α and τ β corresponds to the relaxation times of the springpot elements α and β, respectively (Fig. S3).The exponents α and β are the fractional orders of each springpot, where 1 represents a dashpot and 0 a spring.Furthermore, the parameter G U is the modulus at low temperatures.In contrast, G 0 is the modulus at high temperatures, and i represents the imaginary unit.To describe the thermomechanical spectra of the MC systems (i.e., G* as a function of T), we defined a relationship between τ and T that depends on cooperative and non-cooperative molecular movements (Matsuoka, 1992(Matsuoka, , 1997;;Reyes-Melo et al., 2004).Refer to the Supplementary Information for more details.
The springpot results of the mechanical spectra were compared to the Maxwell Generalized Model (MGM), which requires fitting a Prony series to the results.For the MGM, we used the fitting tools developed by  Springer (2022).Eq. ( 9) and Eq. ( 10) show the representative equations of G′ and G′′, respectively.For the MGM model, in Eq. ( 9), and Eq. ( 10), G ∞ and G n are the equilibrium and relaxation modulus, respectively.n represents the circuit number in the Prony series, and m is the total of circuits in the Prony series.τ n is the relaxation time of each circuit.

Chemical and structural characterization
The results of the SEC analysis in Fig. 2a showed that MC had a wide bimodal molecular weight distribution.The molecular weight of the highest peak (M p ) was 360 kg mol − 1 , and M w was calculated to be 534 kg mol − 1 .Fig. 2b presents the Raman spectrum of MC compared to microcrystalline cellulose (MCC).Representative CH 3 , CH 2 , and CH vibrations of MC were in the spectral range from 3000 to 2700 cm − 1 .The qualitative analysis obtained by 1 H and HSQC spectra for MC agreed with the literature.A diffusion edited 1 H experiment revealed that MC had no significant low molecular weight impurities.

Rheological properties: mechanical spectra
The results of DMA in Fig. 3 indicated in general that the MC systems (2.3-4.0 wt.%) behave as entangled (physically) polymer dispersions.From Fig. 3a and b, G*, G′, and G′′ demonstrated being highly dependent on ω; the viscoelastic components increased as a function of this variable.Fig. 3b provides better insight into the viscoelastic behavior of MC.The MC systems studied satisfied the precondition for physical dispersion stability (G′ > G′′) when ω was higher than 0.01 rad s − 1 (Mezger, 2020).From 2.3 to 4.0 wt.%, the MC dispersions took the appearance of viscoelastic solids due to a combination of different phenomena giving rise to intermolecular forces between the OH and CH 3 groups located in the reducing ends of MC repeating units (Coughlin et al., 2021).Interestingly, from Fig. 3b, the MC systems already exhibited a rheological behavior typical of structured liquids at 25 • C (G′ > G′′); indeed, the first measuring point at ω = 0.01 rad s − 1 showed a cross-point between G′ and G′′.Furthermore, the entanglement of the polymer chain network increased as a function of ω since G′ and G′′ were of higher magnitude at high frequencies and the difference between G′ and G′′ was small.Such rheological behavior is typical of MC aqueous systems and gelling systems with imperfect networks experiencing a gradual relaxation of dangling structures (Arvidson et al., 2013;Moreira et al., 2017;Rubinstein & Colby, 2003).
As expected, the viscoelastic components and viscosity increased with the MC concentration (Moreira et al., 2017).The results of η* and η in Fig. 3c were in good harmony with the rheological behavior observed in Fig. 3a and b.That is, the viscosity measurements as a function of the shear rate showed a small, almost negligible plateau at low deformation rates.The absence of a plateau in the low shear rate region is correlated to the gel-like superstructure of the MC aqueous systems, predominating almost all over the measured ω and γ spectra.This is characteristic of partially cross-linked systems (Winter & Mours, 1999).The Cox− Merz rule was verified to be valid for tested MC concentrations since the experimental data of η* superposition the results of η.Due to the defined criteria G′ > G′′, the shape of the rheograms in Fig. 3c implied that the MC systems tended to behave as gels while at rest, and they had a yield point (Mezger, 2020).However, the MC systems were not true gels as the moduli (G*, G′, and G′′) showed frequency dependency (Funami et al., 2007).Moreover, the shear-thinning behavior observed in Fig. 3c is poorly understood for viscoelastic materials near liquid-solid transitions (Winter & Mours, 1999).
Regarding the results of tan δ, for all the MC concentrations, Fig. 3d shows that the systems dissipated more energy at low frequencies.This initial energy dissipation and further decay indicated that the MC systems were partially entangled.The MC systems experienced molecular frictions that contributed to increasing tan δ.Then, Fig. 3d highlights a molecular relaxation as a broad peak, extending from 0.1 to 100 rad s − 1 , in the case of the 4.0 wt.% sample.These broad peaks in tan δ denoted that more molecular frictions occurred in the MC systems, which may be attributed to the progressive entanglement of polymer chains and a possible broad distribution of relaxation times.The behavior of tan δ in Fig. 3d is typical of viscoelastic materials near liquid-solid transitions, since tan δ was frequency dependent (Winter & Mours, 1999).The latter confirmed that MC aqueous systems with high polymer concentrations (2.3-4.0 wt.%) were not true gels and could be rather described as "structured liquids".This could be attributed to the ability of MC systems to percolate still at higher temperatures than 25 • C.

Rheological properties: thermomechanical spectra
The results in Fig. 4 showed thermal hysteresis loops of the viscoelastic components of the MC systems.Here, the experiments were carried out at a single angular frequency and heating/cooling ramp, but it is worth mentioning that the thermomechanical spectrum depends on such experimental conditions (Arvidson et al., 2013;Coughlin et al., 2021;Mcallister, Schmidt, et al., 2015).It is also important to note that the thermomechanical spectra of samples with concentrations ranging from 2.3 to 4.0 wt.% were not measured due to physical limitations with the Couette geometry.As expected, Fig. 4a and b indicate that G* and G′ increased as a function of the MC concentration.Alternatively, as a function of T, G* and G′ initially decreased, then increased drastically until they reached a near-plateau.
The origin of the initial decay of the moduli in Fig. 4a and b is unknown, but it could be related to the activation energy of the thermogelation process.Alternatively, the sharp increase of the moduli is thought to occur when MC chains self-assemble into (toroidal-like) structures that percolate and collapse, forming a semicrystalline fibril network, a process known as fibril formation (Arvidson et al., 2013;Coughlin et al., 2021;Lott et al., 2013aLott et al., , 2013b;;Ginzburg et al., 2016).Applying the approach described by Miranda-Valdez et al. (2022) to calculate the gelation temperature (T ls ) from δ (the method uses δ as a parameter to define a viscoelastic solid and liquid), the gelation of the MC systems occurred at ~50 • C (during heating).This T ls remained almost unchanged as a function of c.Typically, T ls decreases as a function of c due to the increasing number of cross-links in the entangled gel system (Arvidson et al., 2013).However, because of the high molecular weight of the studied MC and the high MC content in the systems, the effect of the polymer concentration was less remarkable.
After reaching 80 • C, the moduli in Fig. 4a and b decayed as a function of T since the MC systems transition back from viscoelastic solids to viscoelastic liquids (gel-to-sol), showing that MC gelation is a reversible process.Via Cryo-TEM studies, the gel-to-sol process has been ascribed to a process known as fibril dissolution (Arvidson et al., 2013;Lott et al., 2013aLott et al., , 2013b)).The gel-to-sol temperature (T sl ) for the cooling cycle was also estimated from δ.The MC systems showed a T sl of ~38 • C. To calculate T ls and T sl , we preferred to use the method described by Miranda-Valdez et al. (2022) since no clear cross-point was observed between G′ and G′′.Moreover, the viscosity and storage modulus divergence method has been suggested to delay the real gel point due to shear-thinning constraints in a pseudoplastic polymer and the long relaxation times associated with MC thermogelation (Desbrières et al., 2000;Nelson et al., 2022;Winter & Mours, 1999).
Interestingly, the hysteresis loops of G′′ in Fig. 4c resembled the hysteresis responses observed in piezoelectric polymers; they were reminiscent of butterfly curves (Hafner et al., 2019).The asymmetric shape of the butterfly plot in Fig. 4c (before and after the heating-cooling cross-point) indicated that the intermolecular associations giving rise to the gel transition differed in the heating and cooling directions.The hysteresis butterfly-like loop of G′′ agreed with the formation and further dissolution of a metastable phase (fibrils).In the future, such phenomena could be understood by combining DMTA with in situ Raman spectroscopy (Fleissner et al., 2019).We suggest future works to evaluate the absorbance bands of MC functional groups and compare them when an MC system is heated and then cooled.For example, the curves of G′′ contained evidence of molecular relaxations at ~60 • C (heating) and ~30 • C (cooling).These two relaxations were more evident at lower MC concentrations, as the peaks in Fig. S5 exemplify.
The temperatures at which the relaxation peaks appeared are close to the identified T ls and T sl , providing evidence that such relaxation phenomena might be correlated with the gelation and gel-to-sol processes.Another feature to highlight from Fig. 4c is the position of the crosspoint between the heating and cooling curves.The cross-point shifted to higher temperatures as c was higher in the MC system; nevertheless, more studies are needed to investigate the reason for this phenomenon.Owing to the butterfly plot in G′′ as a function of T, the gelation of MC could be described as a non-linear process with elastic instability (Drinčić et al., 2011).However, from our observations, it is worth mentioning that the shape of G′′ depicted no butterfly at higher and lower concentrations than the ones explored with DMTA here, as in Fig. S5.The former agrees with the report of several authors about the different types of gels that MC forms depending on the polymer concentration (Arvidson et al., 2013;Chevillard & Axelos, 1997).
Concerning the results of tan δ, for all the MC concentrations, Fig. 4d shows that energy dissipation followed a damping plateau rather than a damping peak.The damping plateau changed to an abrupt decay as the MC systems approached their gelation temperature.Subsequently, at high T, tan δ formed another plateau of lower magnitude.The inverse process occurred during cooling.The shape of tan δ was in good harmony with the fact that the MC systems are phase change systems with a transition from viscoelastic liquids to viscoelastic solids and vice versa.At low temperatures, MC behaved as a viscoelastic liquid with damping properties similar to a dashpot element.On the contrary, at high temperatures, MC turned into a viscoelastic solid with a reduced damping effect.

Application of fractional calculus: mechanical spectra
We fitted the springpot in Fig. 1 to the DMA results of G* as a function of ω.We chose to use G* for the analysis since the principle of fractional calculus is to characterize viscoelastic behaviors that are not fully elastic nor viscous.Fig. 5a compares the springpot requiring only three parameters against the MGM with a Prony series of 3 and 10 circuits, equivalent to 7 and 21 parameters, respectively.It is unmistakable that the fractional model offered a more proper fit than the MGM by using only the three parameters in Table 1, which have physical meaning for the studied MC systems.The former does not imply that MGM should be replaced.Nevertheless, in complex systems, such as MC, when the traditional models fail to describe the rheological behavior, fractional calculus is a powerful tool to get more insight into viscoelastic properties (Bonfanti et al., 2020).Fig. 5b shows that the springpot adequately fitted all the experimental data.Compared to other modeling techniques, such as the statistical design of experiments, the springpot takes more relevance as it describes polymer systems in terms of viscoelasticity.However, the springpot model deviates from the experimental data at low frequencies.This deviation occurs since the springpot depicts one power-law regime.In the mechanical spectra from Fig. 3b, it is possible to see that the separations between the material function G′ and G′′ are smaller in the range from 0.01 to 0.1 rad s − 1 than in the rest of the spectra; in fact, there are cross-points between G′ and G′′ at 0.01 rad s − 1 .In the end, this behavior of G′ and G′′ at low frequencies is captured by the magnitude of G*.Therefore, at lower frequencies, G* starts depicting a power-law regime different from the one captured by one springpot.Two springpots in series may better fit the mechanical spectra, but here, we decided to use only one springpot to simplify the analysis.
From the springpot parameters in Table 1, G is the characteristic modulus of the sample, τ is the relaxation time, and α is the fractional order of the springpot element.The magnitude of the parameter G as a function of the concentration demonstrated that the MC systems developed a stiffer and elastic matrix when more polymer was present in the aqueous system.The latter was supported by the degree of viscoelasticity α decreasing to values approaching 0, indicating that the MC systems behaved more as viscoelastic solids when the concentration was higher (i.e., c = 4.0 wt.%).The magnitude of α is related to the degree of viscoelasticity of the samples.α initially had a value of 0.36 for the MC system with c = 2.3 wt.%.This indicated that the sample was a viscoelastic solid because α is closer to 0 than 1.Furthermore, when the MC concentration increased to c = 4.0 wt.%, α decreased to 0.24, indicating a gradual increment of the elastic nature of the MC systems.The elastic nature can also be observed from the material functions G′ and G′′ as a function of ω (Fig. 3b).For the lowest concentration (c = 2.3 wt.%), these two functions increased readily with the frequency, while for the higher concentrations, the increment was moderate.The moderate increment of G′ and G′′ as a function of ω for the MC systems with higher concentration speaks for more structured and elastic systems.The  results of α also agreed with the decreasing energy dissipation as a function of ω described for Fig. 3d (the material is more elastic), which was also in harmony with the traditional G′ > G′′ principle to identify if a material exhibits a dominant elastic behavior.Furthermore, the values of τ showed that the response time related to the transient behavior of the MC systems near their liquid-solid transition shortened with increasing the MC concentration.The behavior of τ confirms that with higher MC concentrations, the systems are more elastic and therefore require less time to relax.The classical explanation of the results given for Fig. 3 agreed with the viscoelastic behavior described by the springpot.There was a dominating elastic behavior, but since the MC systems were near their liquid-solid transition, they showed a rheological deviation from true gels.

Application of fraction calculus: thermomechanical spectra
For the thermal behavior, we used an FZM with two springpots since an FZM with one springpot failed to describe the complex rheological behavior observed experimentally with DMTA for MC.This was expected because the gelation of MC involves cooperative molecular mobility during the fibril formation and dissolution processes.In cooperative phenomena (e.g., glass transition), two springpots better described the viscoelastic properties of a polymer system (Reyes-Melo et al., 2004).To characterize the gelation and gel-to-sol, the first springpot α captured the viscoelastic behavior at low temperatures, while the springpot β did for the behavior at high temperatures.Including two springpots gave two characteristic relaxation times for the model; τ α and τ β are dependent on the temperature and follow the relationship exemplified in Eq.S31.Both relaxation times contain their respective pre-exponential factor; τ 0α and τ 0β associated with vibrations of atoms due to their order of magnitude close to 10 − 13 s, as Table 2 shows (Ngai, 1998;Reyes-Melo et al., 2016).Accordingly, Fig. 6 shows the FZM fit for the MC system with 1.8 wt.% concentration.The model worked well with the experimental data but failed to describe the initial decay of G* during heating.Such a phenomenon is difficult to model without including more elements in the fractional framework.Furthermore, we must remark that the FZM fitted the heating and cooling stages separately.There is a model for gelation and another for gel-to-sol.Future works should attempt to unify these two models into one that describes the hysteresis loop.For example, Drinčić et al. (2011) proposed a mass-spring-dashpot model to represent hysteresis loops shaped like a butterfly.Future research could follow a similar approach but use a springpot instead of a dashpot.
For DMTA describing gelation and gel-to-sol, the physical meaning of G U and G 0 is related to the viscous and elastic behavior of the sample.The values of G U and G 0 were consistent with the effect of concentration characterized for DMA using a springpot.In other words, G U and G 0 showed to be sensitive to the polymer concentration and increased their magnitudes proportionally (Table 2).The increasing magnitudes of G U and G 0 as a function of the MC concentration agreed with the observed experimentally in Fig. 4.
In addition to describing the degree of viscoelasticity, α and β could be correlated with molecular mobility during gelation and gel-to-sol phenomena.In the literature, it has been suggested to analyze two fractional parameters by taking the difference between them; a high difference between the parameters implies high molecular mobility (Puente-Córdova et al., 2018).For the concentrations from 0.8 to 1.8 wt.%, |α − β| during heating was consistent (~0.45), and higher than during cooling (~0.10).The consistent differences |α − β| indicated that molecular mobility was less dependent on polymer concentration.One can observe from Fig. 4d that tan δ (related to molecular frictions) depicts a similar hysteresis loop for the three MC concentrations, which agreed with the similar values of |α − β| for all concentrations.At the same time, the latter could agree with experimental findings on the independency of MC fibrillar diameter from polymer concentration, molecular weight, gel temperature, and polymer heterogeneities (Coughlin et al., 2021).However, this is difficult to confirm without cryo-TEM characterization or differential scanning calorimetry of the samples, which at the same time raises a puzzling question.How does the anisotropy induced during DMTA affect fibril formation and dissolution?
Regarding the values for |α − β|, during heating, they denoted higher mobility of the polymer chains than during cooling.The values of |α − β| agreed with the fact that in heating, the MC systems approached the transition from low viscosity, implying higher mobility for the fibril formation process than for the dissolution.While in the cooling stage, the materials left the high viscosity state, where mobility is naturally lower.In general, fractional modeling of methylcellulose systems allowed one to quantify viscoelastic properties in terms of a few material parameters.Further research may improve the prediction of methylcellulose rheology and optimize its ubiquitous use as a thickener, surfactant, or gelling agent.
In summary, in the introduction of this article, we stated a research question.Can fractional calculus describe the molecular mobility during thermogelation as it does for primary relaxation phenomena in polymers?The mechanical and thermomechanical spectra results confirmed that the fractional calculus models could potentially describe the

Table 2
Fractional Zener model (FZM) parameters for the results of dynamic mechanical thermal analysis (DMTA).G U and G 0 stand for low temperature modulus and high temperature modulus, respectively.τ 0α and τ 0β are the pre-exponential factors of the relaxation times of the springpot elements, while α and β are their fractional order.viscoelastic phenomena in MC, including thermogelation.However, establishing a physical meaning for the parameters is still challenging and a problem to solve in future research.We concluded that the thermogelation phenomena in MC involved high molecular mobility, which a few parameters could quantify.For the thermogelation, we did not compare the fractional framework against the Maxwell generalized model since the analysis becomes more complex and falls out of the scope of the article.Previous research suggested that the sol-gel transition in polymers follows a fractional behavior (Warlus & Ponton, 2009).

Sample Direction
Here, when MC systems were near their liquid-solid transition, they exhibited frequency and temperature-dependent responses that a springpot could effectively characterize while traditional frameworks failed in the task.The former is a reason to consider fractional approaches for modeling the viscoelastic phenomena of gel-forming materials.However, it is essential to remark that the fractional Zener model assumes that the measurement was in the steady state, as in this work.In addition, one could assume that if the fractional model fails to model the thermogelation behavior of MC, the molecular mobility follows a different behavior; when translating the fractional Zener model to study other gelling systems or the same MC in another concentration regime, this is key to highlight.Lastly, despite the potential use of the FZM, the model could only describe the heating and cooling branches separately.This work leaves an open question of whether it is possible to unify gelation and gel-to-sol in a single model for thermogelation hysteresis.

Statistical analysis
We finally describe the statistical analysis from the fractional rheology fittings to the rheological results.To do this, Table 3 summarizes the mean absolute error (MAE) of each fitting.The results showed that the fittings yielded a low MAE in frequency.The latter can be compared to the MAE resulting from fitting the MGM model to the DMA data of sample 2.8 wt.%.Fig. 7 shows that increasing the number of parameters in the MGM could reduce the MAE.Still, the error obtained with 201 parameters was four times bigger than the error of the springpot.Fig. 7 provides clear evidence of how the fractional model reduced the number of fitting parameters.On the other hand, the results in temperature showed a larger MAE because the model could not capture the initial decay during heating.The latter was also why the fitting for the cooling cycles had a smaller MAE than during heating.

Conclusions
The study showcased the effectiveness of using fractional rheology to describe various viscoelastic phenomena in materials like methylcellulose.While traditional models could not accurately describe the mechanical spectra of MC, the springpot and fractional Zener model characterized its viscoelastic properties with fewer parameters.Using fractional calculus demonstrated that it is a potential tool for assessing the complex nature of MC thermogelation and gel-to-sol transition.The model showed that thermogelation involves the cooperative mobility of polymer chains, and fractional calculus can model such behavior in the same manner as it does for primary relaxation phenomena in other polymers.This suggests that fractional calculus may become a potential tool for future research on complex viscoelastic systems.
study addresses this gap by demonstrating that a simple springpot can describe the rheological behavior of MC in frequency sweep tests, whereas a modified version of the Zener model, known as the fractional Zener model (FZM), can describe the thermorheological behavior of MC

Fig. 1 .
Fig. 1.Representation of the rheological elements and their constitutive equations.Hooke's law represents the spring element for elastic solids, where σ is the stress, G is the modulus, and γ is the strain.Alternatively, a dashpot represents Newton's law for viscous liquids, where η is the viscosity and D 1 t γ is the first derivative of strain over time.The springpot element combines both behaviors into a single constitutive equation, where τ is the relaxation time.When the fractional parameter α of the springpot constitutive equation is 0, it describes Hooke's law for elastic solids.On the other hand, when α is 1, the springpot expresses Newton's law for viscous liquids.

Fig. 2 .
Fig. 2. Structural characterization of methylcellulose (MC).a) Size exclusion chromatography (SEC) results show an approximation to the bimodal molecular weight distribution (MWD) of MC.From the chromatography data, the following molecular weight averages were calculated approximatively: M n = 353 kg mol − 1 ; M w = 534 kg mol − 1 ; M z = 756 kg mol − 1 ; M p = 360 kg mol − 1 .b) Raman spectra of MC and microcrystalline cellulose (MCC), where the spectrum of MC shows its representative CH 3 , CH 2 , and CH vibrations within the spectral range from 3000 to 2700 cm − 1 .The inset figure in b) depicts the possible chemical structure of MC or MCC."R" can be CH 3 or H for MC, while R is only H for MCC.The numbering indicates the carbon positions, and "n" is the number of repeating units.

Fig. 4 .
Fig. 4. Dynamic mechanical thermal analysis (DMTA) of methylcellulose aqueous systems with different polymer concentrations (0.8-1.8 wt.%) showing the thermal hysteresis loops of the viscoelastic components; experiments were carried out at a constant angular frequency (ω) of 10 rad s − 1 and constant heating rate of 1 • C min − 1 .a) Complex modulus (G*) as a function of the temperature (T).b) Storage modulus (G′) as a function of T. c) Loss modulus (G′′) as a function of T. d) tan δ as a function of T.

Fig. 5 .
Fig. 5. Fitting of the dynamic mechanical analysis (DMA) results at constant temperature T = 25 • C using the springpot model.a) Comparison between the experimental data of the MC system with 2.8 wt.% and the predictions of the springpot and Maxwell generalized model (MGM).The number in the parenthesis from the legend indicates the parameters required for each fitting.b) DMA experimental results and springpot fittings (dashed lines).

Fig. 6 .
Fig.6.Fitting of a dynamic mechanical thermal analysis (DMTA) result using the fractional Zener Model (FZM) with two springpots.The plot shows the complex modulus (G*) response as a function of temperature at a constant angular frequency (ω) of 10 rad s − 1 .The reader must notice that the heating and cooling ramps are modeled separately.

Fig. 7 .
Fig. 7. Mean absolute error (MAE) of fitting the dynamic mechanical analysis (DMA) of sample 2.8 wt.% using the generalized Maxwell model (MGM) and springpot model.Increasing the number of circuits in the MGM and, therefore, the number of parameters reduced the error.

Table 1
Parameters of the springpot used to fit the results of dynamic mechanical analysis (DMA).G stands for the modulus, τ is the relaxation time of the springpot element in Fig.1, while α its fractional order or degree of viscoelasticity.

Table 3
Mean absolute error (MAE) of the springpot and fractional Zener fitting for the mechanical and thermomechanical spectra results, respectively.The isothermal experiments represent dynamic mechanical analyses (DMA).The heating and cooling experiments correspond to dynamic mechanical thermal analyses (DMTA).