Structures and topological defects in pressure-driven lyotropic chromonic liquid crystals

We inject an aqueous solution of 13 wt$\%$ DSCG into a rectangular microfluidic channel of length $l=$50 mm, width $w=$15 mm, and thickness $b=6.5\pm 1$ $\mathrm{\mu }$m. At this concentration, DSCG is in the nematic phase at room temperature $T=22.5\pm 0.5$ $°$C (59, 60). The optical birefringence of 13 wt$\%$ DSCG at rest, $\mathrm{\Delta }n=n_{\mathrm{e}}-n_{\mathrm{o}}$, is $-0.015$ at a wavelength ${\lambda }_w$ = 633 nm, where $n_{\mathrm{e}}$ and $n_{\mathrm{o}}$ are the extraordinary and ordinary refractive indices (61). The corresponding maximum retardance, ${\mathrm{\Gamma }}_{max}=\frac{2\pi }{{\lambda }_w}\mathrm{\Delta }nb$, is $0.98\pm 0.15$ rad. The liquid crystal is initially planar aligned along the direction of the flow (Materials and Methods). The flow is induced by injecting DSCG solutions at a volumetric flow rate $q$, ranging from $q$ = 0.07–25 $\mathrm{\mu }$L/min, controlled by a syringe pump (Harvard PHD 2000). The corresponding Ericksen numbers Er = $\frac{{\eta }_2\overline{\dot{\gamma }}{b}^2}{K_2}$ vary from 306 to 109,267, where Er characterizes the relative importance of the viscous forces to the elastic forces. Here, ${\eta }_2$ is the twist viscosity, $K_2$ is the twist Frank elastic constant, and $\overline{\dot{\gamma }}=q/\left(w{b}^2\right)$ is the average shear rate (5, 59).

Distinct structures emerge in the material upon the onset of the pressure-driven flow, as shown in the snapshots in Fig. 1A, which are imaged in polarizing optical microscopy with a static full-wave-plate optical compensator (560 nm) with the slow axis oriented at $45°$ to the crossed polarizers and in the direction parallel to the flow. Orange colors indicate that the director is parallel to the flow direction ($x$ direction), and blue colors indicate that the director is perpendicular to the flow direction (61). At low flow rate, DSCG is preferentially aligned perpendicular to the flow, adopting a log-rolling state with in-plane orientation angle $\phi =90°$, even though before the onset of flow the director is parallel to the flow direction $\phi =0°$ (Fig. 1B). This realignment of the director is a consequence of the tumbling character of DSCG and the significant anisotropy in the splay Frank elastic constant $K_1$ and twist Frank elastic constant $K_2$, where $K_1/K_2\approx 10$; the director reorients by a twist deformation toward the $y$ axis, instead of deforming in the shear plane, which would involve splay deformations (44, 62). With an increase in flow rate, domains appear where the director increasingly aligns in the direction of the flow. The characteristic size of these domains becomes systematically smaller with increasing flow rate.

To quantify these structures, we use PSIM to obtain a map of the effective optical retardance, as shown in Fig. 1C where the colors represent the value of the optical retardance averaged over the thickness of the channel ($z$ direction). We determine the characteristic size of the domains of varying retardance by calculating the normalized two-dimensional (2D) spatial autocovariance in the $x$ and $y$ directions, $C_x$ and $C_y$, as shown in Fig. 2A (see SI Appendix for details) (63). The observed decrease in domain size with increasing flow rate is reflected in a decay of the autocovariance at increasingly smaller $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$, which denote the shift in the $x$ and $y$ directions, respectively. For the two lowest and the highest flow rates probed, $C_x$ and $C_y$ exhibit a two-step decay suggesting a coexistence of structures of two characteristic sizes. We use a double compressed exponential fit to access the characteristic length scales $L_{x_1}$, $L_{x_2}$, $L_{y_1}$, and $L_{y_2}$, characterizing the average sizes of structures along the $x$ and $y$ directions. For the intermediate range of flow rates, we fit $C_x$ and $C_y$ with a single compressed exponential function, which yields $L_x$ and $L_y$. Details on the fits and fit parameters are provided in SI Appendix, text and Table S1. Remarkably, $L_x$ and $L_y$ decrease monotonically with the average shear rate for $4\hspace{0.17em}{\mathrm{s}}^{-1}<\overline{\dot{\gamma }}<500\hspace{0.17em}{\mathrm{s}}^{-1}$, as displayed in Fig. 2B. The average characteristic size, defined as $L=\sqrt{L_x{L}_y}$, exhibits a power-law dependence with the shear rate, $L\propto {\overline{\dot{\gamma }}}^{-0.19}$. The domains are elongated in the flow direction, with an approximately constant aspect ratio $L_x/L_y=1.8\pm 0.3$, as shown in Fig. 2C. A decrease or an increase in shear rate leads to different characteristics of the domain sizes and the aspect ratio. We focus our discussion on the intermediate range of shear rates.

