Dynamic characterization of cellulose nanofibrils in sheared and extended semi-dilute dispersions

New materials made through controlled assembly of dispersed cellulose nanofibrils (CNF) has the potential to develop into biobased competitors to some of the highest performing materials today. The performance of these new cellulose materials depends on how easily CNF alignment can be controlled with hydrodynamic forces, which are always in competition with a different process driving the system towards isotropy, called rotary diffusion. In this work, we present a flow-stop experiment using polarized optical microscopy (POM) to study the rotary diffusion of CNF dispersions in process relevant flows and concentrations. This is combined with small angle X-ray scattering (SAXS) experiments to analyze the true orientation distribution function (ODF) of the flowing fibrils. It is found that the rotary diffusion process of CNF occurs at multiple time scales, where the fastest scale seems to be dependent on the deformation history of the dispersion before the stop. At the same time, the hypothesis that rotary diffusion is dependent on the initial ODF does not hold as the same distribution can result in different diffusion time scales. The rotary diffusion is found to be faster in flows dominated by shear compared to pure extensional flows. Furthermore, the experimental setup can be used to quickly characterize the dynamic properties of flowing CNF and thus aid in determining the quality of the dispersion and its usability in material processes.


Introduction
Nature has a remarkable way of structuring materials in complex hierarchies resulting in properties that are hard to achieve with synthetic materials 1 . Finding a process to produce exceptional biobased materials on an industrial scale can lead to new sustainable products with potential to replace existing products relying on extensive use of fossil raw material and fresh water. Using cellulose nanofibrils (CNF) as the building block, several examples of nanostructured high performance materials including filaments, 2 films 3 and composites 4,5 have been presented in previous publications. The general strategy to improve the material properties is to tailor the structure on a nanometer-scale through manipulation of the processing conditions.
A major obstacle of achieving full control over the material properties in the different processes, is the fact that the CNFs can vary significantly in size, shape, polydispersity and/or mechanical properties. [6][7][8][9][10][11][12][13][14] This will in turn depend on both the raw material source and method for producing CNF as well as biological activity and thermal degradation. Apart from individual fibrils being affected, it is also likely that their dynamic interactions will be influenced, which in turn affect the way fibrils rotate and align in process related flows.
Therefore, it is important to perform dynamic characterization of the actual dispersions, especially with the same concentration, in the same types of flows as commonly used in the assembly process. The main material process of interest and the starting point for the present study is the preparation of continuous cellulose filaments through hydrodynamic alignment and assembly of CNF demonstrated by Håkansson et al. 2 .
The dynamics of dispersed non-spherical nanoparticles in flows typically can be divided into two parts: advection (the externally forced motion of the particles due to the deformation of the surrounding medium) and diffusion (the internally forced motion of a collection of particles towards thermodynamic equilibrium). Furthermore, the motion of the particle can be divided into translational and rotational motion.
The translational advection of small particles such as CNF is trivial as they will follow the same motion as the fluid. Assuming an ellipsoidal shape of the non-spherical fibrils, the rotary advection due to velocity gradients of the surrounding flow is described by Jeffery 15 .
The diffusion of the particles is typically a result of Brownian (thermal) motion of surrounding molecules, but also other interactions between particles including direct collisions, electrostatics or various complex molecular interactions. The translational diffusive motion, typically leads to a uniform spatial distribution of particles. Translational diffusion is of less interest in CNF material processes, since we assume the spatial distribution in the dispersion to always be uniform and will not be affected by the flow. The important dynamics of CNF instead is described by their rotational motion. Here, the rotary advection due to velocity gradients in the flow (such as shear and extension) will result in a non-isotropic orientational distribution. This process is always in competition with the rotary diffusion process towards an isotropic state. Since rotary diffusion of CNF is a non-trivial process and depends on many factors, it is thus important to characterize this process for different CNF samples in order to understand how the fibrils will behave in process relevant flows.
Dynamic characterization of nanoparticle dispersions to target rotary diffusion is typically done with dynamic light scattering (DLS) using coherent light, commonly achieved with a laser source. By studying the autocorrelation of the scattered light, the translational and rotary diffusion coefficients can be obtained. [16][17][18][19][20] However, due to problems of multiple scattering in the sample, the dispersions are typically required to be dilute (<0.1 vol.%).
Furthermore, the technique can only be used to study the dynamics in isotropic and thus nonflowing systems. Thus, the DLS technique is not suitable to study dynamic CNF interactions in material processes where the concentration typically is around 0.3 vol.%. Rotary diffusion measurements also typically requires hours of measurements to produce a correlation function with acceptable statistics.
More suitable techniques for characterizing CNF at process relevant conditions are small and wide angle X-ray scattering (SAXS/WAXS). At isotropic (non-flowing) conditions, these techniques have previously been used to obtain relevant statistics of the individual particle shapes as well as the structure of fibril networks. [20][21][22][23][24] SAXS has also been widely used to study the orientation of elongated nanoparticles, including CNF, in flows. [25][26][27][28][29] If X-rays would be used to study the dynamics and interactions of fibrils during de-alignment of an initially anisotropic system, long irradiation times are necessary. The damage of the radiation on the sample can thus not be neglected in such a study.
It is possible to estimate the rotary diffusion coefficient in SAXS experiments from measurements at different positions in the flow. 30 This requires an assumption of the orientation distribution to be purely advected downstream such that different positions correspond to different instances of time during the diffusion process. This assumption might not always be applicable due to velocity gradients in the flow.
Recently, there has been advances also in X-ray photon correlation spectroscopy (XPCS), which is equivalent to DLS but utilizing X-rays instead of visible light and thus reducing the problem of multiple scattering. [31][32][33][34][35] This technique has been used to study the advection and diffusion of dispersed nanoparticles in flowing samples. Also XPCS suffers from the downside of radiation damage on biobased samples in non-flowing systems. To avoid radiation damage and to possibly increase the scattering contrast, neutron scattering techniques have also been used on CNF systems in a similar manner as the X-ray techniques. 20,21 However, given the low flux of present spallation sources, such an experiment becomes very time consuming.
Another way of studying the orientation of dispersed nanofibrils is to measure the birefringence of the dispersion using polarized optical microscopy (POM). [36][37][38][39][40][41][42][43] An aligned system of CNF has a higher refractive index for light polarized perpendicular to the fibrils compared with light polarized in the parallel direction. A measurement of how the birefringence of an initially aligned system decays with time is thus a measurement of how the fibrils approach an isotropic orientation distribution, i.e. a measurement of rotary diffusion. The advantages of this method is that a laser with visible light can be used as a light source with negligible damage to the sample. Additionally, regular high-speed cameras can be used to detect the transmitted light. The disadvantage of using POM is that it only provides a relative measurement of the average alignment in the sample. Consequently, it can not provide more detailed information about the orientation distribution of fibrils.
Using static measurements with SAXS and POM, Håkansson et al. 2 concluded that understanding the rotational dynamics of CNF is crucial for controlling the final material performance. In a later work 30 they further indicated that the rotary diffusion process was highly dependent on the instantaneous alignment of the dispersion. In this work, we demonstrate a flow-stop experiment using POM that can be easily used and set up for quick dynamic characterization of birefringent CNF dispersions at process relevant (actually identical) concentrations and flow conditions. By combining the dynamic characterization with SAXS measurements in the stationary flow (before stop), we can also reveal details about the initial orientation distributions giving rise to the dynamical behavior.
The POM flow-stop experiment can be used as a standard characterization tool for determining the quality of CNF dispersions. The results from this characterization method can furthermore be used to properly model the dynamics of dispersed CNF in material processes, which in turn can lead to reasonable optimization strategies to improve the performance of the final material.

