Thixotropy of cellulose nanocrystal suspensions

The thixotropy of cellulose nanocrystal (CNC) water suspensions is intrinsically dependent on the hierarchical structure of the suspension. The diverse hierarchies that comprise individual CNC nanoparticles and mesophase liquid crystalline domains, chiral nematic and nematic structures, contribute selectively to the rheological material response. Here, we combine rheology with polarized light imaging (PLI) to elucidate the thixotropic behavior of CNCs suspended in water. The simultaneous monitoring of PLI and rheological tests enables the observation of mesogens and their orientation dynamics. Creep, dynamic time sweep, ramped hysteresis loop, and thixotropic recovery tests combined with PLI aim to differentiate the contribution of the different hierarchical levels of CNC suspensions to their thixotropy. The range of concentrations investigated comprised biphasic (4 and 5 wt. %) and liquid crystalline phase suspensions (6, 7, and 8 wt. %). The CNC suspensions exhibited complex thixotropy behavior, such as viscosity bifurcations in creep tests and overshoot in ramped hysteresis loop tests. The restructuring and destructuring appeared to correspond to different levels of their hierarchical structure, depending mainly on the phase, in agreement with previous studies. Restructuring was attributed to re-organizations of an individual CNC, e.g., in the isotropic fraction of biphasic suspensions and at the mesogen interfaces in liquid crystalline phase suspensions. However, by increasing liquid crystalline fraction in the biphasic concentrations, restructuring could also involve mesogens, as indicated in the creep tests. For flow conditions above the yield stress, as evidenced by the ramped hysteresis and thixotropy recovery tests, destructuring was dominated by orientation in the flow direction, a process that is readily observable in the form of PLI “ Maltese-cross ” patterns. Finally, we show that a simple thixotropy model, while unable to capture the finer details of the suspension ’ s thixotropic behavior, could be employed to predict general features thereof. © 2021 Author