To identify the nature of the structures, we need to understand the observed changes in the optical retardance, which could be attributed to three possible effects: 1) changes of the out-of-plane orientation angle $\theta$, 2) changes of the order parameter at topological defects, or 3) twist deformations along the $z$ direction (44). To reveal the dominant effect leading to the observed retardance and the nature of the structures formed in DSCG solutions, we perform hydrodynamic simulations of tumbling nematic liquid crystals using a hybrid lattice Boltzmann method, which allows us to access the director field for Ericksen numbers Er = $\frac{{\eta }_2\overline{\dot{\gamma }}{b}^2}{K_2}$ varying from 579 to 8,214, a range which overlaps with the lower flow rates probed in the experiments. At Er < 2,480, the directors are predominantly aligned perpendicular to the flow direction (see SI Appendix, text and Fig. S1 for details). With increasing Er, the directors gradually reorient toward the flow direction. Disclination loops nucleate in the flow, shown as blue lines denoting isosurfaces of order parameter 0.35 in Fig. 3A for Er = 7,438, where the dark rods denote the director field in the plane of the loop. To link the simulations to the experiments, we calculate the effective optical retardance averaged over the channel thickness, as shown in Fig. 3B (SI Appendix). This reveals that the low retardance regions correspond to disclination loops. The majority of the disclination loops are topologically neutral pure-twist disclination loops, where the rotation vector is parallel to the normal direction of the loop (55), as shown in Fig. 3C.

We can rationalize the formation of pure-twist disclination loops by considering the elastic powers of splay, twist, and bend deformations, which can be expressed as $P_{\mathrm{s}\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{y}}={\int }_{\mathrm{\Lambda }}d\mathrm{\Lambda }{\left(\nabla \cdot \mathbf{n}\right)}^2$, $P_{\mathrm{t}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{t}}={\int }_{\mathrm{\Lambda }}d\mathrm{\Lambda }{\left(\mathbf{n}\cdot \nabla \times \mathbf{n}\right)}^2$, and $P_{\mathrm{b}\mathrm{e}\mathrm{n}\mathrm{d}}={\int }_{\mathrm{\Lambda }}d\mathrm{\Lambda }{\left(\mathbf{n}\times \nabla \times \mathbf{n}\right)}^2$, where $\mathrm{\Lambda }$ is a control volume. The simulations indeed show that $P_{\mathrm{t}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{t}}$ is larger than the powers associated with splay and bend modes (SI Appendix, Fig. S2). Therefore, even though the tumbling property of DSCG allows for the formation of wedge disclinations (44), twist-type defects are dominant in pressure-driven flow. To further quantify the prevalence of pure-twist disclination loops, we determine the local winding of the director field along the loop, which is characterized by the twist angle $\beta$ between the rotation vector $\mathrm{\Omega }$ and the local tangent vector $\mathbf{\text{t}}$ of a disclination loop (Fig. 3D, Inset). A value of $\beta =\pi /2$ denotes a local twist winding; a value of $\beta =0$ or $\pi$ denotes a local $\pm 1/2$ wedge winding (55). For a pure-twist disclination loop, the distribution of $\beta$ is a delta function with $\beta =\pi /2$ everywhere along the loop. A wedge-twist disclination loop, by contrast, is characterized by a $\beta$ that continuously varies from 0 to $\pi$ and back to 0 upon one full revolution around the loop (see SI Appendix, text and Fig. S3 for details). The distribution of $\beta$ measured for the disclination loops emerging in the flow of DSCG solutions exhibits a peak at $\beta =\pi /2$, as shown in Fig. 3D, which reveals that pure-twist disclination loops are indeed prevalent compared to wedge-twist disclination loops. This is further evidenced by the finding that the rotation vector $\mathrm{\Omega }$ is almost uniformly parallel to the loop normal $\mathbf{\text{N}}$ (SI Appendix, text and Fig. S4). The emergence of pure-twist disclination loops is a consequence of the smallness of the twist Frank elastic constant, which favors the local buildup of twist distortions in the $z$ direction, from which the loops nucleate (30, 44, 64).