Experimental section
Sample Cellulose nanofibrils (CNF) were prepared similarly as in the work by Håkansson et al. 2 from chemically bleached wood fibers (a mixture of 60% Norwegian spruce and 40% Scots pine, supplied by Domsjo AB, Sweden). The fibrils were prepared by a carboxy methylatation followed by mechanical disintegration, according to a previously reported method 44 . Briefly, carboxy methyl groups were introduced on the surface of the fibrils (degree of substitution of 0.1) to facilitate disintegration of the fiber wall during mechanical treatment. After the chemical pre-treatment, an aqueous dispersion of the fibers was passed through a highpressure homogenizer. The result is a CNF dispersion with a concentration of 6 g/L.
In an additional step, the un-fibrillated and agglomerated fiber bundles were removed from the CNF dispersion as follows. The gel like dispersion was diluted to 3.3 g/L by adding deionized water, mixed thoroughly using a mechanical mixer (12000 rpm for 10 min, Ultra Turrax, IKA, Germany) and sonication (10 min, Sonics Vibracell, USA). This diluted dispersion was then centrifuged at 5000 rpm for 60 minutes, precipitate removed and the Q 1 (CNF) Core flow of CNF supernatant used for further studies. The dry content of the dispersion was determined by gravimetric analysis.

Definitions of particle alignment
The general flow direction in the experiment is defined to be the z-direction. To describe the orientation of a single fibril, we will use the polar angle φ as the angle between the fibril major axis and the flow (z) direction while the azimuthal angle θ is the projected angle in the plane perpendicular to the flow (see fig. 1a). The probability of a single fibril having a certain orientation φ and θ is given by the orientation distribution functions (ODF) Ψ φ and Ψ θ , respectively. To describe the average alignment of the CNF in the flow direction, the order parameter S φ = 1 2 3 cos 2 φ − 1 will be used. The brackets in this expression denote an ensemble average over all fibrils.

