Colloidal rod dynamics under large amplitude oscillatory extensional flow

Abstract

We perform a combined experimental and theoretical investigation of the orientational dynamics of colloidal rods subjected to time-dependent homogeneous planar elongational flow. Our experimental approach involves the flow of dilute suspensions of cellulose nanocrystals (CNC) within a cross-slot-type microfluidic device through which the extension rate is modulated sinusoidally over a wide range of Péclet number amplitudes (Pe₀) and Deborah numbers (De). The time-dependent orientation of the CNC is assessed via quantitative flow-induced birefringence measurements. For small Pe₀ ≲ 1 and small De ≲ 0.03, the birefringence response is sinusoidal and in phase with the strain rate. With increasing Pe₀, the response becomes non-sinusoidal (i.e., nonlinear) as the birefringence saturates due to the high degree of particle alignment at higher strain rates during the cycle. With increasing De, the CNC rods have insufficient time to respond to the rapidly changing strain rate, leading to asymmetry in the birefringence response around the minima and a residual effect as the strain rate approaches zero. These dynamical responses of the CNC are captured in a detailed series of Lissajous plots of the birefringence versus the strain rate. Experimental measurements are compared with simulations performed on both monodisperse and polydisperse systems, with rotational diffusion coefficients D_r matched to the CNC. A semiquantitative agreement is found for polydisperse simulations with D_r heavily weighted to the longest rods in the measured CNC distribution. The results will be valuable for understanding, predicting and optimizing the orientation of rod-like colloids during transient processing flows such as fiber spinning and film casting.

---

1 Introduction

Suspensions of anisotropic particles are prevalent in a wide range of applications, including industrial processing, biomedical engineering, drug delivery, and the development of advanced soft materials. The macroscopic properties of these systems, such as mechanical strength, optical anisotropy, electrical conductivity, and magnetic responsiveness, are often governed by the degree and nature of particle alignment.^1,2 To tailor these properties for specific applications, it is essential to control the orientation of anisotropic particles during flow. Consequently, understanding the hydrodynamic alignment behavior of such particles is of central importance in the design and optimization of processing operations.

While industrial processes typically involve complex, mixed, and transient flows, much of the previous research has focused on the dynamics of colloidal rods under simpler conditions, such as steady shear or extensional flows.^3–8 It has been experimentally shown that the extent of rod alignment (in the absence of inertia) depends on the Péclet number Pe, a non-dimensional parameter that captures the ratio between the deformation rate |E| and the rotational Brownian diffusion coefficient D_r, i.e., Pe = |E|/D_r. Three general trends are observed as a function of Pe. For Pe ≲ 1, the deformation rate is insufficient to overcome rotational Brownian diffusion, leading to rods that remain isotropically oriented, as in equilibrium conditions. At Pe ≳ 1, rods begin to spend an increasing fraction of their time in a preferentially oriented configuration, and as Pe → ∞, they approach a limiting and maximal degree of alignment.^5,9 Experimentally, the extent of rod alignment as a function of Pe has typically been probed using birefringence or small-angle scattering techniques, in combination with torsional rheometers to control shearing flow, or microfluidic devices to allow precise control over both flow type and strength.^{3–6,10–15}

In recent years, considerable interest has been devoted to understanding the dynamics of colloidal rods in complex flow environments that more closely resemble real-life conditions. For example, flow-focusing devices, contraction–expansion geometries, and fluidic four-roll mills have been employed to mimic specific flow fields of interest in a controlled manner.^{8,12–14,16,17} These devices share the common feature of generating regions with some degree of extensional flow, which is essential for replicating deformation modes relevant to many practical applications. In such systems, the flow is steady in time (i.e., in the Eulerian frame), whereas particles experience a transient flow history, and consequently a varying deformation rate, as they are advected through the region of interest (i.e., in the Lagrangian frame). Moreover, colloidal rods in these environments often undergo transient variations in flow type (e.g., during transitions from shear- to extension-dominated regions), which adds significant complexity, making data interpretation challenging and sometimes obscuring the key underlying physical mechanisms.

To overcome such difficulties while still capturing structural information relevant to real-life scenarios, one strategy is to study the dynamics of colloidal rods exposed to a constant type of deformation (or flow type) while transiently varying the magnitude of the deformation rate. In rheology, this is typically done with oscillatory shear tests. In particular, small-amplitude oscillatory shear (SAOS) has been established as a standard method for characterizing linear viscoelastic behavior under infinitesimal deformations, such that the material microstructure remains practically unperturbed from its equilibrium state.¹⁸ However, with regard to measurements of structural evolution, SAOS provides little to no difference compared to equilibrium conditions, as a SAOS measurement, by definition, preserves the material's equilibrium microstructure. In contrast, large-amplitude oscillatory shear (LAOS) is typically employed to probe the material response far from equilibrium, hence capturing features of viscoelasticity and microstructural dynamics under nonlinear deformations.^1,19–23 Considerable attention has been devoted to understanding the connection between microstructure and the macroscopic response of materials under LAOS using scattering techniques such as time-resolved small-angle neutron scattering (t-SANS).^24–27 Typically, such experiments involve the simultaneous measurement of the evolving stress under a sinusoidally varying deformation rate, together with a concurrent determination of the microstructural orientation from the scattering pattern. Using this approach, it has been possible, for instance, to identify specific stress–orientation relationships linking the degree of structural order to the stress as a function of oscillation frequency.²⁵

The presence of a rotational component in LAOS and shear experiments in general leads to dynamics that differ markedly from those observed in extensional flows. For instance, in shear flows, colloidal rods exhibit tumbling and wagging motions, types of dynamics that do not occur in purely extensional flows.^28,29 Additionally, in shear, rods adopt a preferential orientation angle of θ ≈ 45° relative to the flow direction at Pe ≳ 1, which progressively decreases toward θ → 0° as Pe → ∞. In contrast, in extensional flows, θ ≈ 0° remains constant for Pe ≳ 1.^5,28,30–33

Perhaps even more striking is the case of flexible polymers, for which extensional flows induce a drastic conformational transition, from a coil-like state to a flow-sustained stretched configuration (the so-called coil-to-stretch transition), that is much less pronounced in shear flow.^6,34–36 This highlights the need to complement structural analyses performed under LAOS with its counterpart, large-amplitude oscillatory extension (LAOE), to investigate polymer and colloidal dynamics under transient extensional flows.^37–39

In our previous work, we demonstrated that LAOE under planar extension can be generated in an optimized microfluidic cross-slot, focusing on the retrieval of pressure drop and on the modification of the flow field induced by stretching polymers.³⁹ In this study, we extend the use of our LAOE setup to investigate the structural evolution of rigid colloidal rods, specifically cellulose nanocrystals (CNC), using birefringence techniques to probe the extent of rod alignment. Rigid rods, such as the CNC employed here, not only serve as suitable building blocks for engineering well-defined, ordered soft bio-materials but also provide a useful model system for studying inelastic fluids governed solely by particle orientation. This stands in contrast to flexible polymers, which undergo complex flow-induced conformational changes, such as coil-to-stretch and stretch-to-coil transitions, under extensional flows. Nevertheless, the dynamic behavior of rod-like colloidal systems in time-dependent or oscillatory extensional flows remains largely unexplored.⁴⁰

In this work, we employ the LAOE framework to study the dynamics of CNC as a model system for rigid colloidal rods. We demonstrate the usefulness of LAOE for investigating colloidal rod dynamics in transient extensional flows and provide a direct comparison with simulations to validate and support our experimental methodology. We reveal how the dynamical response of the CNC under LAOE depends on the amplitude and period of the applied sinusoidal extension rate, becoming strongly nonlinear at high amplitudes and short periods. We anticipate the results to be valuable for the optimization of particle alignment in processing applications involving transient mixed shear and extensional flows.

2 Material and methods

2.1 Test fluids and rheological characterization

