Phase transitions of cellulose nanocrystal suspensions from nonlinear oscillatory shear

Cellulose nanocrystals (CNCs) self-assemble in water suspensions into liquid crystalline assemblies. Here, we elucidate the microstructural changes associated with nonlinear deformations in (2–9 wt%) CNC suspensions through nonlinear rheological analysis, that was performed in parallel with coupled rheology—polarized light imaging. We show that nonlinear material parameters from Fourier-transform rheology and stress decomposition are sensitive to all CNC phases investigated, i.e. isotropic, biphasic and liquid crystalline. This is in contrast to steady shear and linear viscoelastic dynamic moduli where the three-region behavior and weak strain overshoot cannot distinguish between biphasic and liquid crystalline phases. Thus, the inter-cycle and intra-cycle nonlinear parameters investigated are a more sensitive approach to relate rheological measurements to CNC phase behavior.

source and extraction method (Dong et al. 1998) and are usually 4-10 nm in width and 100-1000 nm in length (Habibi et al. 2010;Klemm et al. 2011;Navarro et al. 2021). Their crystalline nature is the result of removing amorphous segments from the cellulose source of choice. In addition, using sulfuric acid for the extraction results in CNC particles with sulfate half-ester groups on their surface. This leads to electrostatically stabilized CNC aqueous suspensions (Marchessault et al. 1959;Revol et al. 1992;Dong et al. 1996). In static conditions (no flow), dilute aqueous suspensions of CNCs form an isotropic phase. As the CNC concentration is increased above a critical concentration, the suspensions assemble into a biphasic system consisting of isotropic and liquid crystalline ordered domains. The formation of ordered domains with increasing concentration has been described theoretically by Onsager (1949) who elaborated how suspensions of nanorods can spontaneously develop long-range orientational order as a result of the gain in translational entropy overcompensating the loss in rotational entropy. An additional increase in concentration will lead to a complete liquid crystalline phase ). The liquid crystalline (LC) domains comprise chiral nematic and/or nematic CNC structures (Revol et al. 1992;Abitbol and Cranston 2014). Chiral nematic crystals are characterized by a cholesteric pitch, a distance along the helical axis that separates nanorods of the same orientation after a 360° rotation (Habibi et al. 2010). In polarized light microscopy, chiral nematic orders in biphasic and liquid crystalline suspensions can be observed in the form of a 'fingerprint pattern', where the lines indicate the pitch of the chiral order (Revol et al. 1992). The nematic phase refers to a crystal orientational order having the CNCs' long axes approximately parallel, for distances significantly higher than the dimensions of the nanorods (Dufresne 2012). Further increasing the CNC concentration results into a (repulsive) glassy phase (Xu et al. 2018). Their hierarchical structure makes CNC suspensions, and films resulting therefrom, optically active materials which display birefringence. Birefringence can arise from any repeating pattern in a structure, chiral nematic and nematic being two of them (Ureña-Benavides et al. 2011;Dumanli et al. 2014). Subjecting CNC suspensions to flow can cause preferential orientation in the flow direction (Hubbe et al. 2017). The birefringence properties of CNC can be used in films as optical filters for color filters, smart clothes, and in solar-gain-regulating building technologies, free-standing iridescent photonic films with tunable structural colors to detect changes in the environment or in decorative coatings (Haniffa et al. 2016;Fernandes et al. 2017). The self-assembly of CNCs is known to depend on many factors, such as surface charge, ionic strength, aspect ratio, concentration, sonication, temperature, etc. (Dong et al. 1998;Shafeiei-Sabet et al. 2013;Abitbol et al. 2018;Derakhshandeh et al. 2013;Zhang et al. 2018).
Nonlinear rheological analysis, such as Fourier-Transform (FT) Rheology and stress decomposition, can increase measurement sensitivity and have the potential to reveal material response features that are not otherwise observable in linear viscoelastic measurements (Wilhelm and van Dusschoten 2001;Hyun et al. 2011). The nonlinear oscillatory shear analysis has gained increasing attention as a rheological testing method especially in polymer melts and solutions, filled polymeric systems, suspensions and emulsions (Hyun et al. 2011). However, the nonlinear viscoelastic rheological characterization of CNC suspensions is not well established. Chen et al. (Chen et al. 2017) investigated the behavior of CNC particles in PVA and CMC aqueous solutions with respect to linear viscoelastic dynamic moduli in the nonlinear region. They reported that the PVA has stronger interactions with negatively charged CNC particles, resulting in a lower percolation threshold of CNC in the PVA solution during Small Amplitude Oscillatory Shear (SAOS) and a weak strain overshoot behavior in Large Amplitude Oscillatory Shear (LAOS). However, we note that linear viscoelastic dynamic moduli do not accurately represent nonlinear material response and that nonlinear material parameters were not further evaluated. Van den Berg et al. (2018) studied the interfacial layers of CNC particles in (interfacial) SAOS, LAOS and large amplitude oscillatory dilatation (LAOD). Lissajous-Bowditch diagrams in LAOS exhibited strain stiffening, yielding, and flow of the interfacial layers for the high degree of chemical substitution. Another study on CNC in aqueous media with added NaCl revealed that the nonlinear viscoelastic measurements are significantly sensitive to detect changes in the microstructure and yield valuable information about the network structure (Abbasi Moud et al. 2020). Furthermore, performing the nonlinear oscillatory shear analysis enabled detecting the formation of a double network in CNC-PVA/salt hydrogels at high CNC loadings (Moud et al. 2021). Several studies indicate that some materials can exhibit unique nonlinear signatures with a non-quadratic scaling of I 3∕1 with the strain amplitude, I 3∕1 ∝ n 0 , n ≠ 2 , that could be related to microstructural or molecular material properties. Hyun et al. (2007) reported a lower slope of I 3∕1 for long-branched polymers compared to linear polymers in the MAOS region, which could be used as a measure of the degree of branching. There has been evidence that percolated nanocomposite polymer systems can exhibit peculiar nonlinear scaling behavior (Kádár et al. 2017;Gaska and Kádár 2019). Natalia et al. (2020) reported all evidence to date of noninteger power-law expansions for nonlinear viscoelastic properties. Building on previous evidence, Kádár et al. (2020b) determined that the onset of unique nonlinear signatures, dubbed 'nonlinear oddities', in the MAOS region matches the electrical percolation threshold of hierarchical reduced graphene oxide (HrGO)-polypropylene (PP) nanocomposites and subsequent percolated network consolidation. Meurer et al. (2021) applied the same analysis framework for chemical gels in the form of pressure sensitive adhesives. Changes in nonlinear signatures were recorded around the gel-point, however, both inter-cycle and intra-cycle nonlinear behavior differed considerably from those recorded in polymer nanocomposites. Associating such 'nonlinear oddities' to physical gels in the form of suspensions has yet to be elucidated.
Optical visualization techniques combined with rheological characterization can help understand the relationship between microstructure, phase transitions, and rheological properties . One of the commonly used rheological hyphenated methods for optically active suspensions is polarized light imaging (PLI). Shafiei-Sabet et al. (2012) investigated the effect of temperature and ultrasound energy using the rheo-PLI technique. These two factors affect the suspensions' microstructure depending on the concentration and the pitch of the chiral nematic domains. The degree of sulfation of CNC particles also significantly affects the phase transitions from isotropic through the liquid crystal, to a gel appearance (Shafeiei-Sabet et al. 2013). Hausmann et al. (2018) studied the dynamics of CNC alignment with high weight fractions (10-40 wt%) of particles under shear flow coupled with PLI. Kádár et al. (2020a) reported microphase transition sequences of CNC suspensions in shear flow through combined rheology and polarized light optical visualizations. Fazilati et al. (2021) combined rheology with polarized light imaging to elucidate the thixotropic behavior of CNCs suspended in water.
Although the CNC particles were first described in 1949 by Rånby et al. (1949), and numerous factors affecting crystalline phases in suspension identified, there is still a gap in a clear understanding of phase transitions and their relationship to the rheological behavior of CNC suspensions. It is common to use steady shear and (linear viscoelastic) oscillatory shear measurements to associate rheological behavior to CNC phases, however, a phase identification therewith remains a challenge. Our hypothesis is that nonlinear oscillatory shear analysis through increased measurement sensitivity and additional nonlinear rheological parameters is a more sensitive framework to elucidate CNC phases from rheology. To our knowledge, this is the first study to address this. To provide additional insight into the connection between the nonlinear material response of the suspensions and microstructural developments in the hierarchy of the suspensions we analyze the combined rheo-PLI optical visualizations in oscillatory strain sweep tests. While rheological measurements capture the bulk material response corresponding to all levels of the CNC hierarchy, macroscopic PLI data captures only the higher levels of the hierarchy, i.e. self-assembled tactoids/nematic/chiral nematic CNCs distributed in a polydomain structure within the measurement gap, and their orientation dynamics (Fazilati et al. 2021).