Channel geometries
Two different channel geometries were used in this study.  The second geometry is a converging channel (CC) as illustrated in fig. 1c. The channel, with an initial rectangular cross-section 4 × 1 mm 2 , is contracted linearly to a quadratic cross-section of 1 × 1 mm 2 over a length of 10 mm. The flow rate in the CC geometry is denoted as Q.
In both channel geometries, the CNF dispersion is accelerated, leading to an increased alignment of the fibrils. There are however some key differences between the FFC and CC geometries that are illustrated in fig. 2a Firstly, when changing the acceleration in the FFC geometry by changing Q 1 and Q 2 , also the cross-section of the core will change. Using optical techniques to characterize the alignment, e.g. SAXS and POM, there will thus give difficulties in comparing the effect of increasing acceleration since there will also be less amount of CNF in a channel cross-section.
Furthermore, at high accelerations the core flow might fluctuate spatially due to instabilities of the flow, while low accelerations can cause a transition to a flow regime where the core fluid is not detached from the walls. 45 In the CC geometry on the other hand, the acceleration can be increased with Q without changing the amount of CNF in the channel cross-section.
Furthermore, the flow is also more stable at high accelerations.
Secondly, in the FFC geometry, the fibrils in the core are approximately only experiencing a uni-axial extensional flow. 30 An assumption can thus be made about cylindrical symmetry around the centerline, with a uniform distribution of the angle θ (Ψ θ = constant).
Furthermore, Ψ φ can be assumed not to vary over the core cross-section. In this channel it is thus possible to easily obtain the steady state ODF Ψ φ,0 from SAXS experiments using the reconstruction method by Rosén et al. 46 . In the CC geometry, close to the walls, the Due to these differences, the following studies will be done in the two channel geometries: • In the FFC geometry, the influence of the ODF Ψ φ,0 on the rotary diffusion process will be studied at the same conditions as Håkansson et al. 2 , 30 , i.e. the core flow rate is Q 1 = 23.4 ml/h and the sheath flow rate is Q 2 /2 = 13.5 ml/h in each of the two channels. In this geometry, the flow is assumed to only have velocity gradients in the flow direction, i.e. a pure extensional flow.
• In the CC geometry, the influence of acceleration on the rotary diffusion process will be studied by changing the flow rate Q. In this geometry, the combination of both shear and extension will be considered.

X-Ray beam
Detector, for an initial distribution with alignment S φ,0 (z). The rotary diffusion coefficient D r can thus be measured in the POM flow-stop experiment by combining the previous two equations where the steady state intensity level during flow I POM,0 (z) is taken before stop. Further details about the post-processing of the data in the flow-stop experiment is given as supplementary information.