The pure-twist disclination loops form at the boundary between two regions of irreconcilable director alignments: Close to the two channel walls, the directors are aligned in the shear plane, while in a region at the center of the channel the directors adopt a log-rolling state, as shown in Fig. 3E. Due to the high stored energy in the director gradient at the interface between these two frustrated regions, the elastic stress is released by forming topological defects. This complex director field within the channel gap results from the tumbling character of DSCG and the nonuniform shear rate across the gap; DSCG is in a log-rolling state below a critical shear rate $\dot{{\gamma }_c}$, but rotates toward the shear plane above $\dot{{\gamma }_c}$ (44, 54), as schematically shown in Fig. 3F.

To understand the transition to a different regime of flow structures at lower shear rates, we note that twist deformations can, a priori, lead to two types of topological defects: twist walls aligned parallel to the $xy$ plane, and pure-twist disclination loops (55, 64). The energy required to form a twist wall is $E_{\text{twist wall}}\approx 0.4L_{w}b\sqrt{-{\alpha }_2{K}_2\dot{\gamma }}$, where $L_w$ is the length of the twist wall, $K_2\approx 0.4$ pN is the twist Frank elastic constant, and ${\alpha }_2\approx -1.66$ Pa$\cdot$s is the Leslie viscosity coefficient for 13 wt$\%$ DSCG solutions (59) (see SI Appendix for details). The energy required to form a pure-twist disclination loop is $E_{\text{pure}-\text{twist}}=\frac{\pi }{4}{K}_2{L}_p\ln \left(\frac{L}{a}\right)$ (64), where $L_p$ and $L$ are the perimeter and the diameter of a pure-twist disclination loop, respectively, and $a$ is the diameter of the defect core in the nematic phase, which we estimate to be approximately $0.1\hspace{0.17em}\mathrm{\mu }$m. Pure-twist disclination loops rather than twist walls form when $E_{\text{pure}-\text{twist}}<E_{\text{twist wall}}$; this condition is reached for a critical shear rate ${\dot{\gamma }}^{*}\approx 0.8\hspace{0.17em}{\text{s}}^{-1}$. Considering the nonuniform shear rates across the thickness of the channel, twist walls become negligible when the shear rate in the center region [within a nematic coherence length $\approx a$ (30)] reaches $0.8\hspace{0.17em}{\text{s}}^{-1}$, which corresponds to a critical average shear rate ${\overline{\dot{\gamma }}}^{*}=\frac{{\dot{\gamma }}^{*}b}{12a}\approx \mathrm{4}\hspace{0.17em}{\mathrm{s}}^{-1}$. This value indeed denotes the onset of the intermediate shear rate regime in Fig. 2B.