Materials
Cellulose nanocrystals (CNCs) in powder form were purchased from CelluForce, CELLUFORCE NCV100 -NASD90 (Montreal, Canada), having sulfur content of about 0.9% and sulfate content of about 250 μmol/g. According to X-ray photoelectron spectroscopy (XPS) analysis sulfur content was 0.39 ± 0.04 atom% that was performed using a PHI 5000 VersaProbe III Scanning XPS Microprobe at an angle of 45°. Aqueous CNC suspensions in the presence of 0.1 mM NaCl salt showed an average zeta potential value − 32.0 ± 3.4 mV, see Table S1 in the Supplementary information. Suspensions of 2-9 wt% CNC concentrations were prepared with deionized water (Millipore Milli-Q Purification System) and stirred on a bench shaker for 24 h. Prior to all rheological measurements, each suspension was treated for 15 min in an ultrasound bath (Bransonic Ultrasonic Bath 221, Connecticut, U.S.) operating at 50 kHz and 235 W.
CNCs are reported with length between 100 and 300 nm. The diameter is typically between 4 and 10 nm (Qing et al. 2013;Saito et al. 2009). We have determined the diameter of the CNCs to be 4.1 ± 1.0 nm (Navarro et al. 2021;Fazilati et al. 2021) and hence deem the aspect ratios of the CNCs to be between 25 and 75.