SAXS measurements
The  The resulting order parameter S φ,0 (z) obtained from the POM experiments is presented in fig. 6b, where the error-bars correspond to the standard deviation over 18 separate experiments. Compared to the SAXS data, which also has some scatter between the experiments, we find very good agreement using POM. It is thus sufficient to calibrate the POM experiment with only one reference position in the channel.

Dynamic characterization
The decay of the order parameter S φ at z = 2h as function of time t after stop is illustrated in fig. 7a. Note that the light source used in the experiment is not perfectly stable and typically has small intensity fluctuations around its mean intensity. The oscillations seen in fig. 7a should thus not be considered to be caused by oscillations of birefringence. The order parameter S φ (t) is plotted on a logarithmic scale. Therefore, if the dispersion would be dilute and only affected by Brownian motion, we would expect a straight line according to eq. (13). However, this is clearly not the case and the decay seems to occur on multiple time scales where the alignment is initially decaying very rapidly, while at longer times the decay is slower. To analyze this further, we define two rotary diffusion coefficients at short and long times. The data at 0 s < t < 0.1 s is fitted to an exponential function proportional to exp(−6D r,fast t) and the data between at 2 s < t < 4 s is fitted to a function proportional to exp(−6D r,slow t).   Secondly, we find that D r,fast is different depending on downstream position. Initially before the focusing region, the value is around D r,fast = 0.8 rad 2 s −1 . However, in the acceleration zone between 0 < z < h, this effect quickly becomes stronger with a maximum value of D r,fast = 1.27 rad 2 s −1 at z = 0.75h. Even though the alignment of the CNF increases between h < z < 1.5h (seen in fig. 6), the fast diffusion coefficient decays. Further downstream in the channel there seems to be an equilibrium value similar to the one found upstream before the focusing region. In fig. 7c, each pixel is analyzed separately to see the spatial distribution of the fast decay in the channel. In the acceleration zone, the coefficient D r,fast seems to be the same over the whole cross section. However, further downstream the fastest decay rates  The assumption that D r only depends on the order parameter S φ,0 seems not to hold and it becomes even more evident when plotting the two parameters against each other in fig. 8a. The average velocity w and extension rateε in a cross section of the channel can be estimated from the geometry (neglecting the effect of the shear towards the wall): At a given flow rate Q, there is thus an increase in extension rateε, with a maximum at z = 10h. Since A(z) is constant at z > 10h, the extension rate rapidly decays to zero at the end of the contraction.
In fig. 9a, the steady state alignment of the fibrils in the channel is indicated by the intensity I POM,0 at different flow rates Q. As expected, the highest alignment is found at the end of the contraction before the rapid decay of the extension rate. Furthermore, by increasing the flow rate Q at this position, the alignment increases even more as shown in fig. 9b