Having established the emergence of pure-twist disclination loops in the range of intermediate average shear rates now allows us to rationalize the observed power-law dependence of $L\propto {\overline{\dot{\gamma }}}^{-0.19}$. To do so, we consider the nucleation forces and annihilation forces acting on the loop. The nonuniform twist deformation across the gap favors the nucleation; the viscous force imposed by the flow acts to annihilate the nucleated defects. The energy input to nucleate $N$ pure-twist disclination loops is $E_{\text{n}}=N{E}_{\text{pure}-\text{twist}}$. Simulations have determined that the number of nucleated twist defects at a given time induced by shear flow scales as $N\propto {\left(\text{Er}-{\text{Er}}_{\text{c}}\right)}^{0.5}$ (65), where ${\text{Er}}_{\text{c}}$ is the critical Ericksen number above which defects nucleate and $\text{Er}=\frac{{\eta }_2\dot{\gamma }{b}^2}{K_2}$ is the Ericksen number that governs the nucleation of twist-type defects. We here set ${\left(\text{Er}-{\text{Er}}_{\text{c}}\right)}^{0.5}\approx {\text{Er}}^{0.5}$, as we focus on the range where $\text{Er}\gg {\text{Er}}_{\text{c}}$. This then yields the nucleation energy per unit volume as $e_n=E_{\text{n}}/\mathrm{\Lambda }\propto \frac{1}{\mathrm{\Lambda }}{\left(\frac{{\eta }_2\dot{\gamma }{b}^2}{K_2}\right)}^{0.5}{K}_2{L}_p\ln \left(\frac{L}{a}\right)$, where $\mathrm{\Lambda }$ is a control volume. The annihilation of pure-twist disclination loops is driven by the viscous force and resisted by the elasticity, and can be expressed as the sum of the viscous forces and the elastic forces in a control volume $\mathrm{\Lambda }$: $f_a=\frac{1}{\mathrm{\Lambda }}\left(\int {\eta }_2\dot{\gamma }dS+K_2\right)$, where $\int dS\propto L^2$ is the area occupied by a twist loop of size $L$. For our range of intermediate average shear rates, $\text{Er}\gg 1$, indicating that viscous effects dominate, thus $f_a\propto \frac{{\eta }_2\dot{\gamma }{L}^2}{\mathrm{\Lambda }}$. Balancing the nucleation and annihilation forces, $\frac{d{e}_n}{dL}+f_a=0$, gives an expression for the average characteristic size of the disclination loops:This scaling argument indeed yields a power-law exponent for $\overline{\dot{\gamma }}$ fairly close to that observed in experiments, $L\propto {\overline{\dot{\gamma }}}^{-0.19}$. The characteristic loop size is thus governed by a balance between the nucleation force and the annihilation force acting on the loop.

We can likewise understand the aspect ratio of the pure-twist disclination loops, $L_x/L_y\approx 1.8\pm 0.3$, as being due to the asymmetric elastic deformation that results from the nonnegligible splay-bend anisotropy of DSCG. We consider the aspect ratio to be dominated by the elastic relaxation related to the deformation of the director field at the boundary of the disclination loop. This is a justified assumption given that the timescales related to the loop fluctuations induced by the viscous torque are much shorter than those characterizing the elastic deformation. The director field within the plane of a pure-twist disclination loop is described by $\mathbf{n}=\left(\cos \theta \cos \phi ,\cos \theta \sin \phi ,\sin \theta \right)$, where $\theta =0°$. The director field outside the loop is uniform and predominantly along the $x$ direction, so that $\phi$ is a small angle close to $0°$. Inside the loop, the director field is likewise uniform, but $\phi$ > 0. With these assumptions, and realizing that the deformation at the boundary of a pure-twist disclination loop along the $x$ direction involves bend deformations, while that along the $y$ direction involves splay deformations, the nematodynamic equation (30, 64) along the $x$ direction then reads as follows (see SI Appendix for details):and along the $y$ direction:where ${\gamma }_1$ is the rotational viscosity, $\frac{{\partial }^2\phi }{\partial x^2}\propto \frac{1}{{L_x}^2}$, $\frac{{\partial }^2\phi }{\partial y^2}\propto \frac{1}{{L_y}^2}$, and $\frac{{\partial }^2\phi }{\partial z^2}\propto \frac{1}{b^2}$. We interrogate the characteristic length scales in the $x$ and $y$ directions related to the spatial gradient of the director field within a certain time window. The time scales related to the elastic deformation in the $x$ and $y$ directions then scale as $\mathrm{\Delta }{t}_x\propto \frac{{\gamma }_1}{\left(\frac{K_3}{{L_x}^2}+\frac{K_2}{b^2}\right)}$ and $\mathrm{\Delta }{t}_y\propto \frac{{\gamma }_1}{\left(\frac{K_1}{{L_y}^2}+\frac{K_2}{b^2}\right)}$, respectively. With $K_1,K_3\gg K_2$ (5), this simplifies to $\mathrm{\Delta }{t}_x\propto \frac{{\gamma }_1{L_x}^2}{K_3}$ and $\mathrm{\Delta }{t}_y\propto \frac{{\gamma }_1{L_y}^2}{K_1}$. At steady state, $\mathrm{\Delta }{t}_x=\mathrm{\Delta }{t}_y$, which yields the following:The value of 1.9 is in good agreement with the experimentally determined value of $1.8\pm 0.3$.