I. INTRODUCTION
Cellulose nanocrystals (CNCs) are commonly derived from cellulose-rich sources through acid hydrolysis [1]. They have shown potential, among others, reinforcements in composites [2], food applications [3], coatings, transparent films, flexible displays [4,5], substrates for printable electronics [6], and in drug delivery [7]. High aspect ratio, low density, mechanical strength, high modulus of elasticity, unique optical properties, renewable origin, and being nontoxic are properties that are often exploited in the target materials. CNCs are rodlike particles with lengths between 100-300 nm and diameters typically between 4-10 nm [8,9]. An example of atomic force microscopy (AFM) imaging of the CNC particles used in this study is presented in Fig. 1(a). The analysis of the CNC dimensions is challenged by the limited lateral resolution of AFM, the difficulty to capture the individual CNCs on the AFM substrate, and the challenge to define what is an individual particle. Recently, the CNC diameter used also in this study was determined to be 4.1 ± 1.0 nm utilizing the height profiles of AFM topography characterization [10]. CNCs self-assemble in water dispersions and exhibit chiral and nematic liquid crystal orders that are often visualized by birefringence of polarization patterns [11,12]. The optical and mechanical properties of CNCs are dependent on the microstructure and their orientation in the material of interest. CNC suspensions constitute dilute or concentrated isotropic domains, whereas above a critical concentration, CNC tactoids are formed resulting in a biphasic system. Increasing the concentration leads to the fraction of the isotropic regions in the biphasic system progressively diminishing with more CNCs extending into chiral nematic and nematic domains [13][14][15]. This ultimately leads to the formation of a liquid crystalline order in the suspensions and at higher concentrations to a gel state [16]. A schematic overview of the isotropic, biphasic, and liquid crystalline phases is presented in Fig. 1(b). Treatment of the dispersion with ultrasound temperature, altering the content of surface sulfate half ester groups, and ionic strength can change the microstructure [17][18][19]. The formation of liquid crystalline domains c) Author to whom correspondence should be addressed; electronic mail: roland.kadar@chalmers.se a) can be readily observed as birefringence patterns using polarized light optical imaging combined with rheological measurements [15]. Elucidating the rheological properties in relation to the phase behavior of CNC suspensions is critical to exploit the inherent properties of CNC in applications depending on viscosity and flow.
Thixotropic yield stress fluids are used in a multitude of industries, including foods, cosmetics, pharmaceuticals, and paints. Thixotropy refers to the time-dependency of complex fluids, which means the material properties depend not only on the applied shear rate or shear stress but also on the material deformation history [20][21][22][23][24]. Thixotropy and yield stress are important properties when a material is required to flow during dispensing but then must remain set to allow a defined pattern and shape without sagging or slumping after dispensing [20,25]. The microstructural origins of the thixotropic yield stress behavior lie in the balancing between interparticle forces, Brownian motion, and hydrodynamic interactions [26]. A number of rheological tests can be employed to evidence thixotropic behavior. For example, in simple soft materials, the shear viscosity is a function of shear rate or shear stress. Consequently, when subjected to a ramped upward shear sweep followed by a ramped downward shear sweep, the recorded viscosity values are equal. In contrast, thixotropic materials exhibit a hysteresis behavior, i.e., the downward sweep of the flow curve shows a lower viscosity than the upward partly because the microstructure is broken and does not recover within the time scale of the experiment [27][28][29][30][31]. The complex rheological behavior of thixotropic materials can be regarded as a competition between aging (restructuring) and rejuvenation (destructuring) during flow. The microstructure is destroyed at a high imposed shear stress/rate, which is observed as a viscosity decrease in time.
On the contrary, the microstructure rebuilds at zero or low imposed shear stress/rate. In creep tests, by imposing constant shear stress on the sample, aging and rejuvenation can be readily observed. For stresses smaller than the yield stress, restructuring dominates over destruction, and the viscosity increases in time. In contrast, for stresses larger than yield stress, the microstructure is broken down, and the viscosity decreases in time. This contest leads to a viscosity bifurcation around the critical stress, which is the yield stress of the fluid [28,31,32].
The interaction between flow and the CNC suspension phases results in a complex rheological behavior [14]. To further understand the underlying mechanisms in the CNC suspension flow behavior, efforts have been directed to the investigation of CNC suspensions using combined rheological methods, in order to probe time-space scales that complement bulk rheological properties [15]. In particular, due to their phase behavior and resulting optical properties, combining rheometry with polarized light imaging (PLI) is a facile method to characterize CNC aqueous suspensions [17][18][19][33][34][35][36]. Coupled rheo-optical techniques, simultaneously using rheological and optical visualization methods, can enable the connection between the microstructure, flow alignment, and the bulk rheological material response [37,38]. Any change in distance and angle between rodlike particles in the nematic phase or pitch length in the chiral nematic phase results in a variation of scattered beam wavelength according to Bragg's law [12,[39][40][41][42]. In a transparent medium with liquid crystalline ordering or during its realignment in flow, an incident light beam becomes directional and gives an indication of the changed orientation when visualized via polarization patterns [35]. Here, we must distinguish between CNC-CNC interactions and the mesoscale chiral nematic/nematic structures and their interaction in relation to the bulk rheological properties and the observation space scale of PLI. It would be reasonable to assume that bulk rheological properties can be regarded as the result of stochastic interactions between CNCs as well as their mesogens. However, the birefringence patterns at the PLI observation space scale can be regarded as the result of mesogens distributed in a polydomain structure within the measurement gap interacting with the incident polarized light. The transformation from a state consisting of multiple mesophase orientation directions to a structure where they are aligned in the same direction, unidirectional pattern, can be assigned to organization toward nematic alignment, an organization that can be induced by sufficiently high shear rates [15,36]. Shafiei-Sabet et al. evidenced two microstructural transitions with increasing CNC concentration from an isotropic to a chiral nematic liquid crystal and then to the gel state. These transitions were identified based on the viscosity functions and PLI "fingerprint texture." A three-region behavior was observed, usually attributed to CNC concentrations in the biphasic phase, however, also observable in concentrations attributed to the liquid crystalline phase in other studies [33,34]. A single shear thinning behavior was associated with concentrations in the gel phase [17][18][19]. Haywood et al. investigated the combined effects of concentration, shear, and microstructure relaxation on the alignment and optical properties of CNC suspensions and films by using a shear cell coupled with optical polarized light microscopy [33,34] and rheology coupled with small-angle neutron scattering (SANS) [34]. The study noted that CNC suspensions showed greater alignment after shear and less relaxation after the cession of shear in the liquid crystal phase, in comparison to isotropic and biphasic dispersions [33,34]. Hausmann et al. investigated the orientation dynamics, which is a function of concentration and applied shear rate, with the help of a rheo-PLI in the reflection mode, aiming to tailor orientation during 3D printing [36]. The yield stress of the suspension important for maintaining the shape of the 3D printed parts and, therefore, flowinduced stresses are required to overcome the yield stress for orientation. Consequently, the study focused on high CNC concentrations (>10 wt. %) corresponding to the liquid crystalline and glassy state [36]. In this concentration region, the yield stress was found to obey a power law dependence on the CNC concentration with an exponent of 3 [36]. Recently, Kádár et al. [35] examined the microphase transition sequences of CNC during steady shear and oscillatory shear strain sweep tests through coupled rheology and PLI. Their custom transmission mode setup allowed the visualization of vibrant birefringence patterns, and thereon several transitions based on the flow were identified with increasing shear rate. From a (quasi-) linear region, a nonlinear transition region followed by general orientation in the flow direction, as identifiable through the Maltese-cross birefringence pattern, was observed in both biphasic and liquid crystalline phases investigated. A unique radially periodic birefringence pattern corresponding to the biphasic state was observed [35]. The study evidenced the importance of the shear history through the repetition of steady shear tests with a relaxation time between tests. While the birefringence patterns recorded at the beginning of the first tests were not recoverable, this caused minor changes in the viscosity function at low shear rates in the repeated tests. However, a quantification of the CNC suspension thixotropic behavior was not further pursued. The time-dependent behavior of CNC suspensions has been reported in several studies [43][44][45][46]. Derakhshandeh et al. [46] reported the aging behavior of a 10 wt. % CNC suspension, corresponding to the liquid crystalline phase, under both oscillatory and steady shear flow by rheological measurements coupled with light scattering echo (LS-echo) and PLI. Aging was correlated with the increase in storage modulus with time in dynamic time sweep tests and in transient viscosity at constant shear rates. In both experiments, the aging behavior was related to the induced shear stresses being close to or below the yield stress [46]. The LS-echo data also confirmed the absence of slip at the flow domain boundaries [46]. Xu et al. [43] focused on the phase behavior of CNC suspensions as a function of salinity and concentration. Based on their rheological characterization, liquid, gel, and soft glass phase transitions were distinguished. The yielding behavior, as determined from steady shear stress sweeps, was used to distinguish between soft repulsive glasses (high CNC concentration and low salinity) and soft attractive glasses (high salinity). To distinguish between liquid and solid phases, the thixotropy, creep, and yield stress behavior was employed, specifically at the solidliquid lower limit. In particular, the thixotropy of the CNC suspensions was determined based on forward and backward shear stress sweeps.
We focus on the time and shear history dependency of CNC suspensions in water through coupled rheology-PLI. Previously, we had focused on the relationship between phase behavior and alignment in CNC suspensions, with evidence of thixotropic behavior in steady shear data [35]. The study complements previous work with combined rheological techniques [46] by focusing on a broader range of CNC concentrations (4-8 wt. %), comprising the biphasic and liquid crystalline phases. In addition, the work of Xu et al. [43] is complemented by using a combined rheology-in situ PLI in relation to the phase behavior. To the best of the authors' knowledge, this is the first study focusing comprehensively on the thixotropic behavior of CNC suspensions. Four tests that support thixotropy quantification are applied: hysteresis loop, creep, dynamic time sweep, and thixotropy recovery. The birefringence and the rheological measurements of CNC closely track how CNC organization changes under shear flow. In addition, by comparing the thixotropic behavior recorded from rheological measurements with simultaneous macro-scale PLI, we aim to assign the thixotropic behavior to the hierarchical structure, i.e., mesogen polydomain structure and orientation or to individual CNC rearrangements. Finally, depending on the concentration of CNC suspension and subsequent phase assemblies, we explore the variation of particle interaction and network density causing a diverse thixotropic behavior.