Dynamic characterization
The flow stop experiment was performed in the same manner as for the FFC geometry.
However, the intensity I POM was directly used to find D r,fast (with 0 s<t<0.1 s) and D r,slow (at 2 s<t<4 s) according to eq. (3) without converting to an order parameter.
Just as for the FFC geometry, the decay of alignment after stop consists typically of a fast and a slow part described by rotary diffusion coefficients D r,fast and D r,slow . The slow part however was more difficult to measure in regions of low initial alignment, as the intensity quickly comes too close to the background level. Instead, we focus the results on the fast diffusion coefficient D r,fast , which also here varies significantly depending on downstream position. In fig. 10c, the fast diffusion coefficient is obtained at Q = 100 ml/h for each pixel and it is clear that it does not correlate well with the average alignment at the same position in figure 9c. Here, D r,fast actually has a local minimum at the position of highest alignment and highest extension rate (z = 10h) for a given flow rate Q. This seems to be completely different compared to the results for the FFC geometry. Using eq. (6), the coefficient D r,fast is plotted versus the extension rateε for different flow rates Q and positions z in fig. 10a.
The results at the end of the contraction are then compared with the results from the FFC geometry in fig. 10b. Given the same value of Q, it is clear that D r,fast decreases with increasingε. However, for a given position z, the coefficient D r,fast actually increases witḣ ε. At the end of the contraction, at z = 10h (purple markers in fig. 10b), the apparent relationship between D r,fast andε seems to be similar to the one found in the FFC geometry (black markers in fig. 10b).
The reason for the difference between different positions could be related to the average velocity gradients in the cross-stream directions, i.e. the shear rateγ. The shear rate will also naturally increase when moving through the contraction. After the contraction at z > 10h, the shear rate remains high, while the extension rapidly vanishes. Increasing the flow rate Q at z = 13h, the shear rate increases along with D r,fast , which is visible in fig. 10a.
From these results, it seems that D r,fast is even more affected by shear than by extension. At a given Q, when moving through the contraction, the influence from extension becomes higher than from shear and by the end of the contraction, the value of D r,fast is almost what is expected when shear is negligible, i.e. close to the value obtained in the FFC geometry of around D r,fast ≈ 1.3 rad 2 s −1 . Assuming a circular cross-section with the same area as the actual quadratic channel, the average shear rateγ in the fully developed flow downstream of the contraction can be estimated to be (see supplementary information for details): An estimation of the average shear rate at Q = 200 ml/h would thus beγ ≈ 263 s −1 , i.e. almost 15 times higher than the maximum extension rate at the same flow rate. Even though the average shear rate at z = 10h should be of this order of magnitude as well, the regions of high shear typically concentrated close to the walls due to the contraction. Thus, there should be lower shear in the rest of the cross section.
In conclusion, the CC geometry can be used to predict the dynamics of CNF both in pure shear and pure extension. The measured dynamics at the end of the contraction can be assumed to not be influenced by shear and will consequently give a prediction of how the dispersed CNF will behave in the FFC geometry. Measuring further downstream of the contraction at a point where the flow is fully developed, the dynamics can similarly be assumed to not be influenced by extension but instead controlled by shear.

Discussion
Understanding

Theoretical considerations in an ideal system
Assuming only Brownian motion in a dilute dispersion of non-interacting stiff rods, the rotary diffusion coefficient is found through: 49 where k b is the Boltzmann constant, T is the temperature, r p is the particle aspect ratio, η s is the solvent viscosity and a is the half length of the fibril. Using the expected values of k b = 1.38 × 10 −23 m 2 kg s −2 K −1 , T = 293 K, r p = 100, η s = 1 mPa s and l/2 = a = 1 µm, the coefficient is found to be D r,dilute ≈ 2 rad 2 s −1 . In a semi-dilute system of isotropic particles, the particle-particle interactions increase and the individual fibrils can not move as freely as in the dilute case. This hindrance due to the concentration was described with the tube-model by Doi and Edwards 50 leading to the expression: where nl 3 is the amount of fibrils inside a cubic volume with the same side as the fibril length. The constant β depends on the nature of the particle-particle interactions and ranges between 1 and 10 3 . The value must be determined experimentally. For the concentration of fibrils here, assuming mass fraction to be roughly equal to the volumetric fraction, a value of nl 3 ≈ 10 is obtained. Depending on β, it is thus expected to find the rotary diffusion coefficient in the interval D r ∈ [0.02, 20] rad 2 s −1 , which is also the range within which we find the present experimental results. Equations (8) and (9) also highlight the effects of concentration (quadratic dependence) and fibril lengths (cubic dependence), which dramatically influence the value of D r .

Reasons for the multiple time scales
A major factor for the multi-time-scale decay of alignment and the dependence on the velocity gradients could be accredited to the fact the CNF dispersion usually is polydisperse and consists of fibrils with various lengths. This conclusion were drawn in several similar POM flow-stop experiments using a Couette apparatus. 36,37,42 The slow decay will thus correspond to the longest fibrils and the fast decay will correspond to the shortest fibrils that are affected by the velocity gradients. Since the stronger gradient will align shorter fibrils, the fastest decay rate also increases. However, there might be several other factors that also come into play that could also lead to a multiple time scale relaxation process.
With the same tube-model that lead to eq. (9), Doi and Edwards 49 also discussed the influence of particle alignment on the rotary diffusion coefficient. The argument was that there was less hindrance from nearby particles in an aligned system, and that D r would be an increasing function of alignment. Håkansson et al. 30 argued that the electrostatic torque imposed by nearby charged fibrils would also increase as the system aligns, which would result in a D r that increases with alignment. An additional hypothesis is that the fibrils could possibly get stretched from an initial non-straight equilibrium shape, and that the elasticity of the individual fibrils causes them to relax to the equilibrium shape; a process which would present itself as a D r that increases with alignment. Hypothetical inter-fibril bonds with some torsional stiffness 51 created at these concentrations could potentially also act in a similar manner. Given any of these hypothetical effects, it is not surprising that the decay of the order parameter after stop is larger initially when the alignment also is