Our single-shot imaging technique PSIM allows us to resolve the dynamics of the pure-twist disclination loops. We calculate the normalized spatiotemporal autocovariance, $C_t$, which contains the coupled information of two contributions: the fluctuations of the disclination loops characterized by a fluctuation time ${\tau }_1$ and the translation of the disclination loops imposed by the background flow characterized by a translation time ${\tau }_2$. To remove the contribution from the background flow, we need to place ourselves in the frame of reference of the disclination loop. In this Lagrangian framework, we move the region of interest by $\mathrm{\Delta }x=V_f\mathrm{\Delta }t$ at each time lag $\mathrm{\Delta }t$, where $V_f$ is the velocity of the frame of reference. Only if $V_f$ is equal to the center of mass velocity of the disclination loop, $V^{*}$, we access the fluctuations. In our pressure-driven flow, the flow velocity varies across the thickness of the channel. As the location of the pure-twist disclination loops within the channel is unknown, $V^{*}$ is unknown. We thus calculate $C_t$ for different frame of reference velocities, as shown in Fig. 4A for $V_f$ equal to the average velocity across the channel $\overline{V}$, and determine the characteristic time $\tau$ by fitting to a stretched exponential function (details on the fit and fit parameters are provided in SI Appendix, text and Table S2). The relation between the fluctuation time ${\tau }_1$, the translation time ${\tau }_2$, and the characteristic time $\tau$ can be expressed as follows (see SI Appendix for details):where $\frac{1}{{\tau }_2}=\frac{\left|V^{*}-V_f\right|}{L_x}$ with $L_x$ the characteristic length scale along the $x$ direction. This expression indeed well describes the dependence of ${\tau }^{-1}$ on $V_f/\overline{V}$, as shown in Fig. 4B for different flow rates. The frame of reference velocity $V_f$ at which ${\tau }^{-1}$ reaches a minimum denotes the center of mass velocity, $V^{*}$, which is between 1.1 and 1.4$\overline{V}$. This indicates that the disclination loops are located in the bulk flow rather than near the channel walls, in agreement with our simulations. Fitting ${\tau }^{-1}$ with Eq. 5 yields the fluctuation time ${\tau }_1$ and the characteristic length scale along the $x$ direction, $L_x$, as shown in Fig. 4C. The agreement between $L_x$ from this fit and $L_x$ from the normalized spatial autocovariance validates our approach. The fluctuation time ${\tau }_1$ scales as $1/\overline{\dot{\gamma }}$. To understand this dependence, we consider the nematodynamic equations for the director field $\mathbf{n}=\left(\cos \theta \cos \phi ,\cos \theta \sin \phi ,\sin \theta \right)$. Given the high Ericksen number regime of our experiments, we here neglect the elastic contributions (30, 64):where ${\alpha }_2$ and ${\alpha }_3$ are the Leslie viscosity coefficients, and ${\gamma }_1$ is the rotational viscosity. We account for the tumbling character of nematic DSCG solutions by considering small out-of-plane perturbations ${\theta }_1$ and in-plane perturbations ${\phi }_1$ for directors aligned perpendicular to the flow direction: $\theta ={\theta }_1$ and $\phi =\frac{\pi }{2}+{\phi }_1$. Linearizing Eqs. 6 and 7 in terms of these perturbations yields a characteristic fluctuation time of the tumbling nematics: ${\tau }_f\approx \frac{3}{4}\frac{{\gamma }_1}{\sqrt{-{\alpha }_2{\alpha }_3}}\frac{1}{\dot{\gamma }}$ (see SI Appendix for details). Using ${\alpha }_2\approx -1.66$ Pa$\cdot$s and ${\alpha }_3\approx 0.03$ Pa$\cdot$s (59), we observe a good agreement with ${\tau }_1$ from our experiments, as shown by the black line in Fig. 4C. This shows that the fluctuations of the pure-twist disclination loops are a direct reflection of the tumbling dynamics of the director.