A. Materials
Suspensions of 4, 5, 6, 7, and 8 wt. % CNC from CelluForce (NCC ® Montreal, Canada) were dispersed in Milli-Q water (Millipore) homogenized by a bench shaker for 24 h and treated for 30 min in an ultrasound bath (VWR, Radnor, PA, USA) operating at 178 W. Since the suspension rheological properties are greatly affected by the preparation procedure, the mentioned preparation steps have been performed similarly for each sample and before each measurement. CNC topography and phase images were recorded using AFM (NTEGRA, NT-MDT) in the semicontact mode using Tap300AI-G tips. According to the analysis of Navarro et al. [10], we estimate the length of the CNCs to be in the range of 100-300 nm and we use 4 nm as the diameter to estimate their aspect ratio at 25-75, see also Fig. S1 [61] where the difficulty to capture the individual CNCs on the AFM substrate is illustrated.

B. Experimental setup
All rheological measurements were performed on an Anton Paar MCR 702 Twin Drive rheometer, using the transparent glass parallel plate geometry (diameter 43 mm) with a measurement gap of 1 mm. The tests were performed in a single motor-transducer configuration using a modified version of the rheo-optical visualization setup of Kádár et al. [35], based on the P-PTD200/GL accessory. All experiments were performed at ambient temperature, 23°C. A scheme of the experimental setup is shown in Fig. 2. Two linear polarizers were placed above the upper plate geometry and below the lower plate at a 45°relative orientation (0, 45°). HD (1920 × 1080 pixels) video recordings at an acquisition frequency of f acq ¼ 30 fps (frames-per-second) were performed 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). The observation area consisted of the entire geometry circumference, as shown in Fig. 2

C. Methods
Prior to each test, a 4-min waiting time was applied to allow for sample relaxation after setting the measurement gap. Steady shear and oscillatory shear strain sweep tests were performed to help assess the phase behavior of the suspensions using the same procedures described by Kádár et al. [35]. Four different experimental methods were employed to quantify the thixotropic behavior of the suspensions: creep, dynamic time sweep, hysteresis loop, and thixotropy recovery tests.
In creep tests, different constant shear stresses, smaller and larger than the yield stress, are applied and the transient viscosity is measured as a function of time. In this work, the applied shear stresses were in the range of 5-60 Pa.
In dynamic time sweep tests, the complex viscosity (or storage modulus) is monitored in time at constant shear strain amplitude and angular frequency within the linear viscoelastic range (small amplitude oscillatory shear). In our experiments, a strain amplitude of 2% and frequency of 1 Hz were imposed on the sample and the complex viscosity was monitored over time. The imposed shear strain and frequency are constant and sufficiently small to maintain the microstructure, and hence aging is measurable. Thixotropy can be identified by an increase in the complex viscosity in time [47].
Hysteresis loop tests are a typical method to assess the thixotropy of soft materials, where an imposed shear rate is ramped upward followed by a downward ramp. In this work, the imposed shear rate was increased from 10 −3 to 10 2 s −1 and then decreased over the same range with a waiting time of 10 s per data point. Corresponding to each value of the imposed shear rate, the shear stress (or transient viscosity) was monitored.
Thixotropy recovery tests typically comprise three-step shear intervals with the goal of determining the capability of the material to rebuild its structure after destroying it. In the first step, a low shear stress below the yield stress is imposed on the sample. In the second step, a shear stress above the yield stress is employed to break down the material's structure. In the last step, the shear stress is again dropped to the low value of the first step and the rebuilding of the structure is followed for a specified time. In the present experiments, a pre-shear of _ γ ¼ 100 s À1 for 1 min was applied to erase the samples' history [27]. The low and high shear stress intervals for each concentration were selected based on the creep tests. increase the viscosity by restructuring and recover the initial viscosity can be used to quantify the thixotropic behavior.
In all of the rheological experiments described above, the PLI patterns were recorded using the setup described. From the video recordings, space-time diagrams were constructed by extracting one line of pixels at a fixed position out of each still frame and appended to a new image having the x-dimension corresponding to the experimental time [the time between two pixels in the x-direction is Δt ¼ 1/f acq ¼ 0:03(3)s] and the y-dimension corresponding to the length L in Fig. 2, see also [35]. For increased clarity, the height of the space-time diagrams has been extended. The space-time diagrams thus constructed reveal how the birefringence patterns at a specific location vary in time during the various imposed shear flow conditions.
We note that while we cannot exclude the possible occurrence of apparent slip in our experiments, interestingly, CNC suspensions of 10 wt. % were shown by Derakhshandeh et al. [46] to obey the no-slip boundary condition in similar experimental conditions.