Differences between shear flow and extensional flow
The results in this work indicate that there might be significant differences in CNF dynamics between a sheared dispersion compared to the extended dispersion. The first observation in the CC geometry is that the extensional flow is much more effective in aligning the system The results of the steady state ODF Ψ φ 2D are shown in fig. 11. Given the same P e, the system is much less aligned in the sheared case and the ODF is heavily skewed towards positive values of φ 2D . The reason is that the fibril orientation at φ 2D = 0 is marginally stable in the shear flow. If Brownian rotary diffusion pushes the fibril to positive values of φ 2D it will return to the original orientation. However, a push towards negative values of φ 2D , will lead to the fibril flipping half a revolution. The result is that the most probable orientation actually is larger than zero. In the extensional flow, the orientation φ 2D = 0 is asymptotically stable and Brownian rotary diffusion will not lead to any flipping of the fibrils. The result is that the fibrils are aligning more effectively in the extensional flow as well as the ODF being symmetric around φ 2D = 0.
The probable skewing of the ODF in the shear layer in the CC geometry means that we do not correctly measure the maximum alignment since the polarization filters are tilted 45 • to the flow direction. Additional experiments using different angles of the polarization filters could possibly reveal if the ODF indeed is skewed close to the visible walls.
The reason for the faster dynamics in the sheared dispersion of CNF could be explained by the spatial distribution of the aligned fibrils. Both shear and extension will cause the fibrils to preferentially align in the flow direction. However, in the extended system, the transverse distance between fibrils is actually closer than in the sheared system as illustrated in fig. 11 and fibrils located transversely do not have any relative motion. This could cause potential short-range interactions, for example caused by mechanical friction or chemical bonds, to occur more frequently and thus result in a slower rotary diffusion process. Some resistive short-range interactions could thus possibly explain the difference in behavior between a sheared and a extended dispersion.
The idea of a rotary diffusion coefficient that is dependent on the strain rate tensor magnitude is not new and has been widely addressed in semi-dilute suspension of (non-Brownian) macroscopic fibers, [52][53][54] where the rotary diffusion process itself is a result of short-range interactions between fibers. In these cases, it is argued that these the rate of these interactions, such as direct fiber-fiber collisions, increases with the magnitude of the strain rate tensor and thus affect the rotary diffusion process. However, as already mentioned by Håkansson et al. 30 , the value of this type of correction to D r is typically much lower than what is observed for CNF.
A final remark is that the fastest decay rates in this study were actually obtained in absence of extension. If the length distribution is truly the cause of the multi-time-scale decay, then it is maybe possible to capture a signature of the full length distribution by just utilizing a channel flow with high shear rate without any contraction. The experiment described here would thus be an easy characterization technique to estimate the length distribution of CNF in the dispersion.