We dissolve DSCG (TCI America; purity $>98.0\%$) at 13.0 wt% in deionized water. The sample is sealed in a glass tube and heated until it reaches the isotropic phase, indicated by a turbidity change from turbid to transparent. The sample is subsequently cooled to room temperature ($22.5\pm 0.5°$C), where it is in the nematic phase (59, 60). The rectangular microfluidic channel consists of two glass plates separated by $6.5\pm 1\hspace{0.17em}\mathrm{\mu }$m spacers. To ensure a well-controlled initial condition of the liquid crystal, we use a protocol of surface rubbing, which induces a planar alignment of DSCG along the direction of the flow, where both glass plates are rubbed along the cell length direction using diamond particles of diameter $\approx 50$ nm (68).

After injecting the sample in the microfluidic channel, we allow it to relax for 1 h, until it appears black when imaged through crossed polarizers, where one polarizer is placed parallel to the channel direction. The flow is controlled by a syringe pump (Harvard PHD 2000) set to a volumetric flow rate $q$ ranging between $q$ = 0.07–25 $\mathrm{\mu }$L/min. Once the flow has reached steady state, we image the sample at the center line of the channel at a frame rate of 506 frames per second in a $250\times 250\hspace{0.17em}\mathrm{\mu }{\text{m}}^2$ region far from the inlet (20–30 mm) to avoid artifacts due to the injection protocol. At the lowest flow rate (0.07 $\mathrm{\mu }$L/min), we start the measurement 40–50 min after the onset of flow. At intermediate flow rates (0.1–0.5 $\mathrm{\mu }$L/min), we start the measurement after 15–20 min. At higher flow rates (1–25 $\mathrm{\mu }$L/min), we start the measurement after 5–10 min. These times are well within the steady-state regime, as determined in additional experiments.

We adopt a hybrid lattice Boltzmann method to simulate the pressure-driven flow of nematic DSCG solutions. This method has been used in prior studies of passive and active lyotropic nematics (16, 37, 69–71). The nematic is described by a symmetric and traceless tensorial order parameter, defined as follows (64):where $S$ is the scalar order parameter, $\mathbf{n}$ is the unit vector representing the local nematic orientation, and $\mathbf{I}$ is an identity tensor. The governing Beris–Edwards equation of the nematic microstructure (9) reads as follows (72):where $\mathbf{u}$ is the velocity vector and $\mathrm{\Gamma }$ is related to the rotational viscosity of the nematic ${\gamma }_1$ via $\mathrm{\Gamma }=2S_0^2/{\gamma }_1$, with $S_0$ the equilibrium scalar order parameter. The generalized advection term $\mathbf{S}\left(\mathbf{W},\mathbf{Q}\right)$ is defined as follows:where $\mathbf{A}=\left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^T\right)/2$ is the strain rate tensor, $\mathbf{\Omega }=\left(\nabla \mathbf{u}-{\left(\nabla \mathbf{u}\right)}^T\right)/2$ is the vorticity, and $\xi$ is a flow-alignment parameter. Here, we choose $\xi <3S_0/\left(2+S_0\right)$ to enter a flow-tumbling regime (73). The molecular field $\mathbf{H}$ is a symmetric, traceless projection of the functional derivative of the free energy of the nematic. Its index form reads as follows:in which the free energy functional is $F={\int }_{V}fdV$. The density $f$ consists of a short-range Landau–de Gennes component and a long-range elastic component (74):where $A_0$ and $U$ are material constants, $Q_{ij,k}$ denotes ${\partial }_k{Q}_{ij}$, and $L_1$ to $L_4$ are related to the Frank elastic constants through the following (71):where $K_1$, $K_2$, $K_3$, and $K_{24}$ denote the splay, twist, bend, and saddle-splay Frank elastic constants, respectively. Eq. 9 is solved using a finite difference method.