D. Numerical modeling
The experimental results were compared with a structural kinetics thixotropy model. The goal of the modeling was to provide a first estimation of the reach of typical thixotropic theoretical descriptions. Comprehensive overviews of the different models available can be found, for example, in the article by Mewis and Wagner [23]. An exhaustive evaluation of feasible models is outside the scope of this work. The modeling of the thixotropic characteristics of the CNC suspensions is, therefore, approached by choosing a relatively simple model, which quantitatively describes the phenomena observed in the experiments. Thus, a model was chosen that describes an instantaneous shear rate dependence of the apparent viscosity through a power law component as well as long-time thixotropic changes through a structure parameter. The thixotropic response of the CNC suspensions while undergoing shear at different concentrations was modeled through the equation where η is the shear viscosity, _ γ is the shear rate, K 1 and K 2 are fitting constants, and m is the power law index. The parameter λ [ [0, 1] represents the thixotropic state of the material and is described by the rate equation where k 1 and k 2 are the model parameters related to destructuring and restructuring of the material, respectively. Equation (2) and similar combinations of terms to model destructuring and restructuring have been used by several authors to model the thixotropic state of various materials [23]. Equation (1) allows for the combination of a thixotropic component to the viscosity through the constant K 2 as well as an instantaneous dependence on the shear rate via K 1 .
Furthermore, the form of Eq. (2) ensures the boundedness of λ, which simplifies the parameter estimation procedure. The rate Eq. (2) can also be expressed in terms of shear stress as where _ γ ¼ (σ/K 1 þ K 2 λ) 1/m , and can thus be used when shear stress is the input variable rather than the shear rate.
To summarize, the following steps were performed. (i) The model parameters (k 1 , k 2 , K 1 , K 2 , m) were determined for the ramp hysteresis loop tests by minimizing the integrated squared error between the experimental and numerical stress curves in log-space. This particular objective function for the minimization was chosen as it accounts for differences in stress magnitude as well as for different amounts of points in different flow regimes. (ii) Next, the model with the previously fitted parameters was used to predict the thixotropy recovery experiments. (iii) Finally, the creep tests were simulated using the initial parameters, after which an optimization of the results was performed to correct for the factor k 2 based on the observed experimental yield stress range.
We note that the model is a generic continuum description of a generalized thixotropic behavior and that other methods that can capture the self-organization of the CNC into chiral nematic and nematic phases would be more appropriate. However, we note that the theoretical description of liquid crystalline systems remains a challenge [48][49][50][51][52].

III. RESULTS AND DISCUSSION
We begin by discussing the phase behavior of the investigated CNC suspensions followed by a detailed description of the individual thixotropy tests. We focus first on the slow restructuring-destructuring dynamics in relation to the yield stress. Thereafter, the thixotropic behavior is discussed in the context of flow conditions well beyond the yielding of the samples. The relevance of the numerical modeling is discussed at the end of the section. To interpret the results, it is necessary to refer to the hierarchical nature of the CNC suspensions. In general, we refer to CNCs as being the elementary 1D prolate particles while we refer to the self-assembled tactoids/nematic/chiral nematic CNCs as mesogens.

A. CNC phase behavior
All concentrations showed a three-region steady shear viscosity behavior [53,54], which is presented in Fig. 3(a), albeit weakly, especially for the higher concentrations. While assigning the three-region behavior exclusively to a particular phase behavior has yet to be elucidated, it has previously been reported for both biphasic and liquid crystalline phases and confirmed through rheo-SANS experiments [34]. The strain sweep tests [ Fig. 3(b)] show that all concentrations exhibit a rheological gel behavior, i.e., G 0 . G 00 . However, concentrations above 5 wt. % CNC exhibit a significant weak strain overshoot behavior (local increase in G 00 at the transition to nonlinear behavior). Images (still frames from video recordings) of the circumferential birefringence patterns show that while all suspensions display coloring during tests, concentrations 6-8 wt. % displayed large observable isochromatic areas before the beginning of rheological testing ( _ γ ¼ 0) (Fig. 4).
Based on the considerations above, we hypothesize that for the given preparation procedures concentrations, 4-5 wt. % are biphasic CNC suspensions, comprising an isotropic and a liquid crystalline phase, and 6-8 wt. % are liquid crystalline CNC concentrations, where the entire domain is comprised of mesogens. Figure 5 shows the temporal evolution of the transient shear viscosity at different applied shear stresses for all CNC suspension concentrations. Overall, viscosity bifurcations were recorded for all concentrations, a behavior reported for a broad variety of systems [27,32]. However, the critical stress range differed depending on the presumed isotropic (≤5 wt. %) or liquid crystalline (≥6 wt. %) CNC phase, see Table I and also Fig. S2 [61]. For liquid crystalline CNC suspensions, the yield stress range is higher than for biphasic suspensions. The yield stress range determined for 7-8 wt. % CNCs is similar to the yield stress determined by Hausmann et al. [36] for 10 wt. % CNCs and is slightly lower than the yield stress determined by Xu et al. [43] for 11 wt. % CNCs. We also note that extrapolations of the power law determined for yield stress vs CNC concentration [36] would underestimate the yield stress ranges determined for the CNC concentrations reported in our study. The differences are likely due to varying sulfate half ester functional group content [55], aspect ratios, and differences in the sample preparation.

B. Creep tests
In Figs. 5(c) and 5(d), i.e., 6 and 5 wt. % CNC, oscillations can be observed for shear stresses in the intermediate regime between destructuring vs restructuring dominated states. This unexpected oscillation has also been observed in other soft glassy materials, and it refers to instability and coupling between elastic and inertia effects at the solid-liquid transition, which possibly comes along with local structural rearrangement [56][57][58]. Figures S3-S7 in the supplementary material [61] contain representative space-time images corresponding to the creep tests in Fig. 5. Below the yield stress, we can distinguish between the liquid crystalline and biphasic suspensions. Liquid crystalline concentrations did not exhibit any observable birefringence pattern dynamics within the observation time, e.g., see σ , 30 Pa in Figs. S3-S5 [61]. The biphasic 5 wt. % CNC suspension exhibited a gradual birefringence pattern shift (Fig. S6 [61]). However, the biphasic 4 wt. % CNC suspension (Fig. S7 [61]) exhibits a similar lightbrown, blue coloration for σ 20 Pa, with no significant , could suggest that structural changes associated with the thixotropic behavior are the result of restructuring mainly involving the primary CNC particles, rather than their mesogens, that can be captured by the bulk rheological properties but not by the birefringence patterns at the observation scale. The increase in CNC concentration is known to increase the liquid crystalline phase [59], and thus it could be inferred that for 5 wt. % the thixotropic behavior is due to significant structural changes in the mesophase. This could be attributed to either a change in orientation of the chiral nematic/nematic units to structural adjustments of CNCs within chiral nematic/nematic units or to a combination of both. Interestingly, for 5 wt. % CNC and intermediate shear stresses, σ ¼ 20 Pa, the color patterns were unchanged, which indeed corresponds to almost constant nominal transient viscosities in Fig. 5(d). For shear stresses larger than the yield stress, the birefringence patterns change in time, as the conditions to significantly alter the sample structure are met. The CNCs are negatively charged in conditions used here and the CNC interaction is repulsive. Increase in the concentration of the likely charged CNCs has been reported to lead to a transformation from a viscous liquid to a (repulsive) glass phase [43]. For samples exhibiting yield stress and G 0 . G 00 in the high concentration range, particle dynamics are considered arrested due to particle friction, jamming, and squeezing. The probability of van der Waals interactions is low due to the dominant repulsion between the particles. The structure at rest is considered to be determined by the gellike network and its density. The functional groups that reside on the CNC surfaces contribute to the network structure via hydrogen bonding. Shear can lead to progressive decrease in hydrogen bonding density that cannot be recovered after the flow due to the repulsion dominating the restructuring [45]. Hence, when shear is applied, destructuring occurs. Transition from low to high color wavelengths in the spacetime diagrams could, therefore, be regarded as indicative of changes in the characteristic distance between CNC within the mesophase, especially at the interface between mesogens [43], as well as the orientation in the flow direction of the mesophase chiral nematic/nematic units [35]. Above the yield stress region, space-time diagrams can also appear in the form of   Fig. S8 [61]), and nonuniform patterns that are circumferentially very similar (see 6 and 7 wt. % in Fig. S8 [61]). In Fig. 6, the space-time diagrams of the largest imposed shear stress are compared. The birefringence color progression with time is overall relatively similar between the tests and suggests similar flow-induced conditions that favor orientation in the flow direction. Overall, the bifurcation behavior is in broad agreement with similar systems dominated by mostly repulsive interactions between the primary elements/particles, such as Carbopol gels, repulsive latex suspensions, and nonadhesive emulsions [60].