Other explanations
The hypotheses of what is causing the experimental observations rely on the assumption that the flow really stops instantly and that the spatial distribution of fibrils is homogeneous over the channel cross section. Secondly, there could maybe be potential migration mechanisms involved in the system of CNF that result in local variations of the fibril concentration over the cross section. If there for example is a mechanism pushing fibrils towards the center of the channel, it will lead to a lower concentration of fibrils closer to the walls. This in turn might lead to a lower POM signal in the shear layers and potentially higher mobility of the fibrils. To analyze this further, the local concentration variations over the cross-section should be carefully examined.
Finally, when the dispersion is subject to shear, the resulting ODF will probably be slightly bi-axial. This means that the alignment will depend on the viewing angle in the plane perpendicular to the flow direction. It might also be possible that the diffusion process of this ODF will seem different depending on the viewing angle. In the converging channel geometry, the diffusion process in the shear layers is observed from different directions depending on if we look close to the lower/upper walls (viewing direction is along the vorticity direction) or if we look in the middle of the channel (viewing direction is along the gradient direction).
Consequently • Previous indications that D r of CNF would depend on the initial ODF are found to be an artefact of other causes since the same ODF can give two different values of D r,fast .
• The coefficient D r,fast is rather dependent on the local velocity gradients both in the streamwise direction (extension) and in cross-stream directions (shear). Higher extension rateε and higher shear rateγ both result in a higher value of D r,fast . High shear in the absence of extension gives a higher value of D r,fast than what is achieved in a pure extensional flow. Thus, a system of fibrils aligned with shear flow will de-align quicker than a system aligned by extensional flow. We hypothesize that there could be resistive short-range interactions that are more likely to occur in the extensional flow and thus slowing down the rotary diffusion process. Since the initial ODF does not seem to influence D r,fast , it is believed that electrostatic fibril-fibril interactions in the dispersion due to the surface charge of the fibrils do not play a major role in the rotary diffusion process.
In order to understand all the underlying processes, we propose future experiments using CNF dispersions with well defined fibril sizes and shapes and performing parametric variations of the concentration, the fibril charge, the pH of the solvent, the temperature or even adding different functional groups to the cellulose molecule.
Additionally, it would be interesting to use computational fluid dynamics (CFD) to determine the actual average shear ratesγ and extension ratesε at the different positions in the channels. Combining experiments and simulations in this way could also increase our knowledge of how rotary diffusion of CNF is affected by velocity gradients.
Even though the rotary diffusion process of CNF is not fully understood, the flow-stop experiment presented here is believed to be a valuable tool for dynamic characterization of CNF dispersions used for material production. In order to create a material with highly aligned fibrils, it would be important to try and minimize the fastest de-alignment time scales in the dispersion. The measurement with this device can consequently be used for a quick determination of the quality and usability of the nanocellulose prior to the material production process.
As a final remark, the presented flow-stop experiment is not limited to the study of cellulosic materials, but can be used to study any dispersion of non-spherical nanoparticles/macromolecules that become birefringent when aligned, e.g. liquid crystals, proteins or polymers.
where d is the distance that the light has traveled through the sample, λ is the wavelength of the light and I 0 is the intensity of the incoming light. This expression holds as long as I I 0 and the angle that the sample has rotated the light is much less than 90 • , which is assumed to be the case in the present experiments. Furthermore, the difference in refractive indices ∆n is related to the order parameter S φ through 47 : where ∆n max corresponds to the difference in refractive index at perfect alignment (S φ = 1).
The order parameter can consequently be found in POM experiments using: Rotary diffusion coefficients can thus be obtained through the decay of the intensity signal I POM (t) by combining equations 13 and 12 and using the stationary order parameter at position z as reference case: In order to only observe the influence of the light from the CNF, a background subtraction is done using measurements with water flowing through the channel. We also divide the subtracted intensity with the same background intensity, to remove any effects of the laser light not being completely uniform over the field of view, i.e.
I POM,raw (z, t) = I POM,CNF (z, t) − I POM,background (z) The intensity in this case however does not decay to zero, as the optical properties of the equilibrium state of the CNF dispersion seems to be different from water. As this equilibrium level I eq. = I POM,raw (z, t → ∞) is not consistently positive or negative (after the background subtraction), we can not draw any conclusions about its origin. Therefore, the mean intensity during the last second of the experiment (around 10 s after stop) is set as the equilibrium level I eq. and the true intensity is given by I POM = I POM,raw − I eq. . All negative values of the intensity after this subtraction are set to 10 −16 in order to convert the intensity to an order parameter according to eq. 12.
Each POM flow-stop experiment is initialized by starting the syringe pumps. After The steady state intensity during flow I POM,0(z) is taken as the mean intensity in a two second period, 3-5 s after the camera starts recording, just before the flow is stopped. The exact stopping time is found during post-processing as the time when the intensity drops below the standard deviation of the mean intensity.