Cellulose nanocrystals (CNCs; CelluForce, Montreal, Canada) were dispersed at volume concentrations (vol/vol) of 0.05%, 0.1%, and 0.2% in a glycerol–water mixture containing 80 vol% glycerol. This solvent mixture behaves as a Newtonian fluid with a shear viscosity of (η_s = 74 mPa s). Fig. 1(a) shows the size distribution of the CNC rods obtained from atomic force microscopy measurements reported by Calabrese et al.³² The rods are polydisperse, with a number-averaged contour length of 〈l_c〉 = 260 nm, which is significantly smaller than the volume-averaged contour length, . An average diameter of 〈d〉 = 4.8 nm and a persistence length of l_p = 30〈l_c〉 were also estimated by Calabrese et al.³²

Fig. 1 CNC characterization. (a) Histogram of the characteristic contour length l_c distribution using N = 10 bins, as obtained from AFM counting of 976 isolated particles (see ref. 32). (b) Distribution of the rotational diffusion coefficient D_r, calculated from the data in (a) using eqn (1), with l_c corresponding to that of each bin.

Within the Doi–Edwards framework, concentration regimes can be identified based on the rod number density ν (i.e., the number of rods per unit volume) and the rod length l_c. Very generally, for νl_c³ ≲ 1, rods are isolated and effectively non-interacting, whereas for νl_c³ > 1, inter-rod interactions become relevant. Using 〈l_c〉 = 260 nm and a reported cellulose density of ρ = 1500 kg m⁻³,⁴¹ we estimate 1.2 ≲ ν〈l_c〉³ ≲ 5 for the concentration range explored here (i.e., 0.05%, 0.1%, and 0.2%). While these estimates suggest the possible presence of inter-rod interactions, it is important to note that they do not account for polydispersity and should therefore be interpreted as indicative. In previous studies, Bertsch et al., using CNCs from the same manufacturer, reported the onset of interactions at concentrations of ≈0.5% based on scattering experiments.^42,43 Furthermore, Calabrese et al., using the same stock CNC dispersion as employed here, identified the onset of interactions, defined as the concentration at which the diffusion coefficient becomes concentration dependent, at approximately 0.3%.⁴ Based on these previous results, and further supported by our steady-state birefringence experiments discussed in Section 3.1, we consider the CNCs at the concentrations used here to be dilute and effectively non-interacting.

For non-interacting rods, the rotational diffusion coefficient is expected to be concentration independent and is given by:^9,44

where k_B is the Boltzmann constant, T = 295 K is the absolute temperature, η_s is the solvent shear viscosity, and d_eff is the effective diameter of the rods, which accounts for the bare rod diameter d and the contribution of the electric double layer δd as d_eff = d + δd.⁴³ We consider δd = 22.6 nm, as reported in deionized water, to be a representative value of the electric double layer thickness.⁴³ Based on the number-averaged length (〈l_c〉 = 260 nm), we estimate D_r ≈ 6.73 s⁻¹ while based on the volume-averaged length , we obtain D̄_r ≈ 1.44 s⁻¹. The distribution of D_r (Fig. 1(b)) used in the polydisperse simulations to be described later (Section 2.6) is derived from the contour length distribution shown in Fig. 1(a).

The steady shear rheology of the test fluids was measured at 22 °C using a stress-controlled rotational rheometer (MCR502, Anton Paar, Austria) equipped with a double gap Couette geometry. The 0.05% and 0.1% samples exhibit a nearly constant viscosity close to that of the Newtonian solvent (labeled as 0%, see Fig. 2(a)). At 0.2%, the sample shows slight shear-thinning behavior within the investigated shear rate range. Average viscosity values for the 0.05%, 0.1%, and 0.2% CNC suspensions in the shear rate range of 10 ≤ [small gamma, Greek, dot above] ≤ 1000 s⁻¹ are η = 79 mPa s, η = 82 mPa s, and η = 88 mPa s, respectively.

Fig. 2 Rheological responses of the CNC test fluids in shear and uniaxial extension. (a) Shear viscosity as a function of the applied shear rate. (b) Representative curves of the thinning filament diameter during capillary thinning in a CaBER experiment. The legend in (a) indicates the CNC concentration in vol%. The inset in (b) displays a representative snapshot of the filament for the 0.1% CNC suspension, where the white scale bar represents 2 mm.

The behavior of the test fluids under uniaxial extension was characterized using a capillary breakup extensional rheometer (CaBER) device (Thermo-Haake, Germany). Circular plates with a diameter of 6 mm were separated from an initial gap of 2 mm to a final gap of 6 mm within a strike time of 20 ms. The diameter D(t) of the resulting fluid filament was monitored as it thinned over time (see Fig. 2(b), where t = 0 s corresponds to the time at which the plates reached their final separation). We did not observe a detectable exponential filament decay characteristic of the elasto-capillary thinning regime, which is consistent with the rigid CNC rods being effectively inextensible and having negligible entropic elasticity.^45–47

2.2 Microfluidic setup

A microfluidic optimized shape cross-slot extensional rheometer (OSCER) geometry was used to generate a planar extensional flow.^48,49 The OSCER has a half-height of H = 1 mm and a characteristic half-width of W = 100 µm for the inlet and outlet channels (Fig. 3(a)). The high aspect ratio of the geometry H/W = 10 ensures a uniform flow across most of the channel's height and hence a good approximation to a two-dimensional (2D) planar flow field. The flow in the central part of the OSCER geometry (|x/W|, |y/W| ≤ 15) approximates pure planar homogeneous elongation with hyperbolic streamlines surrounding a free stagnation point located at x = y = 0.^50–52 The geometry is fabricated from stainless steel using wire electrical discharge machining and is sealed with glass slides on both the top and bottom surfaces, enabling optical interrogation of the flow in the region surrounding the stagnation point.

Fig. 3 Overview of the experimental setup. (a) Schematic illustration of the optimized shape cross-slot extensional rheometer (OSCER) geometry. The geometry has a height of 2H and features two pairs of opposing inlet and outlet channels aligned along the x and y axes, respectively, each with a width of 2W. (b) Representation of the normalized set extension rate profiles [small epsi, Greek, dot above]_x,set/[small epsi, Greek, dot above]_0,set and [small epsi, Greek, dot above]_y,set/[small epsi, Greek, dot above]_0,set in the x and y directions, respectively, shown over one normalized period t/T.

2.3 Flow control

2.3.1 Steady flow conditions. The fluids are driven through the microfluidic channel using low-pressure syringe pumps (Nemesys S, 29 : 1 gear ratio, Cetoni GmbH, Germany), equipped with proprietary 5 mL borosilicate glass syringes. The syringes are connected to the OSCER device via PTFE tubing (inner diameter 1 mm, Darwin Microfluidics, France). Careful measures are implemented to eliminate air bubbles in the microfluidic system, which can impact the response of the system to the transient flows imposed during LAOE.³⁹

Two of the pumps are set to inject the fluid at equal volumetric flow rates Q into the opposing inlet channels aligned along the y direction, and two additional pumps are set to simultaneously withdraw fluid at equal and opposite rates from the outlet channels along the x direction. The average flow velocity in each channel of the geometry is U = Q/(4WH) and the resulting set extension rate is [small epsi, Greek, dot above]_set = 0.1U/W.^49,50,52,53 Under this arrangement, the x axis corresponds to the extensional axis and the y axis is the compressional axis, and we denote the average strain rates along the extensional and compressional axes as [small epsi, Greek, dot above]_x and [small epsi, Greek, dot above]_y, respectively. For a Newtonian fluid under steady 2D planar elongation, [small epsi, Greek, dot above]_x = −[small epsi, Greek, dot above]_y = [small epsi, Greek, dot above]_set. However, for transient LAOE flows in the OSCER device, system damping can result in a phase shift and reduced amplitude of [small epsi, Greek, dot above]_x and [small epsi, Greek, dot above]_y relative to [small epsi, Greek, dot above]_set.³⁹