Optical microscopy
Transmission optical microscopy was carried out with a Carl Zeiss A1 (Oberkochen, Germany) optical microscope equipped with linear polarizers in a cross-polarizer setup. Each sample was placed and squeezed in the space between two microscope slides after 15 min sonication in an ultrasound bath.

Rheological characterization
All rheological measurements were performed on an Anton Paar MCR 702 TwinDrive rheometer (Graz, Austria) at 23 • C . Before measurement, each sample was allowed to relax for 240 s after moving to the gap position (1 mm for all tests). Two sets of configurations were used: (i) a single motor-transducer setup for rheo-PLI and (ii) a separate motor-transducer setup for nonlinear oscillatory shear analysis.

Rheo-PLI
Rheo-PLI experiments were performed in a single motor-transducer configuration using a custom rheooptical visualization setup based on the P-PTD200/ GL accessory with a parallel-plate geometry of (2R =)43 mm diameter. Two linear polarizers were placed at 45 • relative orientation above and below the parallel-plate setup, respectively. HD video recordings (1280×720 px) were performed at 60 fps from below, perpendicularly to the shear plane, using a Canon DSLR camera setup comprising a Canon L-series 100 mm macro lens combined with a Canon 25 mm extension tube (Tokyo, Japan). More details about the setup can be found elsewhere (Kádár et al. 2020a), while a schematic diagram of the setup can be found in Fig. S1 (Supplementary information). Space-time diagrams were constructed from the video recordings by extracting one line of pixels at a fixed position out of each still frame and appended to a new image having the x-axis corresponding to the experimental time and the y-axis corresponding to the length L (Kádár et al. 2020a), see also Fig. S1. Both oscillatory and steady shear measurements were performed on the setup. The strain sweep measurements were performed at 0.6 and 4 rad/s with the strain amplitude ranging from 0.01 to 1500%. In turn, the steady shear measurements were conducted within the shear rate range of 0.001 to 100 s −1 using a custom procedure for steady-state detection (Kádár et al. 2020a).

Nonlinear oscillatory shear analysis
In the separate motor-transducer configuration, linear and nonlinear oscillatory shear measurements were performed using a parallel-plate geometry of (2R =)50 mm in diameter. Strain sweep measurements were performed within a strain amplitude range of 0.01 to 1500 % at the following frequencies: 0.6, 1, 2, and 4 rad/s. Performing measurements at different angular frequencies could reveal essential information about the microstructure, e.g. through the dependency of I 3∕1 on (Kádár et al. 2020b). The nonlinear data analysis of the shear stress output signal was performed in the framework of Fourier-Transform Rheology and stress decomposition analysis. A graphical overview thereof can be found in Fig. 1. One simple way to distinguish between linear and nonlinear oscillatory shear is to consider the measured time dependent shear stress response and its corresponding Fourier spectra. In linear viscoelastic oscillatory shear, a sinusoidal shear strain input (t) = 0 sin( t) , where 0 is the strain amplitude and the angular frequency, will result in a sinusoidal shear stress output (t) = 0 sin( t + ) , shifted with the phase angle . Accordingly, the Fourier spectrum of the shear stress contains only the fundamental intensity I 1 of the imposed , see Fig. 1a. In contrast, for nonlinear deformations a sinusoidal strain input translates into a distorted shear stress output which leads to higher harmonics in the corresponding Fourier spectra, in addition to the fundamental intensity, I n , n = 3, 5... , (Hyun et al. 2011) see Fig. 1b. This can also be translated into the behavior of the linear viscoelastic dynamic moduli, the storage and loss modulus G ′ , G ′′ , in dynamic strain sweep tests, which is how the linear viscoelastic limit is usually quantified. For strain amplitudes where the shear stress output is sinusoidal the elastic and loss modulus are independent of the strain amplitude. However, as the applied strain amplitude is increased and the shear stress output becomes distorted, the storage and loss moduli become functions of the applied strain ampli- Fig. 1c. Thus, it should be noted that in this case the linear viscoelastic dynamic moduli contain an unreported intrinsic error owing to the distorted shear stress signals. This nonlinear component of the material response can be quantified using nonlinear material parameters. The third relative higher harmonic, I 3∕1 ≡ I 3 ∕I 1 , can be used to characterize the shear stress signal nonlinearities, as I 3 contains the largest nonlinear contribution to the Fourier spectra (Hyun et al. 2011). For small strain amplitudes (small amplitude oscillatory shear, SAOS) I 3∕1 is scattered with general scaling behavior of I 3∕1 ∝ −1 0 corresponding to nonlinearities due to instrumentation noise (Hyun et al. 2011). In this framework, the onset of the nonlinear behavior is marked by a steep increase in I 3∕1 . It can be shown through theoretical reasoning (Hyun and Kim 2011) that I 3∕1 ∝ 2 0 and this strain amplitude region is called medium amplitude oscillatory shear (MAOS). For even higher strain amplitudes, i.e. large amplitude oscillatory shear (LAOS), the quadratic scaling is lost (Fig. 1c). In addition to quantifying directly the nonlinearity of the shear stress signal, we briefly note that by taking as example the data presented in Fig. 1c the dynamic range of I 3∕1 spans approximately 3 decades compared to less than one decade in G ′ , G ′′ .
In addition to material parameters based on shear stress Fourier spectra, which essentially quantify how distorted the shear stress signal is, nonlinear material parameters can also be derived based directly on the intra-cycle distortion of the shear stress. The elastic ( (t)∕ 0 vs. (t)∕ 0 ) and viscous ( (t)∕ 0 vs. ̇(t)∕̇0) Lissajous-Bowditch (LB) diagrams can therefore be used to assess the nonlinear material response (Ewoldt et al. 2008). An elliptic shaped response indicates a linear viscoelastic behavior, see the light red and light blue where G ′ M is the minimum-strain and G ′ L is the largestrain modulus, and where ′ M is the minimum-rate and ′ L the large-rate viscosity, to quantify the intra-cycle material behavior based on LB diagrams. These parameters were used to define the strain-stiffening (S) and shear-thickening (T) dimensionless ratios (Ewoldt et al. 2008): where S > 0 indicates intra-cycle strain-stiffening, S < 0 intra-cycle strain-softening, T > 0 intra-cycle shear-thickening and T < 0 intra-cycle shear-thinning. (1)