C. Dynamic time sweep tests
The evolution of complex viscosity by applying a constant oscillatory shear strain is presented in Fig. 7. At high CNC concentrations, a significant decrease can be observed for t , 1000 s. Generally, at high CNC concentrations, i.e., liquid crystalline phase concentrations, particles experience substantial interparticle friction, thereby being blocked within their nearest neighboring cages [14,43]. Since all concentrations exhibit some degree of rheological gellike behavior (Fig. 3) and considering that the stresses induced in the dynamic time sweep tests are lower than the lowest stresses imposed in the bifurcation tests, we can assume a similar "locked-in" behavior even in the isotropic fraction of the biphasic concentrations, albeit weak. Thus, the applied low shear strain amplitude and frequency are not sufficient to rearrange the large-scale mesoscale chiral nematic/nematic structures, but it is sufficient for short-time diffusion and local particle rearrangement. Consequently, the decrease at the beginning of shear can be due to local rearrangement in the jamming zones and is consequently less significant for the lower CNC concentrations that have a large isotropic phase [46]. The subsequent increase in complex viscosity with time is due to restructuring, similar to the described creep tests performed below the yield stress. This is also confirmed by the space-time diagrams as they did not indicate microstructural changes over the test period, with the induced stresses not sufficiently high to induce birefringence changes even in the 5 wt. % CNC biphasic suspension. Overall, a qualitatively similar complex viscosity time dependence was recorded for all CNC concentrations. Figure 7(a) shows that the normalized complex viscosities, i.e., (jη * j À jη * j)/σ jη * j , where jη * j is the average complex viscosity and σ jη * j is the standard deviation, for different concentrations collapse on the same curve. Therefore, it could be inferred that a similar aging mechanism takes place, likely due to structural arrangements at the mesogen interfaces that are not altering the mesogen structures sufficiently to be observed through the birefringence observations.