2.3.2 Time-dependent flow conditions. In this study, we investigate the response of our dilute CNC suspensions to pulsatile LAOE.³⁹ For this, we impose sinusoidal extension rate profiles at the inlets and outlets of the microfluidic OSCER device described by:

[small epsi, Greek, dot above]_set(t) = [small epsi, Greek, dot above]_x,set(t) = −[small epsi, Greek, dot above]_y,set(t) = [small epsi, Greek, dot above]_0,set[1 + sin(2πt/T)], (2)

where [small epsi, Greek, dot above]_0,set is the amplitude of the extension rate oscillation, and T = 1/f is the pulsation period with f as the pulsation frequency. In this study, we cover a broad range of 0.1 ≤ [small epsi, Greek, dot above]_0,set ≤ 100 s⁻¹ and 1 ≤ T ≤ 50 s. Note that we impose a unidirectional pulsatile flow, where the extension rates oscillate with an amplitude of [small epsi, Greek, dot above]_0,set around constant background flow of [small epsi, Greek, dot above]_0,set. The resulting extension rate profiles (as commanded to the pumps) along the x and y axes of the OSCER device are shown in Fig. 3(b). The constant background flow is maintained for several seconds to stabilize the flow before the superimposed oscillatory modulation begins. We use a multi-function DAQ device (USB-6009, National Instruments, TX) to synchronize the pumps and velocimetry system or polarization camera with a global trigger signal at the beginning of the pulsation cycle.

2.4 Microparticle image velocimetry

Applying time-dependent flows in microfluidic systems can lead to significant deviations between the set input signal ([small epsi, Greek, dot above]_set(t)) and the actual extension rate profile within the microfluidic chip.^39,54 Therefore, for each imposed oscillatory extension rate profile, we measure the flow field in the OSCER using micro-particle image velocimetry (μ-PIV, TSI Inc., MN).^55,56 For this, the OSCER device is mounted on an inverted microscope (Eclipse Ti, Nikon, NY) equipped with a 4× air objective (PlanFluor, Nikon, NY) with a numerical aperture of NA = 0.13. The fluids are seeded with 0.02 wt% of 2 µm red fluorescent tracer particles (Fluor-Max, Thermo Scientific, Germany) with excitation and emission wavelengths of 542 nm and 612 nm, respectively. The tracer particles are excited with a dual-pulsed Nd:YLF laser (527 nm) using a volumetric illumination technique. The relative depth over which the tracer particles contribute to the velocity field measurement is δ_z ≈ 0.16H.⁵⁷ The flow in the xy midplane of the OSCER geometry is recorded using a high-speed camera (Phantom MIRO, Vision Research, Canada), which operates in frame-straddling mode, synchronized with the laser. The frame rate of the camera and the separation of laser pulses are adjusted based on the applied flow rate to achieve an average particle displacement of approximately 4 pixels between consecutive images.

For steady flow measurements, 50 image pairs are recorded and ensemble-averaged over the image sequence. For measurements under pulsatile LAOE, the camera frame rate is set to capture at least 125 image pairs per cycle. At least two full oscillation cycles are recorded per PIV measurement. PIV analysis (TSI Insight 4G, TSI Inc., MN) is performed to obtain the velocity components u and v in the x and y directions, respectively. Subsequent analysis uses a custom MATLAB (R2024a, The MathWorks, MA) algorithm.

We fit the extension rate [small epsi, Greek, dot above]_x(t) along the extension axis measured by PIV with a sinusoidal function, [small epsi, Greek, dot above]_fit(t) = [small epsi, Greek, dot above]_off + [small epsi, Greek, dot above]₀ sin(2πt/T + φ). From this fit, the offset [small epsi, Greek, dot above]_off, the amplitude [small epsi, Greek, dot above]₀ of the measured signal, and the phase shift φ between the measured strain rate profile and the set signal are extracted. Hence, the maximum of the temporal strain rate profile is [small epsi, Greek, dot above]_max = [small epsi, Greek, dot above]₀ + [small epsi, Greek, dot above]_off, which is used for data normalization.

2.5 Birefringence imaging

Birefringence imaging was performed following the methodology described by Onuma and Otani.⁵⁸ A high-speed polarization camera (CRYSTA PI-1P, Photron Ltd, Japan) is fitted with a 10× air objective (PlanFluor, Nikon, NY) with a numerical aperture of NA = 0.3, and monochromatic light with a wavelength of λ ≈ 520 nm. We measure the retardance R to probe the extent of fluid anisotropy due to the rod alignment during flow and convert it to birefringence as Δn = R/(2H). A background value of R ≈ 0.5 nm was determined for the test fluid at rest as the lower detection limit of the device. This noise level is subtracted from all experimental retardation data. The upper detection limit of the device depends directly on the wavelength as R = λ/4 ≈ 130 nm, and is not exceeded in our experiments.⁵⁸

Under time-dependent flow conditions, a frame rate between 60–250 fps is used based on the set extension rate profile to capture at least 250 images per cycle during at least two full oscillation cycles. This data is related to the measured strain rate profile determined via the PIV. We spatially average the birefringence 〈Δn〉 and the orientation angle 〈θ〉 around the central stagnation point along the stretching axis (|x| ≤ 0.2 mm and |y| ≤ 0.02 mm) in each image acquired by the polarization camera.

2.6 Simulations

We consider a dilute suspension of inertialess rods with a finite aspect ratio, homogeneously dispersed in a Newtonian medium. The adjective dilute implies that each particle is hydrodynamically independent. Each particle has an aspect ratio r and its orientation is defined by a unit vector p parallel to the axis of revolution of the spheroid. In these conditions, the time change of the single spheroid orientation follows the Jeffery's equation:⁵⁹

ṗ = Ω·p + λ(E·p − E: ppp), (3)

where E is the rate of strain tensor and Ω the vorticity tensor. λ represents the spheroid form factor, and it is defined as the ratio (r² − 1/r² + 1). To describe the evolution of the orientation of a suspension, it is necessary to employ the orientation distribution function ψ(p,t) defined such that the probability to find a particle oriented within a solid angle dp of p is ψdp. The time evolution of ψ is governed by the following equation:⁶⁰The orientation distribution evolves because of the torque induced by the external flow field and Brownian motion. From the orientation distribution, it is possible to calculate the second-order orientation tensor defined as:This tensor contains average information on the overall orientation state of the system. Following Fuller,⁶¹ we can express the average anisotropy factor with respect to the flow direction as a function of the orientation tensor as:where a_ij are the components of a₂. We solve numerically eqn (4) since no simple analytical solution is possible. We employ a finite volume method as previously performed by Férec et al.^62,63 (see Férec et al.⁶² for more details about the mesh and the treatment of the periodic boundary conditions). A central scheme is implemented to interpolate properties between nodes. Once the orientation distribution is numerically solved, the orientation tensor is evaluated, and the anisotropy factor is consequently obtained.

In the model, two input parameters are required: λ and D_r. These two parameters can be obtained experimentally from AFM particle characterization (Fig. 1). For a monodisperse system, the average values of λ and D_r (or D̄_r) are obtained using the CNC number- and volume-weighted average lengths.

In order to account for polydispersity, we defined subsets of the suspensions and calculated (for each of the subsets) the average D̄_r starting from the AFM characterization as shown in Fig. 1. Since the system is dilute, each subset of particles can be simulated independently from the rest of the subsets. Hence, the N_i particles in the subset i have an average rotational diffusion D̄_r,i (calculated according to eqn (1)), an average form factor [small lambda, Greek, macron]_i and average length L̄_i. By solving eqn (4) with D̄_r,i and [small lambda, Greek, macron]_i, we can calculate the evolution of average AF for the subset i, AF_i. Once all the subsets are solved, we calculate the following averages:

These averages consider the impact on birefringence given by the number of particles A₀, the size of the particle A₁, the area that the particles span A₂, and the volume that the particles A₃ cover during motion. A similar approach was employed to explain the effect of polydispersity in graphene suspensions under simple shear flow, where it was found that the dichroic response was strongly influenced by the area spanned by each particle during rotation.³⁰ It is worth noting that an alterative averaging could be employed here. One can average the microstructural state of the system by averaging, with equivalent weights used in eqn (7)–(10), the components of a₂.⁶⁴ For our dataset, the difference between the two approaches is less than the experimental error (<3%).

3 Results and discussion

First, we probe the flow behavior and birefringence of the CNC test fluids under steady flow conditions in Section 3.1 to validate the homogeneity of the flow field with the covered strain rate range. The experimental results are compared to simulations of monodisperse and polydisperse systems. Second, we employ pulsatile flow conditions and test the dynamic behavior of the rods under LAOE in Section 3.2, focusing on the temporal evolution of the birefringence.

Note that in all experiments, we define the Reynolds number, describing the ratio of inertial to viscous forces during flow, as Re = ρUD_h/η, where the fluid density, ρ = 1215 kg m⁻³, the hydraulic diameter of the rectangular inlet and outlet channels D_h = 4WH/(W + H) = 0.364 mm, and fluid viscosity η. Based on the maximum velocity reached in the experiments, the maximum Reynolds number reached in this study is Re ≈ 0.23.

3.1 Steady flow conditions

3.1.1 Flow field characterization. Under conditions of steady extensional flow, the flow field measured for all of the test fluids remains Newtonian-like over the entire range of investigated strain rates (i.e., 0.1 ≤ [small epsi, Greek, dot above]_set ≤ 100 s⁻¹, Re < 0.23), maintaining symmetry about the central stagnation point and the principal flow axes. This is exemplified in Fig. 4(a) for the 0.1% CNC sample at two representative strain rates. The velocity components v (in y direction) and u (in x direction) are extracted along their respective axes, indicated by the dashed lines in Fig. 4(a). Over the range of flow rates studied, u increases, and v decreases linearly along their respective axes (representatively shown for [small epsi, Greek, dot above]_set = 100 s⁻¹ in Fig. 4(b)). The average strain rate in the x direction ([small epsi, Greek, dot above]_x = ∂u/∂x) and the y direction ([small epsi, Greek, dot above]_y = ∂v/∂y) are obtained via μ-PIV by averaging the velocity gradients over the spatial domain |y/W| ≤ 6 and |x/W| ≤ 6 (indicated by the black dashed lines in Fig. 4(b)).

Fig. 4 Flow velocimetry obtained with the 0.1% CNC dispersion under steady flow conditions. (a) Normalized velocity field with superimposed streamlines at [small epsi, Greek, dot above]_set = 1 s⁻¹ (left) and [small epsi, Greek, dot above]_set = 100 s⁻¹ (right). (b) Normalized velocity components along the x and y directions representatively shown for [small epsi, Greek, dot above]_set = 100 s⁻¹. Black dashed lines represent linear fits over |y/W| ≤ 6 and |x/W| ≤ 6, which are used to calculate the extension rates [small epsi, Greek, dot above]_x and [small epsi, Greek, dot above]_y, respectively. (c) Extension rates along the x and y directions as a function of the set elongation rate. (d) Streamwise velocity profiles normalized by the centerline velocity, measured across an outlet 7 mm downstream of the stagnation point at various imposed [small epsi, Greek, dot above]_set.

The magnitudes of the average extension rates along the compression and elongation axes are the same (i.e., [small epsi, Greek, dot above] ≡ |[small epsi, Greek, dot above]_x| = |[small epsi, Greek, dot above]_y|) and increase linearly as a function of the set extension rate [small epsi, Greek, dot above]_set (Fig. 4(c)).

Across all investigated strain rates, we observe parabolic profiles of the streamwise velocity across the outlet and inlet channels (Fig. 4(d)) without the emergence of any flow modifications, which have been reported for flows of polymeric systems even at c < c*.^39,50,53,65 Taken together, the fact that [small epsi, Greek, dot above]_x = −[small epsi, Greek, dot above]_y = [small epsi, Greek, dot above]_set and the absence of flow modification at both the outlet and inlet imply that the CNC dispersions behave analogously to a Newtonian fluid, consistent with previous observations for rigid colloidal rods at c ≲ c*.^6,32

Based on the homogeneity of the flow field during planar extension for the used CNC systems in this study, we define [small epsi, Greek, dot above] ≡ [small epsi, Greek, dot above]_x and [small epsi, Greek, dot above](t) ≡ [small epsi, Greek, dot above]_x(t) as characteristic strain rate imposed on the fluid along the stretching axis.

3.1.2 Birefringence. Under steady flow conditions, we observe the emergence of a relatively large area of birefringence around the central stagnation point in the OSCER device as the extension rate is increased (Fig. 5(a)). This relatively large region of strong birefringence around the stagnation point can be attributed to the small fluid strain required for the rods to attain a preferential orientation. This behavior contrasts sharply with that of flexible polymer solutions, which exhibit a highly localized birefringent strand along the extension axis, since polymer chains must experience a substantially larger fluid strain before undergoing stretching and backbone alignment.

Fig. 5 Birefringence of the CNC dispersions under steady flow conditions. (a) Time-averaged birefringence normalized by the CNC concentration for the same conditions as in Fig. 4(a). The orientation angle θ is indicated by the solid segments. The white boxes around the central stagnation point in (a) correspond to regions of |x| ≤ 0.2 mm and |y| ≤ 0.02 mm, where the spatially-averaged data is calculated. (b) Spatially-averaged orientation angle 〈θ〉 (top) and birefringence 〈Δn〉 (bottom) for the three samples as a function of the extension rate [small epsi, Greek, dot above]. Gray areas in (b) indicate the polarization camera's maximum retardance limit. (c) Spatial averaged birefringence 〈Δn〉 scaled with the CNC volume fraction ϕ as a function of Péclet number Pe, using D̄_r = 1.44 s⁻¹. The red solid line shows a fit to the data according to eqn (11).

In the microfluidic OSCER device, the CNC aligns perpendicularly to the flow direction along the compressional y-axis due to the deceleration of the fluid, as indicated by the solid black segments in Fig. 5(a). In contrast, the CNC aligns parallel to the flow direction along the elongation axis in the x direction. In previous studies, similar orientation patterns have been reported for CNC systems and other anisotropic particles.^12,32,66,67 We spatially average the orientation angle 〈θ〉 in the central part of the OSCER around the stagnation point (white boxes in Fig. 5(a)) to quantify the CNC alignment along the extension axis. We observe alignment in the flow direction with 〈θ〉 ≈ 0° over the entire range of tested [small epsi, Greek, dot above] (Fig. 5(b), top), consistent with a previous report.³²

The average birefringence 〈Δn〉, also spatially-averaged in the central part of the OSCER around the stagnation point (white boxes in Fig. 5(a)), increases with increasing extension rate [small epsi, Greek, dot above] for all investigated samples (Fig. 5(b), bottom). At low strain rates, the birefringence increases almost linearly with a slope of 1 for [small epsi, Greek, dot above] ≲ 1 s⁻¹. Increasing the strain rate results in a progressive saturation of the birefringence. For [small epsi, Greek, dot above] ≳ 10 s⁻¹, the birefringence approaches a plateau for 0.05% and 0.1% CNC (Fig. 5(b), bottom). For the 0.2% CNC suspension, the birefringence at the saturation plateau lies beyond the maximum retardance limit of the polarization camera.