The hydrodynamic flow is governed by a momentum equation (73, 75, 76):where $\rho$ is the density, $\eta$ is the isotropic viscosity, and $P_0=\rho T-f$ is the hydrostatic pressure with $T$ being the temperature. The additional stress accounting for the nematic anisotropy is defined as follows (73, 75, 76):Eq. 14 is solved simultaneously via a lattice Boltzmann method over a D3Q15 grid (75). The simulation is performed in a rectangular box $\left[X,Y,Z\right]=\left[\mathrm{150,150,50}\right]$ in simulation units with periodic boundary conditions in the $x$ and $y$ directions, and no-slip and planar anchoring condition in the $z$ direction. A body force $g\in \left[2.5\times 10^{-5}-5\times 10^{-5}\right]$ is applied to generate a pressure-driven flow. Additional details on this method can be found in ref. 76. Typical simulation parameters are $\mathrm{\Gamma }=0.1$, $\eta =0.33$, $A=0.01$, $U=3.5$ corresponding to $S_0\sim 0.62$, $\xi =0.6$ giving rise to ${\alpha }_3/{\alpha }_2=-0.08$, and $\left[L_1,L_2,L_3,L_4\right]=\left[0.01,\mathrm{0,0.03247},0.01333\right]$ leading to $K_1:K_2:K_3:K_{24}=1:0.33:3:1$.

The optical design of PSIM is based on the combination of off-axis shearing interferometry and circular polarization microscopy, as shown in Fig. 5. A supercontinuum laser serves as the illumination source with center wavelength at 633 nm set by a bandpass filter with a 10-nm pass band. The excitation is transmitted through a quarter waveplate to produce circular polarization light before impinging on the birefringent sample. The scattered light is then collected by a microscope and transmitted through a second quarter wave plate to transform the electric field’s polarization states back to linear polarization. The electric field component parallel and perpendicular to the slow axis of the quarter wave plate encode specimen birefringence information. To quantify the ratio of these two electric field components, the image is duplicated by a diffraction grating, each polarization component is selected by a linear polarizer at the Fourier plane, and they are recombined after a linear polarizer set at $\mathrm{45}°$ to form an interferogram on the CMOS camera.

In the interferogram, interference fringes will only appear in the region that has birefringent signals. The retrieval of the optical retardance is thus related to the extraction of the interference fringe’s amplitude. Similar to quantitative phase microscopy (77), we extract the amplitude of the retrieved +1 order $E$, along with the amplitude of the DC term $A$ from the interferograms with a digital holography algorithm. The retardance $\mathrm{\Gamma }$ is extracted as follows:This straightforward polarization parameter retrieval algorithm avoids the amplification of noise while quantitatively mapping the retardance from a single interferogram.

During image processing, we down-sample each frame by a factor of 10 ($\mathrm{1,200}\times \mathrm{1,200}$ pixels to $120\times 120$ pixels) to reduce the data size. This significantly increases the processing speed, but negligibly affects the information retrieved from the images as the down-sampled pixel size is still comparable to the diffraction limit (the diffraction limit is $1.54\hspace{0.17em}\mathrm{\mu }$m, the down-sampled pixel size is $2.08\times 2.08\hspace{0.17em}\mathrm{\mu }{\text{m}}^2$). The time interval between two consecutive frames is 1.97 ms.