Morphological characterization
Polarized optical microscopy (POM) was used to capture micrographs of the suspensions investigated to observe CNC phases independently from the rheological characterization, (Fig. 2). No birefringence is observed for 2 wt% CNC, which is therefore categorized as isotropic. As concentration increases, the micrographs show birefringence which is deemed to be due to the formation of liquid crystalline domains. The biphasic textures are visible for 3-5 wt% CNC suspensions, where a liquid crystal fraction coexists as bright droplets with a dark isotropic phase. After reaching the critical concentration from biphasic to liquid crystalline, there is no visible dark isotropic phase and they merge into one domain for 6-9 wt% CNC. Due to the absence of clear 'fingerprint' textures in the micrographs, we consider the evidence for whether we have nematic or chiral nematic structures in the CNC samples investigated, inconclusive. To summarize, using the given preparation conditions, the 2 wt% CNC suspension is isotropic, 3-5 wt% CNC are biphasic, and 6-9 wt% CNC are liquid crystalline. We note that biphasic compositions in the range of 3-6 wt% have been identified for aspect ratios of 61 ± 40 (Honorato-Rios et al. 2018).

Steady and oscillatory shear tests
Here we present steady shear and linear viscoelastic oscillatory shear tests, which are the most common types of tests employed for the characterization of CNC suspensions. The birefringence data from rheo-PLI steady shear tests is also discussed to underline the correlation between the birefringent properties of the suspensions under shear and their phase behavior. Rheo-PLI data from the oscillatory shear measurements is discussed after the nonlinear analysis. The steady shear viscosity functions of the CNC suspensions investigated are shown in Fig. 3. From a qualitative point of view, the steady shear viscosity of the suspensions could be grouped into two categories. In the first category, a zero-shear viscosity Newtonian plateau was apparent followed by a shear-thinning region. Such viscosity functions are characteristic of isotropic phase suspensions, see 2 wt% CNC. The second group corresponds to 3-9 wt% CNC, where the suspensions showed evidence of a three-region steady shear viscosity function, which is characteristic of liquid crystal systems ). Onogi and Asada (1980) were the first to assign the three-region behavior in polymer liquid crystals, however, there is overall quite a large discrepancy among subsequent studies on the estimation of CNC phases based on the three-region viscosity curve.  Haywood et al. (2017) where the threeregion behavior was observed in both biphasic and liquid crystalline phases. We note that Haywood et al. (2017) were able to distinguish more than on shear thinning region beyond what could be visually easily assessed by fitting the low and high shear rate regions with the Power-law and comparing the flow indices. Similar observations were made by Fazilati et al. (2021) for suspensions differing only slightly in the preparation steps compared to the present study. In the earliest estimation, Onogi and Asada (1980) hypothesized that at low shear rates in Region I, the polydomain chiral structures are not altered, and domains are able to flow over each other, which results in the first shear-thinning region. At intermediate shear rates, Region II is observed, sometimes called the Newtonian region, where the orientation of the initially unaffected polydomain starts to change. Haywood et al. (2017) identified through Rheo-SANS measurements that the anisotropy in this region remains relatively constant, with an isotropic ring still visible, which could indicate that only part of the rods align in the flow direction, with a significant part of the material remaining randomly oriented. The liquid crystal alignment occurs in Region III, which is a secondary shear-thinning region of the viscosity curve. We also note that for the higher concentrations, above 5 wt%, the relative change in the viscosity is comparatively small. This phenomenon could be explained by considering the Brownian motion of the particles and van der Waals interactions (Dufresne 2012). At a low CNC concentrations of aqueous suspensions, the rate of particle Brownian rotation is considerable and prevents alignment of the particles (Derakhshandeh et al. 2013). As concentration increases, CNC particles can (locally) orient and form nematic and chiral nematic structures. Similar observations were made by Pignon et al. (2021) who investigated CNC suspensions (2-12.2 wt%) using two combined techniques, rheo-SAXS and rheo-SALS.
Overall, from the present data, including also previous findings on very similar systems (Kádár et al. 2020a;Fazilati et al. 2021), assigning the three-region behavior to a particular phase remains ambiguous. While power law flow indices, see Eq. S1 in the Supplementary information, of low and high shear rates show a similar behavior with the results of Haywood et al. (2017), see Fig. S2 and Table S2 in the Supplementary information, the three-regions are visually apparent only between 5 and 9 wt% CNC concentrations and even therein not fully consistent, e.g. see 7-8 wt% in Fig. S2. However, a relatively small change in the viscosity for concentrations 6-9 wt%, compared to lower concentrations, could indicate the suspensions are in the LC phase. Phase transitions in CNCs can also be inferred based on the birefringence patterns of the suspensions in flow. Fig. 4 presents the evolution of CNC suspensions flow with shear rate as PLI images corresponding to the tests in Fig. 3. The isotropic 2 wt% CNC suspension did not indicate changes in color through all applied shear rates. However, at ̇= 23 s −1 , a Maltese-cross pattern was observed. This pattern indicates a parallel orientation of one of the main axes of the optical indicatrices of a birefringent material, which is CNC, in the plane of polarization of the incident light (Völker-Pop 2014). As concentration increases, the Maltese-cross appeared at lower shear rates. For the biphasic 3 wt% CNC, it Fig. 4 Still frames showing the birefringence patterns at selected shear rates from the steady shear measurements in Fig. 3 was observed at ̇= 5 s −1 , while for 4 wt% CNC at ̇= 0.3 s −1 . For concentrated biphasic (5 wt% CNC) and all LC concentrations, a Maltese cross could be observed at ̇= 0.1 s −1 . Additionally, as the fraction of the liquid crystal phase increased, a colorful pattern was observed both as large isochromatic areas at the beginning of the tests and in the Maltese-cross patterns. We note that while a color-less suspension prior to the beginning of the tests ( ̇= 0 ) followed by a colored Maltese-cross pattern generally signals a biphasic suspension and large isochromatic areas at ̇= 0 generally suggest liquid crystalline phase suspensions, the observed birefringence patters cannot be used to accurately determine CNC phase transitions with increasing concentration.
Dynamic strain sweep measurements, at two different angular frequencies, 0.6 and 4 rad/s, are presented in Fig. 5. The rheological behavior of isotropic 2 wt% CNC was dependent on the applied angular frequency. For = 0.6 rad∕s the suspension exhibited gel-like behavior, G � ≅ G �� , however, for = 4 rad/s, the same concentration exhibited liquid-like behavior, G ′′ > G ′ . This can indicate that at 2 wt% CNC the suspension contains a weak gel-like CNC network. All higher concentrations, 3-9 wt% CNC, exhibited gel-like behavior, G ′ > G ′′ , independently of the applied . The biphasic 5 wt% CNC suspension exhibited a weak strain overshoot (WSO) whereby G ′′ locally increases at the transition to the nonlinear viscoelastic regime. The WSO is further observed for all higher concentrations, > 5 wt% CNC. In general, a higher load of CNC particles would lead to stronger associations that can resist restructuring at large deformations causing the overshoot in G ′′ . We note that while the -dependent liquid/gel-like behavior is characteristic solely of the isotropic 2 wt% CNC, the onset of WSO behavior is well withing the biphasic concentration range and also characterizes the liquid crystalline concentrations.