Since the birefringence intensity is proportional to the number of aligned rods in the optical path, we normalize the averaged birefringence data 〈Δn〉 of Fig. 5(b) by the corresponding CNC volume fraction ϕ of the samples (ϕ = vol%/100). The resulting 〈Δn〉/ϕ curves are plotted as a function of the Péclet number, computed using the rotational diffusion coefficient for dilute rods as Pe = [small epsi, Greek, dot above]/D̄_r, with D̄_r = 1.44 s⁻¹ (Fig. 5(c)). The collapse of the birefringence curves onto a single master curve indicates that the use of the rotational diffusion coefficient for dilute suspensions (given by eqn (1)) is appropriate. If inter-rod interactions were significant, the effective rotational diffusion coefficient would decrease with increasing concentration, and a concentration-dependent rotational diffusion coefficient would be required to achieve a comparable collapse. This behavior has been demonstrated for colloidal rods in general, and specifically for the CNCs used here.^3,4 Furthermore, the overall collapse of the curves provides clear evidence of self-similar behavior across concentrations, consistent with the rods experiencing the same flow-induced dynamics irrespective of concentration. The data in Fig. 5(c) also highlight the linear increase in birefringence with a slope of 1 for Pe ≲ 1 followed by a progressive saturation at large Pe ≳ 10.

This trend in birefringence as a function of Pe can be described, similarly to Santos et al.,³¹ by the following empirical relation:

By fitting eqn (11) to the experimental data shown in Fig. 5(c), we obtain Δn_max/ϕ ≈ 3.8 × 10⁻² as the plateau value of the normalized birefringence at high Pe, a = 6 as the inflection point of the sigmoidal curve, and b ≈ 1 as the slope parameter. The value of Δn_max/ϕ = 3.8 × 10⁻², often referred to as the intrinsic birefringence, is in good agreement with those reported for cellulose-based structures.⁶⁸

3.1.3 Simulations. Fig. 6 compares the experimental birefringence data 〈Δn〉/ϕ (see Fig. 5(c)) with numerical simulations for monodisperse (Fig. 6(a)) and polydisperse systems (Fig. 6(b)). For the monodisperse case, two simulations were performed using D̄_r = 1.44 s⁻¹ and D_r = 6.73 s⁻¹, based on the volume-averaged contour length and the number-averaged contour length 〈l_c〉 = 260 nm, respectively. The degree of CNC alignment in the simulation is quantified by the anisotropy factor AF, which ranges from 0 to 1. We scale AF by α_s = 1/(Δn_max/ϕ) ≈ 26, i.e., as AF/α_s, to enable a direct comparison with the experimental birefringence data. This scaling is applied uniformly to all numerical data for both monodisperse and polydisperse systems.

Fig. 6 Simulation results for (a) a monodisperse system and (b) polydisperse systems showing the anisotropy factor AF scaled by a constant factor α_s. Experimental data of normalized birefringence are superimposed. In (a), two simulations for different diffusion coefficients are shown corresponding to calculations of eqn (1) using the volume-averaged contour length (D̄_r = 1.44 s⁻¹) and the number-averaged contour length 〈l_c〉 = 260 nm (D_r = 6.73 s⁻¹). (b) Simulations for different averages A_0–3 using N = 10 bins based on Fig. 1(b).

In both cases, the anisotropy factor AF increases with strain rate [small epsi, Greek, dot above], similarly to the experimental birefringence measurements. The main difference between the two simulations is that, for the monodisperse simulations with D_r = 6.73 s⁻¹, AF curve is shifted to higher extension rates compared to the AF obtained using D̄_r = 1.44 s⁻¹.

In the polydisperse simulations, the different averages A₀ to A₃ for AF are calculated using N = 10 bins based on the distribution of D_r (Fig. 1(b)), as described in Section 2.6. Importantly, as the weight function shifts from A₀ to A₃, giving higher weight to longer rods, the average contour length increases. For all polydisperse simulations, the average AF increases with strain rate and eventually reaches a plateau at high strain rates. Scaling the numerical AF data with α_s = 26 yields excellent agreement of A₃ with the experimental data across the entire strain rate range. The other averages also exhibit a qualitatively similar increase in AF. However, a progressive shift of AF toward higher extension rates is observed as the weight function shifts to lower powers of the rod length, consistent with the expected increase in the effective value of D_r. Overall, the birefringence experiments under steady flow conditions show good agreement with both the monodisperse simulation, based on the volume-averaged contour length (D̄_r = 1.44 s⁻¹), and the polydisperse simulation with A₃.

3.2 Rod dynamics under LAOE

3.2.1 Time-dependent birefringence. We next investigate the time-dependent response of the CNC systems to LAOE over a broad range of set strain rate amplitudes [small epsi, Greek, dot above]_0,set and oscillation periods T.

For time-dependent measurements, we define Pe₀ = [small epsi, Greek, dot above]₀/D̄_r using the extension rate amplitude [small epsi, Greek, dot above]₀ determined from the time-dependent strain rate in the PIV measurements. Moreover, we define the Deborah number De = f/D̄_r = (D̄_rT)⁻¹, which relates the diffusion timescale to the pulsation period in the time-dependent measurements.

Fig. 7 presents representative results for the temporal average birefringence, 〈Δn〉(t), of the 0.1% sample, compared side-by-side with the true extensional rate in the OSCER, [small epsi, Greek, dot above](t), retrieved experimentally through our PIV analysis. Note that the deviation between the true extensional rate, [small epsi, Greek, dot above](t), and the sinusoidal set strain rate at the pump, [small epsi, Greek, dot above]_set(t), arises as the syringe pumps approach their operational lower flow rate limit.³⁹ We emphasize that all measurements remain below the polarization camera's maximum retardance detection limit (≡Δn ≈ 7 × 10⁻⁵).

Fig. 7 Time-dependent strain rate [small epsi, Greek, dot above](t) (top sub-panels) and average birefringence 〈Δn〉(t) (bottom sub-panels) in the experiments under pulsatile LAOE. Example data is shown for the for 0.1% sample over two normalized pulsation cycles t/T at various [small epsi, Greek, dot above]_0,set for T = 2 s and T = 20 s. The residual birefringence 〈Δn〉_res that does not decay to zero during the cycle at higher frequencies is indicated in the bottom panels for the measurements at T = 2 s.

At low strain rates and large period durations, e.g., [small epsi, Greek, dot above]_0,set = 1 s⁻¹ and T = 20 s in Fig. 7, the measured birefringence 〈Δn〉(t) follows the sinusoidal waveform of the imposed strain rate [small epsi, Greek, dot above](t) along the extension axis. This indicates that the CNC orientational response to the evolving extension rate is linear and that the CNC has sufficient time to reach a steady-state orientation. It is, however, important to note that, although the birefringence responds linearly to the imposed extension rate, we are probing the out-of-equilibrium behavior of the CNC. Moving beyond this linear regime by increasing the set strain rate amplitude, [small epsi, Greek, dot above]_0,set, and/or decreasing the period duration, T, leads to a nonlinear response of the birefringence, manifested in two distinct fingerprints: the saturation of the birefringence and the onset of residual birefringence, 〈Δn〉_res, which we examine in more detail below.

Saturation of the birefringence. Increasing the set strain rate amplitude [small epsi, Greek, dot above]_0,set while keeping the period duration fixed at T = 20 s results in a progressive deviation from the sinusoidal waveform of the birefringence response. Although the strain rate still exhibits a sine, we observe a saturation of the birefringence profile 〈Δn〉 (t) at large [small epsi, Greek, dot above](t) during the cycle, e.g., at [small epsi, Greek, dot above]_0,set = 10 s⁻¹ and T = 20 s. This saturation of the birefringence at relatively high [small epsi, Greek, dot above](t) indicates that the extent of CNC alignment is approaching its maximum, analogous to the behavior observed under steady flow conditions (see Fig. 5(b) and (c)). Hence, the birefringence does not increase linearly with [small epsi, Greek, dot above](t) once it exceeds a certain value.
Residual birefringence. Upon increasing the pulsation frequency and applying a low strain rate amplitude, e.g., [small epsi, Greek, dot above]_0,set = 1 s⁻¹ and T = 2 s in Fig. 7, the birefringence response follows again the imposed strain rate profile. However, in contrast to the response at low frequencies (e.g., T = 20 s in Fig. 7), the birefringence does not decay back to zero as [small epsi, Greek, dot above](t) → 0 but instead exhibits a residual birefringence, 〈Δn〉_res, during the cycle (see red arrows in the bottom panels of Fig. 7), indicating a nonlinear response of the birefringence to [small epsi, Greek, dot above](t).