D. Hysteresis loop tests
The hysteresis is obvious in the flow curves (Fig. 8), where the shear stress, σ(t), is recorded with respect to the shear rate, which is first increased (filled symbols) and then decreased back to the starting value (hollow symbols). During the sweep up, the microstructure is disrupted at high shear rates; therefore, structural interactions are altered in the down sweep and hysteresis is induced in the flow curves. Overall, the shapes of the hysteresis loops suggest a process that is initially dominated by the breakdown of the initial structure resulting in stress overshoot; a relative measure of thixotropy in the hysteresis loop tests is the area within the loop [23]. This can be computed by integrating the stress with respect to the shear rate on the ramp-up and subtracting the integrated stress with respect to shear rate on the rampdown, i.e., The computed hysteresis loop areas are proportional to the CNC concentration, see Table S1 [61] and Fig. 13. For CNC concentrations in the liquid crystalline phase, ≥6 wt. %, the overshoot is more complex, and for the two highest concentrations, a dual overshoot could be inferred. Overall, at higher concentrations, the gel strength and thereby gellike network density are greater and because the hysteresis is more obvious as the more frequent connections in the network are disrupted during the sweep up. Interestingly, compared to the steadystate viscosity curves, there is no clear three-region behavior identifiable for the biphasic concentrations but could be related to the upward ramp dynamics of the liquid crystalline concentrations. This can mean that the time scale associated with the development of steady-state equilibrium microstructure, as associated with the three-region behavior, is significantly longer for the biphasic concentrations than the arbitrary scale imposed by the 10 s time interval per fixed shear rate. We note also that increasing the shear rate leads to that the time required for steady-state (hti) approaches the imposed experimental time scale as ht i / 1/ _ γ. However, a three-region behavior could be inferred on the downward ramp, likely because the flow conditions resemble comparatively more the steady-state conditions in Fig. 3(a). Figure 9 illustrates the space-time diagrams corresponding to the flow curves in Fig. 8, while 3D surface plots thereof can be found in Fig. S9 [61]. Additionally, in Fig. 10, still frames extracted from the video recordings at selected shear rates for the tests in Fig. 8 are shown. For shear rates below 10 −2 s −1 , no significant pattern distortions were observed. However, significant radial intensity distortions are present for CNC concentrations >4 wt. %, most readily observable in the color intensity plots in Fig. S9 [61]. This could be attributed to structural rearrangements coupled to the shear rate gradient in the radial direction. In addition, this could be a factor contributing to the overshoot behavior in the flow curves. Increasing the shear rate, and consequently increasing the imposed shear stress on the sample, a transition from random orientation toward flow-induced alignment of the domains was recorded. By domains, here, we refer to CNC particles (isotropic phase) and tactoids for the biphasic concentrations and mesogens in the liquid crystalline phase. Domains predominantly oriented in the flow direction are signaled by the onset of the Maltese-cross-like birefringence pattern, see Fig. 10. The onset of the Maltese-cross pattern as a function of CNC concentration is listed in Table II, showing a monotonic increase of the critical shear rate with increasing concentration. As all suspensions investigated showed coloring in the Maltese-cross patterns, the assembled structures were likely composed of units that were predominantly oriented in the flow direction [15,35]. In the space-time diagrams, the Maltese-cross-like birefringence can be associated with loci where the pattern color is nearly uniform. Interestingly, further increasing the shear rate results in a change in color of the Maltese-cross pattern, which is a behavior that we had not observed previously [35] but is similar to the observations of Hausmann et al. [36].
The characteristic colors of the Maltese-cross patterns with increasing shear rate on the upward ramp (increasing shear rate) are summarized in Table II. We note that the color transitions for the biphasic concentrations appear different than the ones reported in the creep tests summarized in Fig. 6, while for the liquid crystalline concentrations, they appear more similar. This is because the maximum induced shear stresses in the ramp tests of biphasic samples were below those applied in the creep tests (ca. 30 Pa for 4 wt. % and 40 Pa for 5 wt. %), whereas the stresses induced in the liquid crystalline concentrations are comparable to those imposed in the creep tests. While the concentrations were classified as biphasic (4 and 5 wt. %) and liquid crystalline (6-8 wt. %), the flow-induced birefringence patterns appear unique in coloring succession for each concentration. Overall, increasing the shear rate shifted the characteristic color wavelength toward higher wavelengths. We briefly note that from standard RGB images, wavelength estimations are not reliable and, therefore, a direct quantification thereof from the video recordings is not possible. In the case of biphasic 4 and 5 wt. %, the increasing shear rate first shifts to wavelengths corresponding to a red color but the final state is returned to a lower wavelength (violet). The increase in wavelength could be an indication that the breakdown of network structures can lead to an increase in the spacing between the CNC within a nematic structure or, due to the large gap encompassing a polydomain structure, to the spacing between the nematic units. We also note that if at rest the mesoscale structures are dominated by chiral nematic units, some of the wavelength changes during orientation could also be related, as conjecture, to chiral nematic-nematic transitions. By decreasing the shear rate (downward ramp), the Maltese-cross was maintained, and the order of the color transformation was a reverse of the upward ramp, with the exception that the birefringence patterns do not return to the initial state. Radial heterogeneities are also observed during ramp-down: see birefringence patterns between 5.2 and 0.3 s −1 for 8 wt. %, 2.5 and 1.2 s −1 for 7 wt. %, and 5.2 and 0.3 s −1 for 6 wt. %. We note that, especially for the biphasic concentrations, as the gellike network is broken due to the orientation, there would be significant possibilities for individual CNCs to break off from the network, thus decreasing the level of interaction that the mesophase domains have with their surroundings as they orient on the upward curve or partially deorient on the downward curve.