CNC phase bahavior from nonlinear oscillatory shear analysis
The nonlinear material response as expressed by the third relative higher harmonic, I 3∕1 , from the FT analysis of dynamic strain sweep measurements is presented in Fig. 6. The instrumentation noise in the SAOS region is represented by semi-transparent datapoints. For the isotropic 2 wt% CNC suspension, Fig. 6a, the nonlinearities were detected at the highest recorded strain amplitudes, as expected considering that the suspension viscosity leads to limiting torques for the geometry used (Wojno et al. 2019). A strong dependence in the MAOS region was observed. For = 0.6, 1 rad/s two scaling regions, with I 3∕1 ∝ n 0 , where (i) n ≈ 1.8 for 0 ∈ [3, 10] % and (ii) n ≈ 0.5 for 0 ∈ [12, 60] % can be distinguished. A decrease in I 3∕1 magnitude was recorded with increasing > 1 rad/s. For = 2 rad/s two scaling regions remained identifiable: (i) n ≈ 2 for 0 ∈ [6, 23] %, (ii) n ≈ 0.85 for 0 ∈ [24, 95] %. However, for = 4 rad/s, the behavior approached one scaling region, where n ≈ 1.9 for 0 ∈ [15, 95] %, thus approaching the quadratic scaling theoretically expected. The tendency towards quadratic scaling with increasing , for nonlinear signatures where I 3∕1 = I 3∕1 ( ) , has also been observed in percolated polymer nanocomposites. This behavior has been attributed to a weakly percolated network being disrupted by shear and confirmed through electrical conductivity measurements (Kádár et al. 2020b). For the CNC samples investigated, the existence of a weak gel network is confirmed by the linear viscoelastic dynamic moduli. Significantly low CNC concentrations lead to lower interaction between ions in aqueous suspension, which can be destroyed Fig. 6 Third relative higher harmonic, I 3∕1 , from dynamic strain sweeps for = 0.6, 1, 2, 4 rad/s: a 2 wt%, b 3 wt%, c 4 wt%, d 5 wt%, e 6 wt%, f 7 wt%, g 8 wt%, h 9 wt% CNC while applying shear flow. We note here that such dependence seems to be inherent to the gel/percolation concentration, Kádár et al. (2020b), Meurer et al. (2021), which for other CNC systems may occur at higher concentrations than the isotropic-biphasic transition. In such a case, given its expected weakly elastic rheological behavior, a single quadratic scaling in the MAOS region would be expected for concentrations in the isotropic case, thus still being distinguishable from the other phases. For the biphasic 3-4 wt% CNC suspension, Fig. 6b-c, a very weak angular frequency dependence was observed and was practically absent at 5 wt% CNC, Fig. 6d. Two scaling regions could be inferred for the biphasic 3-5 wt%: for 3 wt% CNC (i) n ≈ 0.6 for 0 ∈ [1, 5] % and (ii) n ≈ 1.5 for 0 ∈ [5, 10] %; for 4 wt% CNC: (i) n ≈ 1 for 0 ∈ [0.7, 3.8] % and (ii) n ≈ 1.6 for 0 ∈ [3.8, 10] %. The biphasic 5 wt% CNC suspension did not exhibit angular frequency dependency, but two distinct regions are observable at = 0.6 rad/s, however their characteristic slopes are nearly identical. This could be interpreted as a result of approaching the transition to the LC phase. For the liquid crystalline suspensions, 6-8 wt% CNC, Fig. 6e-g in the range of strain amplitudes, where multiple scaling regions are generally identifiable, there is a clear quadratic scaling MAOS region, n ≈ 2 . Interestingly, a distinct − dependent region at the transition to LAOS could be inferred for 6-8 wt% CNC. Based on analysis above, we hypothesize that, two scaling regions (regions with different n exponents, I 3∕1 ∝ n 0 ) with a significant angular frequency dependence is related to isotropic phase (2 wt%) and originates from a weak gel-like behavior. In contrast, the biphasic 3-5 wt% CNC, exhibited lowerdependency, with two identifiable scaling regions. We note that the onset of the gel behavior does not necessarily correspond to the isotropic-biphasic transition, as the formation of tactoids can occur at lower or higher concentrations depending on the aspect ratio (Honorato-Rios et al. 2018). In addition, while the two scaling regions remained distinct for = 0.6 rad/s, their scaling was increasingly similar, signaling the transition to LC. Furthermore, a mostly quadratic MAOS region characterized the LC suspensions, 6-9 wt% CNC, followed by a a weak −dependent region at the MAOS-LAOS transition. In turn, -dependency in the MAOS region could indicate, that LC 9 wt% approaches a glassy state.
The elastic and viscous Lissajous-Bowditch (LB) diagrams for selected concentrations representative of the different CNC phases are given in Fig. 7. The LB diagrams are summarized in a Pipkin space which provides a "rheological fingerprint" that is the intra-cycle response as a function of the applied and 0 . At low strain amplitudes, the stress response is affected by instrumentation noise, which is most observable for the isotropic phase due to the lower torque values measured, Fig. 7a. The elliptic curves, representing the linear viscoelastic response of the material, are observed at small strain amplitudes, see Fig. 7a for ≤ 1 rad/s and 0 = 1 %. As the strain amplitude increases, the LB loops appear to be increasingly distorted, which indicates the occurrence of intra-cycle nonlinearities. From qualitative point of view, the isotropic LB diagrams, 7b, stand out compared to the biphasic and the liquid crystalline phase LB diagrams, 7c, that are comparatively more similar to each other. The shape of the loop also depends on the imposed , most readily observable for 2 wt% CNC, Fig. 7a, as expected. This is in agreement with the LVE dynamic moduli and I 3∕1 results, pointing towards the existence of a weakly gelled structure. We note that the viscous LB diagrams show self-intersecting loops for the biphasic and liquid crystalline phases. Self-intersecting loops have also been observed by Moud et al. (2021) for highly concentrated hydrogels with cellulose nanocrystals and was attributed to an overshoot in shear stress response. The strain-stiffening (S) and shear-thickening (T) dimensionless ratios are compiled in Fig. 8 for all investigated concentrations. As expected, at low strain amplitudes, within the linear viscoelastic region, S, T ≅ 0 . The lowest CNC concentration, 2 wt%, due to its low viscosity, was affected by a low torque, resulting in scattered datapoints at low 0 , Fig. 8a.
All suspensions exhibit a dominant intra-cycle strain-stiffening nonlinear behavior, where S > 0 , independently on applied . There is also a notable increase in strain-stiffening between the lowest and the highest CNC concentrations, i.e. compare Fig. 8a with 8 h. We note that the positive sign in S is an apparent stiffening because at large strain amplitudes both G ′ L , G ′ M soften, see Eq. 3 and Fig. S3, however the nonlinear response is dominated by G ′ M softening more quickly than G ′ L (Khandavalli and Rothstein 2015). For the isotropic 2 wt% CNC S is a function of the applied in the nonlinear region, mirroring the G ′ , G ′′ and I 3∕1 behavior. Interestingly, for the 6-9 wt% CNC suspensions there is a slight but persistent S < 0 region, see I) in Fig. 8e-h, i.e. a local intra-cycle strain softening behavior, around 0 = 2% before transitioning to S > 0 with increasing 0 .
For all investigated concentrations, a dominant nonlinear viscous intra-cycle shear-thinning behavior is observed, where T < 0 . Interestingly, the -dependence, that has been observed in G ′ , G ′′ and I 3∕1 , is mostly limited to the strain stiffening ratio S in the stress decomposition analysis. At 3 wt% CNC, we noted a slight T > 0 region, see II) in Fig. 8b, i.e. a local intra-cycle shear-thickening behavior region around 10 1 %, whereafter T drops, resulting in intracycle shear-thinning, T < 0 . By comparing the relative magnitude of ′ L and ′ M , see Fig. S4 in the Supplementary Information, it can be concluded that they both experience a local thickening before thinning at large strain amplitudes. The strain-rate thickening behavior, where T > 0 , at medium strain amplitudes was caused by stronger nonlinearity of ′ L than ′ M . In turn, a negative sign in T at increased strain amplitude was dominated by ′ M thinning more quickly than ′ L (Khandavalli and Rothstein 2015). This behavior is more readily observable for concentrations >3 wt% CNC, Fig. 8c-h, with the T > 0 region covering approximately 1 decade of strain amplitude for the highest investigated concentration, 9 % CNC, Fig. 8h. Local intra-cycle shear-thickening behavior has also been observed in other complex fluid systems (Khandavalli and Rothstein 2015; Sun et al. 2011). A dual behavior, where at lower strain amplitude intra-cycle shear-thickening region is observed, could indicate jamming of the gel-like microstructure, while intracycle shear-thinning at higher strain amplitude could be associated with the orientation or yielding of CNC particles in the flow direction (Wagner and Brady 2009;Kamibayashi et al. 2008). Since the behavior is more prominent with increasing CNC concentration from the biphasic region onward, it is apparent that the liquid crystalline domains play a major role in the microstructural jamming. The 6-9 wt% CNC suspensions exhibited an apparent dependence on the imposed angular frequency at very high strain amplitudes ( 0 > 100% ). We assume that this could be an artifact affected by the self-intersecting visible in LB diagrams at the same . Based on stress decomposition analysis, we hypothesize that, similarly to I 3∕1 , a significant angular frequency dependence is related to the isotropic phase (2 wt% CNC) and relates to the weak gel-like network present at that concentration. In turn, a dual behavior exhibited by the T parameter, where at medium strain amplitudes T > 0 and at large strain amplitude T < 0 , can indicate the presence of a biphasic phase for 3-5 wt% aqueous CNC suspensions. Furthermore, the intra-cycle strain softening behavior 10 2 10 1 1 1 0 1 10 2 10 3 9 wt% CNC Fig. 8 (continued) at medium strain amplitude, where S < 0 , effectively marks the transition from the biphasic to the liquid crystalline phase for 6-9 wt% CNC. Finally, we must also note that the present analysis cannot exclude the existence of slip artefacts (Hubbe et al. 2017). Wall slip has been evidenced in nonlinear oscillatory shear tests for crosslinked hydrogels (Kamkar et al. 2021). Interestingly, in the scientific literature concentrations higher than investigated in this study, e.g. 10 wt% CNC, were found not to exhibit wall slip (Derakhshandeh et al. 2013). Danesh et al. (2021) showed that concentrations 1-5 wt% CNC obeyed the noslip boundary condition, while the addition of NaCl caused wall slip.