Based on the development of residual birefringence at low frequency, two contrasting regimes can be distinguished, depending on the pulsation period and the relaxation timescale of the rods. For sufficiently long pulsation periods, T ≫ 1/D̄_r, hence De ≪ 1, the rod dynamics closely resemble those observed under steady-state conditions, as the rods have enough time to adopt an average orientation comparable to that in steady flow during both the ascending and descending parts of the imposed sinusoidal strain rate cycle. Importantly, in this regime, the rods have sufficient time to fully relax and reach an anisotropic distribution at [small epsi, Greek, dot above] < D̄_r, with the birefringence signal decaying to the noise floor. In contrast, for relatively short pulsation periods, T ∼ 1/D̄_r, the pulsation timescale becomes comparable to or faster than the intrinsic rod relaxation timescale. As a consequence, the rod dynamics are highly transient and strongly dependent on their orientational history. This manifests itself as a progressively increasing asymmetry in the birefringence response with a decreasing pulsation period, together with a finite residual alignment persisting even at [small epsi, Greek, dot above] = 0.

3.2.2 Simulations. Fig. 8 compares the experimental LAOE results of the 0.1% sample against the corresponding numerical simulations for monodisperse and polydisperse systems. Representative data is shown for various T and [small epsi, Greek, dot above]_0,set over one normalized pulsation cycle t/T.

Fig. 8 Comparison of the numerical simulations and experiments under LAOE. Normalized time-dependent average birefringence signal 〈Δn〉/ϕ for the experiment and anisotropy factor AF/α_s for the simulations under pulsatile LAOE. Exemplary experimental data are shown for the 0.1% CNC sample, and the simulations show data for a monodisperse system with D̄_r = 1.44 s⁻¹ and for the polydisperse systems. Data is shown for various T and [small epsi, Greek, dot above]_0,set over one normalized pulsation cycle t/T.

At low strain rate amplitudes and large pulsation periods, e.g., [small epsi, Greek, dot above]_0,set = 1 s⁻¹ and T = 50 s in Fig. 8, the numerical data AF/α_s follows a sinusoidal profile throughout the cycle, similar to the time-dependent average birefringence 〈Δn〉/ϕ observed in the experiment. In this case, the monodisperse simulations show excellent agreement with the experimental data. Similar to the steady flow conditions, A₃ provides the best agreement for the polydisperse simulations, while the other averages progressively underestimate the birefringence. This trend among the different averages is observed for most tested combinations of strain rate amplitudes and pulsation periods under LAOE.

Increasing the set strain rate amplitude [small epsi, Greek, dot above]_0,set while keeping the period fixed at T = 50 s (Fig. 8) leads to a progressive deviation from the sinusoidal waveform of AF/α_s in both monodisperse and polydisperse simulations, in qualitative agreement with the experimental birefringence response. At [small epsi, Greek, dot above]_0,set = 100 s⁻¹, the polydisperse simulation with A₃ captures the experimental data very well. In contrast, the monodisperse simulation exhibits a faster relaxation and rise of AF/α_s during the descending and ascending phases of the cycle compared to the experiment. A qualitatively similar evolution of AF/α_s during the sinusoidal strain rate modulation is observed for both monodisperse and polydisperse systems across all simulations with increasing strain rate amplitude.

At T = 2 s, we also observe a transition from a sinusoidal to a more saturated (i.e., flat-topped) waveform in simulations as the strain rate amplitude is increased. In addition, similar to the experimental birefringence, a distinct residual in the anisotropy factor AF_res/α_s appears in all simulations as the temporal strain rate approaches zero during the cycle. The magnitude of AF_res/α_s increases from A₀ to A₃ in the polydisperse simulations and also rises with increasing strain rate amplitude at a fixed LAOE period.

Overall, the temporal evolution of AF/α_s during the LAOE cycle qualitatively matches the experimental observations. Both the saturation at high strain rates and the residual birefringence when the flow temporarily ceases are well captured in the simulations. Among them, the polydisperse simulation with A₃ shows the best agreement with the experimental data, indicating the important influence of the longer rods on the overall birefringence signal.

3.2.3 Residual birefringence during LAOE. As shown in Fig. 7 and 8, a residual birefringence 〈Δn〉_res (or AF_res/α_s in the simulations) emerges during LAOE as the pulsation frequency is increased. To generalize the results of 〈Δn〉_res as a function of Pe₀ for different CNC concentrations and LAOE periods, we employ a Pipkin-space representation. To further explore the critical conditions for the emergence of residual birefringence during the LAOE cycle, we identify the values of Pe₀ and De at which a non-zero residual birefringence, exceeding the noise threshold, first appears. Fig. 9 shows a Pipkin space representation of the experimental results across various Pe₀–De combinations. Colored symbols indicate measurements where a measurable residual birefringence is observed, while gray symbols correspond to measurements where the residual birefringence remains within the lower detection limit of the polarization camera. At low pulsation frequencies, e.g., De < 0.02, a residual birefringence emerges only above Pe₀ ≳ 20. As the LAOE frequency, and thus De, increases, the onset of residual birefringence shifts toward lower critical Péclet numbers. For De ≈ 0.05, the critical Pe₀ ≈ 1. At even higher pulsation frequencies (De ≥ 0.7), all experiments exhibit a measurable residual birefringence, regardless of how small the strain rate amplitude is within the explored Pe₀–De space. In the SI (see Fig. S1), we provide the plots of 〈Δn〉_res as a function of Pe₀ for different CNC concentrations.

Fig. 9 Pipkin space representing all experiments for the different CNC systems at different Pe₀ and De. Colored symbols indicate measurements where a residual birefringence is measurable, i.e., larger than the polarization camera's lower detection limit. Gray symbols show measurements where the residual birefringence is within the polarization camera's lower detection limit.

3.2.4 Lissajous curves. Another common approach to visualize the material response to a sinusoidally varying deformation is through Lissajous plots. Fig. 10 presents Lissajous curves of the normalized average birefringence, 〈Δn〉′ = 〈Δn〉(t)/〈Δn〉_max, plotted against the normalized extension rate, [small epsi, Greek, dot above]′ = [small epsi, Greek, dot above](t)/[small epsi, Greek, dot above]_max, for representative values of T and [small epsi, Greek, dot above]_0,set for the 0.1% CNC sample. Experimental data in Fig. 10 are compared with the normalized anisotropy factor, AF′ = AF(t)/AF_max, obtained from polydisperse simulations with A₃ averaging, which showed the best agreement with experiment under both steady and time-dependent flow conditions (see Fig. 6 and 8). The data are displayed on log–log scales to emphasize the behavior as [small epsi, Greek, dot above] → 0. Overall, we observe good qualitative agreement between simulations and experiments. In particular, the simulations capture both the hysteretic behavior and the residual birefringence. Since inter-rod interactions are neglected in the simulations, the presence of hysteresis and residual birefringence suggests that these features do not necessarily originate from inter-rod interactions, but can already arise from the dynamics of single rods under time-dependent flow. The deviation between simulations and experiments can arise from multiple factors, including the presence of inter-rod interactions (neglected in the simulations). Since we do not have experimental evidence for significant inter-rod interactions (discussed in Section 3.1.2), we conclude that the most probable cause of this deviation is the “non-ideality” of our CNC system, as the CNCs are not perfect cylindrical rods and exhibit surface heterogeneity.

Fig. 10 Lissajous curves of the normalized average birefringence 〈Δn〉′ from the experiments and the normalized anisotropy factor AF′ from the polydisperse simulation based on A₃ are plotted as a function of the normalized extension rate [small epsi, Greek, dot above]′. Experimental data are shown for the 0.1% sample. Data is represented in log–log scale at representative values of T and [small epsi, Greek, dot above]_0,set.