SAXS experiments
As an X-ray beam is passing through the dispersion of CNF, the fluctuations of the electron distribution (due to the solid particles) gives rise to scattered light 24 . The scattered light intensity contains structural information of the particles and can therefore also be used to  quantify the alignment of the CNF. The present procedure to analyze the SAXS data is similar to Håkansson et al. 30 and the additional considerations to obtain S φ is described in detail by Rosén et al. 46 .
The scattered light intensity is probed for different scattering vectors q with direction given by 24 : where e scat. and e inc. refer to unit vectors of the scattered and incident light, respectively.
The length of the scattering vector is given by 24 : where λ is the wavelength of the X-ray and 2ϕ is the angle between e scat. and e inc. . In small angle scattering experiments, the assumption sin ϕ ≈ ϕ holds and we can probe the scattering intensity on a flat detector, where pixel coordinates x d and z d corresponds to scattering vector components 24 : for the sample-to-detector distance D. To remove any contribution from the solvent, a background subtraction was done with scattering from pure water in the channel. An example of a scattering intensity image is showed in figure 12a. The dead pixels on the detector and pixels behind the beamstop were first needed to be corrected. This was done by using the fact that the scattering pattern is symmetric and should ideally contain the same information at (±q x , ±q z ). A pixel without information in the image was thus corrected by replacing its value by a value from any of the other three "mirror" pixels. The resulting intensity image after this procedure is illustrated in figure 12b. The pixel coordinates are transformed to the cylindrical coordinates q = q 2 x + q 2 z and χ = tan −1 (q x /q z ). We then construct a histogram of I SAXS (q, χ) in the range q ∈ [0.25, 0.5] nm −1 ( fig. 12c) with certain number of bins in qand χ-directions. The number of bins in qand χ-direction was chosen to be significantly lower than the number of pixels in the qand χ-directions, respectively. This is ensuring that each histogram bin contains information.
x z y φ θ χ Figure 13: Orientation of an elongated particle described by the polar angle φ and the azimuthal angle θ; the projected angle observed in the xz-plane is denoted χ.
For each subrange of q, the isotropic contribution to the scattering is removed by assuming that the highest aligned case (at z = 1.5h) has no fibers aligned perpendicular to the flow. This means that all intensity values in the q-subrange are subtracted by I iso (q) ≈ I SAXS (q, χ = ±π/2). Each subrange in q is then normalized, the subranges are averaged and we obtain the steady state orientation distribution function (ODF) of Ψ χ,0 ( fig. 12d). The measurement thus gives the average projected angle χ of particles in the viewing (xz) plane, as illustrated in figure 13. This angle is related to φ and θ through the relationship: tan χ = tan φ cos θ. (19) To obtain the ODF of the polar angle φ, we make the assumption that the azimuthal angle θ is distributed uniformly in the flow. The steady state ODF Ψ φ,0 is then obtained with the same method as described by Rosén et al. 46 . The corresponding steady state order parameter S φ,0 is found through: with normalization Estimating the average shear rate in a fully developed channel flow The velocity w at a radial position r of a fully developed flow in a circular pipe with radius R and streamwise pressure gradient dp/dz is given directly from the incompressible Navier-Stokes equations as: where µ is the dynamic viscosity of the fluid. The flow rate Q can be related to the pressure gradient through: Q = 2π 0 R 0 w(r)rdrdφ = 2π 4µ dp dz R 0 (rR 2 − r 3 )dr = π 2µ dp dz The average shear rateγ avg.,pipe in the pipe will be given by: γ avg.,pipe = 1 πR 2 2π 0 R 0 dw dr rdrdφ = 2 R 2 R 0 1 2µ dp dz r 2 dr = R 3µ dp dz .
Using the result in eq. (23), this can be rewritten to: Assuming the same average shear rate in a quadratic channel with side h and same crosssectional area, i.e. R = h/ √ π, we finally obtain: The two-dimensional Smoluchowski equation The steady state two-dimensional ODF Ψ φ 2D is given by the solution of the stationary 2D Smoluchowski equation: 49 In the simulations, the particles are assumed to be prolate spheroids with aspect ratio (length/width) r p = 100.
In a shear flow, the Péclet number is given by P e =γ/D r , whereγ is the shear rate and D r is the rotary diffusion coefficient. The angular velocityφ 2D is given by Jeffery 15 as: In a planar extensional flow, the Péclet number is given by P e =ε/D r , whereε is the extension rate. In this case the angular velocity is given by: 15