Rheo-PLI in oscillatory shear
The PLI space time-diagrams corresponding to strain sweep tests at = 0.6 rad/s are shown in Fig. 9, with individual cycles from selected strain amplitudes shown in Fig. 10. The same dataset for = 4 rad/s can be found in Figs. S5 and S6 in the Supplementary information. The onset of nonlinear viscoelastic behavior is overlaid on the images by marking the strain amplitudes for which S, T are non-zero.
We briefly note that a quantitative image analysis would not be reliable due to variable illumination conditions during and between the tests and that the onset of nonlinear behavior is a strong function of the data analysis method, with the onset differing between the linear viscoelastic dynamic moduli, I 3∕1 and S, T (Meurer et al. 2021). Overall, however, linear viscoelastic behavior generally corresponds to Fig. 9 Space-time optical visualizations of the birefringence pattern development and phase transitions in dynamic strain sweep measurements at = 0.6 rad/s no observable birefringence effects, e.g. see 0 = 1% in Fig. 10, except at the higher end of the concentration range, 8-9 wt%, where there is a significantly higher density of the mesophase while approaching the glassy state. Interestingly, within the T overshoot region, see especially 0 = 50% in Fig. 10, although intra-cycle birefringence variations are observable there are no significant color distortions. This Fig. 10 Intra-cycle birefringence patterns at selected strain amplitudes for = 0.6 rad/s, corresponding to the space-time diagrams in Fig. 9. The corresponding elastic and viscous LB diagrams, see also Fig. 7, are also included supports the premise that the local intra-cycle shearthickening behavior is determined by microstructural jamming with no major structural re-arrangements. For increasing strain amplitudes and concentrations, increasingly significant intra-cycle color distortions are observed, corresponding to S > 0 and T < 0 , Fig. 10. The distortions are induced from the free surface of the flow domain (lower end of the images) and propagate inwards, as the local shear rates are proportional to the radial position in plate-plate flow (Kádár et al. 2020a). It is apparent that the color progression within a cycle could mimic the elastic LB diagrams, e.g. compare 5 wt% CNC 0 = 100% and 0 = 500% in Fig. 10: a more gradual color progression appears to corresponds to smoother more elliptic shaped elastic LB, whereas a more abrupt progression corresponds to more square-like elastic LBs. Likewise, for high strain amplitudes and CNC concentrations, intra-cycle color asymmetries are recorded that could be contributing factors to the oscillations recorded around the minimum strain in elastic LB diagrams and the self-intersecting loops in the viscous LB diagrams, see also Fig. 7. Interestingly, at the onset of the MAOS region in I 3∕1 no significant intra-cycle color variations were recorded, suggesting that the nonlinear behavior originate from both the isotropic parts of biphasic suspensions and interfaces between mesogens in the LC phase, as well as, with increasing strain amplitudes, changes in the orientation of the liquid crystalline domains (Fazilati et al. 2021).

Conclusions
The linear and nonlinear viscoelasticity of isotropic (2 wt% CNC), biphasic (3-5 wt% CNC) and liquid crystalline ( ≥ 6 wt%) phases in CNC water suspensions were investigated in this study. The phases were determined based on a combined analysis of polarized optical microscopy, linear viscoelastic oscillatory shear, steady shear viscosity functions and combined rheo-PLI experiments. The results were then compared with Fourier-transform Rheology and stress decomposition analysis to determine the influence of the particular CNC phase on the nonlinear material parameters. The relationship between the CNC phase and the rheological parameters investigated is summarized in Table S3 and Fig. S7 in the Supplementary information. The intra-cycle and inter-cycle nonlinear material parameters investigated in this study were capable of identifying all CNC phases, with changes in I 3∕1 scaling behavior with 0 distinguishing all investigated phases, the shear-thickening overshoot in T identifying biphasic suspensions, and a local slight strain-stiffening in S appeared to be associated only with the liquid crystalline phase. Therefore, based on the present results it can be concluded that nonlinear materials parameters from FT rheology and stress decomposition could be a more sensitive framework to identify CNC suspensions phases from rheological tests.