E. Thixotropy recovery tests
Thixotropy recovery tests were performed to determine the restructuring ability after microstructure destruction and the main results are shown in Fig. 11. The shear stresses for intervals I-III were chosen based on the creep tests as follows (σ I , σ II , σ III ): (20, Fig. 3(a).
Based on the phase assignment of the suspensions, liquid crystalline concentrations, ≥6 wt. % CNC, exhibit a more pronounced transient behavior. To compare the thixotropic behavior between the CNC concentrations, the rate of viscosity recovery, which is defined as the percentage of the ratio of recovered viscosity, η 3 À η 2 , to the destructured viscosity, η 1 À η 2 , see also the graphical illustration in Fig. 11(b), was calculated for each CNC concentration and is listed in Table III.
Increasing the CNC concentration leads to a progressive dynamic evolution in qualitative terms. For the lowest concentration, 4 wt. %, the transient viscosity settles quickly for a stable rate of change in all the intervals, at the expense of a relatively low rate of viscosity recovery within the experimental time. This would be expected for biphasic samples (low concentration) as the presence of shear would be expected to distort the existing particle/mesogen network thus creating space facilitating orientation in interval II, after which during recovery in interval III the isotropic phase particularly would need extensive time to recover. In contrast, for the liquid crystalline phase concentrations that reach an aging equilibrium in interval I, the structural adjustment in interval II is impeded by the strong mesogen interactions. Consequently, the recovery after interval II is more substantial. The presence of surface charges in CNC suspension causes repulsive interaction, and after flow cessation, this interaction facilitates the recovery state. At high concentrations, the viscosity recovery ratio is more significant because of particle friction, jamming, and squeezing in the mesophase.
The space-time diagrams (Fig. 12) corresponding to thixotropy recovery tests in Fig. 11 indicate how the birefringence changes in the three-interval recovery tests of CNC suspensions. The space-time visualizations predominantly show color changes in interval II, σ . σ y , for the biphasic concentrations, while in intervals I and III for all concentrations and interval II in the liquid crystalline phase CNC concentrations (>5 wt. %), no significant changes in the color patterns were recorded. The findings are in agreement with the creep tests. For applied stresses that are below the yield stress, i.e., intervals I and III, the thixotropic behavior is detected by the transient viscosity data but is not observable at the PLI macroscopic observation scale. This confirms that restructuring is likely dominated by structural rearrangements unobservable at the observation length scale. For the biphasic CNC suspensions,

F. Considerations to the numerical modeling
The set of shear ramp experiments (hysteresis loop tests) in Fig. 8 were used to fit the parameter set (k 1 , k 2 , K 1 , K 2 , m) where λ i is the calculated value in step i, resulting from applying the constant shear rate _ γ i for a time Δt i , starting from the previous value λ iÀ1 . Using this equation, the value of λ, and thus η and σ, were calculated explicitly for the experiment for a given parameter set. The fitting was performed by generating a set of starting guesses for the parameters using Latin hypercube sampling. The viscosity and stress were then calculated for the given stepwise constant shear rate using Eq. (5), and the parameters were optimized for each starting guess by comparing the squared difference in log-space between the experimental and modeled stresses. The parameter set resulting in the smallest error was then considered optimal. For each simulation, λ 0 ¼ 0:5 was assumed as the initial value. This value was naïvely chosen, assuming that the thixotropic structure was neither fully built up nor broken down.
The resulting simulated and experimental stresses are shown in Fig. S10 [61] as a function of time and in Fig. 8 as a function of the shear rate. Comparing the modeling data to the experimental results of the hysteresis loop tests (Fig. 8), the model appears to capture the stress overshoots on the upward ramp for the biphasic concentrations in Figs. 8(e) and 8(f). However, for the liquid crystalline CNC suspensions, where the overshoot exhibits more complex dynamics with possible dual overshoots on the upward ramp, the model is unable to follow the experimental data. The same applies to the nonmonotonic stress dynamics on the downward curve. However, the simulated hysteresis loop area matches closely the experimental data, see Fig. 13 and Table S1 [61].
The power law thixotropy model parameters were then directly used to predict the thixotropic recovery under the imposed three shear stress intervals in Fig. 11. The rate equation for λ expressed in terms of the shear stress [Eq. (5)] was therefore solved. Given the preshearing procedure, the initial value λ 0 is taken as the steady-state value for constant shear rate, reading The set of previously estimated model parameters were also applied to the creep test in Fig. 12. The initial results are shown in Fig. S11 [61] and the model parameters are listed in Table S3 [61]. For this simulation, the rate equation expressed in terms of the applied constant shear stress [Eq. (3)] was solved for the structure parameter λ, and the resulting viscosities and shear rates were then calculated. A numerical simulation was performed due to the lack of an analytical solution for the rate equation expressed in terms of shear stress. A viscosity bifurcation was apparent in the simulations, however, not at the same stress level as in the experiments. We, therefore, investigated if the parameter set (k 1 , k 2 , K 1 , K 2 , m) could be corrected using the bifurcation observed in the experiments at the yield stress range. At stress levels below the yield stress, or bifurcation stress, the viscosity is increasing, i.e., λ(t) ! λ 0 , dλ/dt ! 0. Similarly, for stresses above the bifurcation stress level, the following condition should hold: λ(t) λ 0 , dλ/dt 0. At the bifurcation stress level, i.e., the limiting case in between the two cases above, the following should hold λ(t) ¼ λ 0 , dλ/dt ¼ 0, which using Eq. (3) can be expressed as the condition where σ y is the yield stress. The condition was thus used to correct the value of the parameter k 2 , with the data presented in Fig. 5. It was assumed in the calculations that the initial value was 0 , λ 0 , 1, such that it can increase or decrease, depending on the stress level. Therefore, λ 0 ¼ 0:5 was chosen. With the adjustment, a reasonable qualitative and quantitative prediction was achieved. A summary of model parameter variation with CNC concentration can be found in Fig. S12 [61].
While the simple thixotropy model employed is not capable of capturing the finest details of CNC phase behavior and interactions, the results show that it can be a reasonable approximation in most cases. Also considering the breadth of the modeled thixotropy tests, and that the model parameters were fitted only for the ramp tests and slightly adjusted for the creep tests, the results show that such modeling tools are useful for more practical applications, e.g., CNC as rheology modifiers.

IV. SUMMARY AND CONCLUSIONS
The thixotropy characterization of CNC suspensions using simultaneous rheology and PLI suggests that restructuring and destructuring occur at different length scales depending on the starting phase. In creep tests, viscosity bifurcations were recorded for concentrations between 4 and 8 wt. %. Below the yield stress liquid crystalline phase concentrations and at the lower biphasic phase concentration that were investigated (4 wt. %), restructuring was observed in the increase of the transient shear viscosity with time but not in the PLI space-time diagrams. This suggests that in the respective cases, restructuring is associated with adjustments of individual CNC rather than of mesogens. In liquid crystalline suspensions, restructuring likely occurs at the interface between the mesogens. For the biphasic concentrations, restructuring could be broadly associated with the isotropic phase; however, since the biphasic 5 wt. % CNC suspension does exhibit birefringence changes with time, it suggests that a more significant fraction of the liquid crystalline domains combined with sufficient neighboring "space" that is due to the isotropic domains being less dense can lead to birefringence dynamics during restructuring. Above the yield stress, destructuring at the PLI observation length scale exhibited complex color dynamics, while at the highest shear stresses recorded, the color progression was relatively similar for all concentrations, suggesting a similar orientation behavior for sufficiently high shear stresses. For all concentrations, dynamic time sweep tests showed similar normalized transient complex viscosity behavior and no significant color changes recorded in the PLI. This is likely to be due to that the induced stresses are below the lowest stresses considered in the bifurcation tests. Ramped hysteresis loop tests showed an overshoot behavior with CNC suspensions in the liquid crystalline phase that was significantly more complex than in the biphasic suspensions. Interestingly, the PLI pattern dynamics associated with increasing and then decreasing shear rates seem to be mainly CNC concentration-dependent with a less clear distinction between the phases. This could be attributed to the different stress ranges achieved during the tests. Liquid crystalline phase concentrations showed evidence of the possible influence of the shear rate gradient on the FIG. 13. Comparison between the hysteresis loop areas determined based on Eq. (4) for the experimental data in Fig. 7 and numerical modeling of the thixotropic behavior according to Eqs. (1) and (2) (data in Table S1 [61]). downward ramp with the occurrence of radial PLI color pattern distortions. This was attributed to mesophase jamming as it deorients. Thixotropic recovery tests confirmed largely the previous tests. In interval II, for stresses above the yield stress, orientation in the flow direction dominated the microstructural evolution with significant PLI color progression in the biphasic CNC suspensions and less so with increasing concentration in the liquid crystalline phase suspensions. Similarly, in interval III, for stresses below the yield stress, the thixotropic recovery is likely dominated by rearrangements at the interface between mesogens in the liquid crystalline phase concentrations as no recordable color changes are recorded at the PLI observation scale. The assignment of restructuring and destructuring to the hierarchical levels of the CNC suspensions follows broadly previous descriptions thereof [14]. This leads to a recommendation to use combined rheo-PLI setups as a simple but versatile method to investigate the rheology of CNC suspensions in relation to their thixotropic behavior. This may be significant, especially for the design of chemically modified CNC suspensions, where their design in relation to their thixotropy and phase behavior is significant for surface coatings or 3D printed CNC-based materials. Finally, we showed that a simple thixotropy model, while unable to capture the finest details of the CNC behavior, qualifies as a predictive tool for some applications of CNC suspensions.