Although the Lissajous curves encode the same information as the time-resolved data shown in Fig. 7 and 8, they emphasise specific nonlinear features of the response, such as the asymmetry of the birefringence during the pulsation period. In particular, the Lissajous plots in Fig. 10 highlight the increasing hysteretic behaviour with pulsation frequency and, hence, with the Deborah number De. While this behaviour is visible in Fig. 10 for representative values of [small epsi, Greek, dot above]_0,set and T, its strong dependence on De becomes evident when exploring a broader parameter space in [small epsi, Greek, dot above]_0,set and T, as shown in Fig. S2 of the SI. The strong De-dependence of the hysteretic loop is further highlighted by the fact that hysteretic behavior appears only above a critical value of De (Fig. S2 of the SI). We note that both the residual birefringence 〈Δn〉_res and the hysteretic loop in the birefringence are signatures of the orientational history dependence of the CNCs, and are therefore both strongly dependent on De. Additionally, plots in Fig. 8, and Fig. S2 of the SI, display distinctive nonlinear loop shapes, reminiscent of those reported in large-amplitude oscillatory uniaxial extension (LAOE) of complex fluids. For example, Dessi et al.⁶⁹ reported nonlinear stress–strain behavior in elastomers attributed to extension-induced thickening. Similar nonlinear loops have been observed in soft polymer networks and polymer melts.^70,71 More recent studies on dilute polymer solutions under LAOE have revealed a nonlinear dependence of stress on strain rate at large deformation amplitudes.³⁹ In contrast to these systems, where nonlinearity originates from the macroscopic bulk response (e.g., the stress), the curved Lissajous figures observed in the present work (see Fig. 10) primarily arise from the CNCs responding nonlinearly to the imposed extensional rate, as a result of both saturation of the alignment at high Pe₀ and the fast pulsation time scale dictated by De.

4 Conclusions

This study investigates the dynamics of model colloidal rod-like systems under large amplitude oscillatory extension (LAOE). A microfluidic OSCER device was used to generate a homogeneous planar extensional flow, through which the extension rate was modulated sinusoidally using precision programmable syringe pumps. Our experimental rod-like colloidal system is based on dilute suspensions of well-characterized cellulose nanocrystals (CNC). The time-dependent strain rate within the OSCER was characterized using micro-particle image velocimetry. The orientation of the colloidal rods in response to the time-dependent strain rate was determined using flow-induced birefringence imaging, performed with a high-speed polarization camera.

Under steady flow conditions, the birefringence increases with strain rate until saturating above Pe ≳ 10 as the CNC rods become highly aligned. This behavior is consistent with previous studies on steady planar extension of CNC suspensions and is also well captured by both our monodisperse and polydisperse numerical simulations, with the rotational diffusion coefficient matched to the experimental system.

At sufficiently low extension-rate (or Péclet number) amplitude ([small epsi, Greek, dot above]_0,set) and low pulsation frequency, the birefringence response is linear and therefore remains sinusoidal. However, at sufficiently high strain-rate amplitude and/or high pulsation frequency, the birefringence response becomes nonlinear. This nonlinearity is evidenced by two main fingerprints: (i) at high [small epsi, Greek, dot above]_0,set, intra-cycle saturation of the birefringence occurs at large instantaneous extension rates. This effect is associated with the rods approaching their maximum alignment early in the cycle. (ii) In addition, nonlinearity is manifested by a residual birefringence that does not decay to zero during the cycle at sufficiently high pulsation frequencies, and hence at relatively high Deborah numbers. These high-De effects arise because the rods do not have sufficient time to equilibrate under the rapidly varying strain rate.

The simulations show that the polydisperse case A₃, which is most strongly weighted toward the longest rods in the CNC length distribution, provides the best overall qualitative agreement with the experimental data, capturing both the linear and nonlinear fingerprints of the birefringence response.

In these experiments, the birefringence response could only be related to the instantaneous extensional rate, and no direct information on the stress could be obtained. In future experiments, it would be highly interesting to perform concomitant measurements of the evolving stress and birefringence during LAOE, in order to assess to what extent analytical frameworks developed for LAOS are applicable to LAOE (e.g., the stress orientation relationship by Rogers et al.²⁵).

We believe that these original results should provide a firm foundation for understanding, predicting, and optimizing colloidal particle alignment in processing applications involving complex transient extensional and mixed flows.

Author contributions

SMR: conceptualization, methodology, investigation, visualization, writing – original draft, writing – review & editing. VC: writing – original draft, writing – review & editing. AQS: writing – review & editing, supervision, funding acquisition. GN: conceptualization, software, writing – original draft, writing – review & editing. SJH: conceptualization, writing – review & editing, supervision, funding acquisition.

Conflicts of interest

There are no conflicts to declare

Data availability

Supplementary information: figures providing (1) additional detail on the reported “residual birefringence” and (2) experimental Lissajous plots over an expanded range of imposed flow conditions than shown in the main text. See DOI: https://doi.org/10.1039/d5sm01122a.

The data supporting this study's findings are available within the article.

Acknowledgements

SMR, VC, AQS, and SJH gratefully acknowledge the support of Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan, along with additional financial support from the Japanese Society for the Promotion of Science (JSPS, Grant No. 24K07332, 24K17736, and 24K00810). VC also acknowledges support from the “la Caixa” Foundation (ID 100010434), fellowship code LCF/BQ/PI25/12100026. The authors thank Dr Yuto Yokoyama from OIST for fruitful discussions.

Notes and references

1. B. D. Leahy, D. L. Koch and I. Cohen, J. Rheol., 2017, 61, 979–996.
2. D. Kumar, A. Shenoy, S. Li and C. M. Schroeder, Phys. Rev. Fluids, 2019, 4, 114203.
3. C. Lang, J. Kohlbrecher, L. Porcar, A. Radulescu, K. Sellinghoff, J. K. G. Dhont and M. P. Lettinga, Macromolecules, 2019, 52, 9604–9612.
4. V. Calabrese, S. Varchanis, S. J. Haward and A. Q. Shen, Macromolecules, 2022, 55, 5610–5620.
5. V. Calabrese, A. Q. Shen and S. J. Haward, Biomicrofluidics, 2023, 17, 021301.
6. V. Calabrese, T. Porto Santos, C. G. Lopez, M. P. Lettinga, S. J. Haward and A. Q. Shen, Phys. Rev. Res., 2024, 6, L012042.
7. K. Osawa, V. K. Gowda, T. Rosén, S. V. Roth, L. D. Söderberg, J. Shiomi and F. Lundell, Flow, 2025, 5, E12.
8. T. Rosén, N. Mittal, S. V. Roth, P. Zhang, F. Lundell and L. D. Söderberg, Soft Matter, 2020, 16, 5439–5449.
9. M. Doi and S. F. Edwards, J. Chem. Soc., Faraday Trans. 2, 1978, 74, 560.
10. R. Ghanbari, A. Terry, S. Wojno, M. Bek, K. Sekar, A. K. Sonker, K. Nygård, V. Ghai, S. Bianco and M. Liebi, et al., Adv. Sci., 2025, 12, 2410920.
11. M. Detert, T. Porto Santos, A. Q. Shen and V. Calabrese, Biomacromolecules, 2023, 24, 3304–3312.
12. P. T. Corona, N. Ruocco, K. M. Weigandt, L. G. Leal and M. E. Helgeson, Sci. Rep., 2018, 8, 15559.
13. V. Lutz-Bueno, J. Zhao, R. Mezzenga, T. Pfohl, P. Fischer and M. Liebi, Lab Chip, 2016, 16, 4028–4035.
14. A. Rodriguez-Palomo, V. Lutz-Bueno, M. Guizar-Sicairos, R. Kadar, M. Andersson and M. Liebi, Addit. Manuf., 2021, 47, 102289.
15. P. T. Corona, B. Berke, M. Guizar-Sicairos, L. G. Leal, M. Liebi and M. E. Helgeson, Phys. Rev. Mater., 2022, 6, 045603.
16. N. Mittal, F. Ansari, K. Gowda, V. C. Brouzet, P. Chen, P. T. Larsson, S. V. Roth, F. Lundell, L. Wagberg, N. A. Kotov and L. D. Söderberg, ACS Nano, 2018, 12, 6378–6388.
17. Y. Liu, K. Zografos, J. Fidalgo, C. Duchêne, C. Quintard, T. Darnige, V. Filipe, S. Huille, O. Du Roure, M. S. N. Oliveira and A. Lindner, Soft Matter, 2020, 16, 9844–9856.
18. J. D. Ferry, W. M. Sawyer and J. N. Ashworth, J. Polym. Sci., 1947, 2, 593–611.
19. K. Hyun, S. H. Kim, K. H. Ahn and S. J. Lee, J. Non-Newtonian Fluid Mech., 2002, 107, 51–65.
20. M. Wilhelm, Macromol. Mater. Eng., 2002, 287, 83–105.
21. A. S. Khair, J. Fluid Mech., 2016, 791, R5.
22. M. De Corato and G. Natale, Macromolecules, 2019, 52, 4907–4915.
23. M. Das, L. Chambon, Z. Varga, M. Vamvakaki, J. W. Swan and G. Petekidis, Soft Matter, 2021, 17, 1232–1245.
24. N. H. Cho, J. Shi, R. P. Murphy, J. K. Riley, S. A. Rogers and J. J. Richards, Soft Matter, 2023, 19, 9379–9388.
25. S. Rogers, J. Kohlbrecher and M. Lettinga, Soft Matter, 2012, 8, 7831–7839.
26. B. Lonetti, J. Kohlbrecher, L. Willner, J. Dhont and M. Lettinga, J. Phys.: Condens. Matter, 2008, 20, 404207.
27. G.-R. Huang, L. Porcar, R. P. Murphy, Y. Shinohara, Y. Wang, J.-M. Carrillo, B. G. Sumpter, C.-H. Tung, L. Ding and C. Do, et al., Appl. Crystallogr., 2025, 58, 637–658.
28. J. Vermant, H. Yang and G. G. Fuller, AIChE J., 2001, 47, 790–798.
29. M. De Corato and G. Natale, Macromolecules, 2019, 52, 4907–4915.
30. N. K. Reddy, G. Natale, R. K. Prud’homme and J. Vermant, Langmuir, 2018, 34, 7844–7851.
31. T. P. Santos, V. Calabrese, M. W. Boehm, S. K. Baier and A. Q. Shen, J. Colloid Interface Sci., 2023, 638, 487–497.
32. V. Calabrese, S. J. Haward and A. Q. Shen, Macromolecules, 2021, 54, 4176–4185.
33. N. M. Mehdipour, N. Reddy, R. J. Shor and G. Natale, Phys. Fluids, 2022, 34, 083317.
34. D. P. Pope and A. Keller, Colloid Polym. Sci., 1978, 256, 751–756.
35. P. G. De Gennes, J. Chem. Phys., 1974, 60, 5030–5042.
36. D. E. Smith, H. P. Babcock and S. Chu, Science, 1999, 283, 1724–1727.
37. Y. Zhou and C. M. Schroeder, Phys. Rev. Fluids, 2016, 1, 053301.
38. Y. Zhou and C. M. Schroeder, Macromolecules, 2016, 49, 8018–8030.
39. S. M. Recktenwald, T. P. John, A. Q. Shen, R. J. Poole, C. P. Fonte and S. J. Haward, J. Non-Newtonian Fluid Mech., 2025, 342, 105421.
40. D. Kumar, Orientational dynamics of anisotropic colloidal particles in a planar extensional flow, arXiv, 2023, preprint, arXiv:2307.00688 DOI:10.48550/arxiv.2307.00688.
41. R. Wagner, A. Raman and R. Moon, Cellulose, 2010, 7, 27.
42. P. Bertsch, S. Isabettini and P. Fischer, Biomacromolecules, 2017, 18, 4060–4066.
43. P. Bertsch, A. Sánchez-Ferrer, M. Bagnani, S. Isabettini, J. Kohlbrecher, R. Mezzenga and P. Fischer, Langmuir, 2019, 35, 4117–4124.
44. M. Doi and S. F. Edwards, The Theory of Polymer Dynamics, Oxford University Press, 1988, vol. 73, pp. 289–323.
45. S. L. Anna and G. H. McKinley, J. Rheol., 2001, 45, 115–138.
46. V. Calabrese, A. Q. Shen and S. J. Haward, Macromolecules, 2024, 57, 9668–9676.
47. C. Lang, J. Hendricks, Z. Zhang, N. K. Reddy, J. P. Rothstein, M. P. Lettinga, J. Vermant and C. Clasen, Soft Matter, 2019, 15, 833–841.
48. M. A. Alves, AIP Conf. Proc., AIP, 2008, vol. 1027, pp. 240–242.
49. S. J. Haward, G. H. McKinley and A. Q. Shen, Sci. Rep., 2016, 6, 33029.
50. S. J. Haward, M. S. N. Oliveira, M. A. Alves and G. H. McKinley, Phys. Rev. Lett., 2012, 109, 128301.
51. F. J. Galindo-Rosales, M. S. N. Oliveira and M. A. Alves, RSC Adv., 2014, 4, 7799.
52. S. J. Haward, Biomicrofluidics, 2016, 10, 043401.
53. S. J. Haward, S. Varchanis, G. H. McKinley, M. A. Alves and A. Q. Shen, J. Rheol., 2023, 67, 1011–1030.
54. S. M. Recktenwald, C. Wagner and T. John, Lab Chip, 2021, 21, 2605–2613.
55. S. T. Wereley and C. D. Meinhart, Annu. Rev. Fluid Mech., 2010, 42, 557–576.
56. S. T. Wereley and C. D. Meinhart, Microscale Diagnostic Tech., Springer-Verlag, 2019, pp. 51–112.
57. C. D. Meinhart, S. T. Wereley and M. H. B. Gray, Meas. Sci. Technol., 2000, 11, 809–814.
58. T. Onuma and Y. Otani, Opt. Commun., 2014, 315, 69–73.
59. G. B. Jeffery, Proc. R. Soc. London, Ser. A, 1922, 102, 161–179.
60. L. G. Leal and E. J. Hinch, J. Fluid Mech., 1972, 55, 745–765.
61. G. G. Fuller, Optical rheometry of complex fluids, Oxford University Press, 1995.
62. J. Férec, M. Heniche, M.-C. Heuzey, G. Ausias and P. J. Carreau, J. Non-Newtonian Fluid Mech., 2008, 155, 20–29.
63. G. Natale, G. Ausias, J. Férec, M.-C. Heuzey and P. J. Carreau, J. Rheol., 2016, 60, 1069–1083.
64. S. S. Rogers, P. Venema, L. M. Sagis, E. Van Der Linden and A. M. Donald, Macromolecules, 2005, 38, 2948–2958.
65. P. N. Dunlap and L. G. Leal, J. Non-Newtonian Fluid Mech., 1987, 23, 5–48.
66. D. Kiriya, R. Kawano, H. Onoe and S. Takeuchi, Angew. Chem., Int. Ed., 2012, 51, 7942–7947.
67. M. Trebbin, D. Steinhauser, J. Perlich, A. Buffet, S. V. Roth, W. Zimmermann, J. Thiele and S. Förster, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 6706–6711.
68. K. Uetani, H. Koga and M. Nogi, ACS Macro Lett., 2019, 8, 250–254.
69. C. Dessi, D. Vlassopoulos, A. J. Giacomin and C. Saengow, Rheol. Acta, 2017, 56, 955–970.
70. H. K. Rasmussen, P. Laillé and K. Yu, Rheol. Acta, 2008, 47, 97–103.
71. A. G. Bejenariu, H. K. Rasmussen, A. L. Skov, O. Hassager and S. M. Frankaer, Rheol. Acta, 2010, 49, 807–814.

---