Generalized Scalar Auxiliary Variable Exponential Integrator for A Modified Landau-de Gennes Theory for Smectic Liquid Crystals

Abstract: The Smectic-A (SmA) phase is modeled by a modified Landau-de Gennes (mLdG) model proposed by Xia et al. [Phys. Rev. Lett., 126 (2021), 177801], in which a tensor order parameter $\mathbf{Q}$ for the orientational order is coupled with a real scalar $u$ characterizing the positional order. In this paper, we propose and analyze a novel, highly efficient, and unconditionally energy-stable numerical scheme for this coupled system by combining the generalized scalar auxiliary variable-exponential integrator (GSAV-EI) approach with a relaxed correction strategy. In particular, we reformulate the exponential time differencing time discretization into an equivalent quasi-implicit backward Euler-type structure, a pivotal step that eliminates the restrictive CFL mesh-ratio conditions of the original GSAV-EI method and enables a rigorous fully discrete error analysis. Theoretically, we rigorously establish the unconditional energy stability with respect to a modified discrete energy and the uniform boundedness of the numerical solutions $\mathbf{Q}$, along with optimal error estimates in both time and space. Comprehensive numerical experiments are presented to demonstrate the accuracy, efficiency, and structural preservation of the algorithm, as well as its capability in capturing complex topological defect dynamics.

• The paper introduces a GSAV-EI method that achieves unconditional energy dissipation for the modified Landau‐de Gennes model of smectic-A liquid crystals.
• It combines exponential integrator techniques with a relaxation strategy to decouple the coupled dynamics of the Q-tensor and smectic order parameter.
• Numerical experiments validate first-order temporal and second-order spatial convergence, accurately capturing defect dynamics and layer periodicity in both 2D and 3D.

Generalized Scalar Auxiliary Variable Exponential Integrator for a Modified Landau-de Gennes Theory for Smectic Liquid Crystals

Introduction and Motivation

The paper addresses the numerical analysis and simulation of the smectic-A (SmA) liquid crystal phase using a modified Landau-de Gennes (mLdG) theory, which couples an orientational $Q$-tensor field with a positional scalar order parameter $u$. The SmA phase, characterized by both nematic order and one-dimensional positional layering, requires extensions of the classical Landau-de Gennes theory to faithfully capture its complex defect and phase transition structures. The authors focus on constructing structure-preserving, unconditionally energy-stable numerical schemes for the coupled mLdG system to compute dynamics and topological transitions in SmA liquid crystals accurately and efficiently.

A significant innovation is the integration of the generalized scalar auxiliary variable-exponential integrator (GSAV-EI) approach with a relaxation strategy (R-GSAV-EI). Unlike prior exponential time differencing (ETD) or SAV schemes, this formulation combines rigorous energy stability with removals of restrictive Courant–Friedrichs–Lewy (CFL) constraints and admits a fully discrete error analysis. The work further provides optimal error estimates, maximum bound principles (MBP), and demonstrates the method's efficacy through a series of high-dimensional, physically relevant numerical experiments.

Mathematical Model

The mLdG energy functional for SmA phases is formulated as

$$E(\mathbf{Q}, u) = \int_\Omega \left[ f_n(\mathbf{Q}, \nabla \mathbf{Q}) + f_s(u) + f_{int}(\mathbf{Q}, u) \right] dx$$

where $\mathbf{Q}$ is a symmetric, traceless tensor capturing nematic alignment, and $u$ is a scalar representing smectic layering. The free energy includes nematic elastic and bulk terms, a scalar polynomial for smectic layering, and higher-order couplings encoding the interaction of orientational and positional degrees of freedom. The resulting $L^2$-gradient flow yields a coupled system of high-order, strongly nonlinear PDEs for $(\mathbf{Q},u)$.

A key physical requirement is unconditional energy dissipation:

$\frac{dE}{dt} \le 0$

and the preservation of admissible set constraints (symmetry, tracelessness, and, crucially, a maximum bound on $\mathbf{Q}$).

R-GSAV-EI Schemes and Theoretical Properties

Reformulation and Scheme Construction

The authors introduce a scalar auxiliary variable $s$ tracking the essential nonlinear part of the discrete energy $u$0, and recast the coupled system using a modified energy $u$1 and an exponential stabilization factor $u$2. This admits a linear, decoupled, unconditionally energy-stable time-marching scheme, with key updates:

• Exponential integrator type updates for $u$3 using operator splitting with stabilization,
• A relaxation-corrected update for the auxiliary variable $u$4, tightly tying the SAV to the original energy,
• Careful spectral norm analysis, operator splitting, and quasi-implicit discretization for optimal time stepping without CFL restrictions.

Energy Stability and Maximum Bound Principle

A central rigorous result is Theorem 1: the scheme guarantees unconditional energy dissipation for any step size and parameter regime, as the discrete modified energy $u$5 strictly decreases at each time step (except possible relaxation truncation errors controlled by a relaxation parameter). The MBP is established for the $u$6-tensor, ensuring all iterates remain within physically admissible bounds and obviating spurious blow-up or ordering defects.

Error Analysis

A comprehensive error analysis is conducted, featuring:

• Fully discrete $u$7, $u$8, and $u$9 norm a priori estimates,
• Uniform-in-time bounds for all solution components,
• Lipschitz stability for all nonlinear operators,
• Discrete Gronwall arguments demonstrating first-order convergence in time and second-order convergence in space, matching the underlying scheme's formal order,
• Regularity bootstrapping via analytic semigroup arguments to show robust performance in both 2D and 3D regimes.

Numerical Experiments

Numerical results illuminate the practical behavior and performance of the R-GSAV-EI scheme on representative 2D and 3D SmA liquid crystal dynamics.

Smectic Layer Evolution and Defect Dynamics

The method robustly captures the complex time evolution of smectic layers, including phase front propagation, defect nucleation and recombination, grain boundary formation, and relaxation to stable, defect-rich target patterns, as shown in the evolution of the $$E(\mathbf{Q}, u) = \int_\Omega \left[ f_n(\mathbf{Q}, \nabla \mathbf{Q}) + f_s(u) + f_{int}(\mathbf{Q}, u) \right] dx$$0-tensor field Figure 1.

Figure 1: Evolution of the smectic-A liquid crystal system, showing Q-tensor magnitude snapshots that reveal smectic layer formation and defect dynamics.

The method dynamically reflects topological transitions and the complex interplay of orientational and positional order, even during severe rearrangement events such as collision and recombination of grain boundaries.

Energy Dissipation, Temporal Accuracy, and Robustness

A detailed comparison of supremum norm and energy evolution Figure 2 demonstrates exact alignment of energy dissipation profiles and supremum norms between large and fine time step runs, highlighting minimal temporal error and unconditional stability.

Figure 2: Supremum norm and total energy over time, quantifying strict energy dissipation and convergence under both relaxed and standard schemes.

The relaxation-enhanced scheme eliminates nonphysical energy drift present in the standard SAV-EI formulation and closely follows the physical energy, substantiating the theoretical results on accuracy and stability.

Layer Periodicity and Physical Parameter Control

Simulations probing variations in $$E(\mathbf{Q}, u) = \int_\Omega \left[ f_n(\mathbf{Q}, \nabla \mathbf{Q}) + f_s(u) + f_{int}(\mathbf{Q}, u) \right] dx$$1 (wavenumber) and $$E(\mathbf{Q}, u) = \int_\Omega \left[ f_n(\mathbf{Q}, \nabla \mathbf{Q}) + f_s(u) + f_{int}(\mathbf{Q}, u) \right] dx$$2 (coupling strength) Figure 3 confirm the method's ability to resolve parameter-driven changes in intrinsic layer periodicity and director field arrangement.

Figure 3: Steady-state profiles of the density field $$E(\mathbf{Q}, u) = \int_\Omega \left[ f_n(\mathbf{Q}, \nabla \mathbf{Q}) + f_s(u) + f_{int}(\mathbf{Q}, u) \right] dx$$3 for various $$E(\mathbf{Q}, u) = \int_\Omega \left[ f_n(\mathbf{Q}, \nabla \mathbf{Q}) + f_s(u) + f_{int}(\mathbf{Q}, u) \right] dx$$4 and $$E(\mathbf{Q}, u) = \int_\Omega \left[ f_n(\mathbf{Q}, \nabla \mathbf{Q}) + f_s(u) + f_{int}(\mathbf{Q}, u) \right] dx$$5 demonstrating physically accurate layer spacing and director orientation.

3D Topological Dynamics

Large-scale 3D computations Figure 4 showcase the scheme's scalability and stability, capturing boundary nucleation, defect line recombination, and relaxation to complex concentric 3D smectic patterns without numerical breakdown.

Figure 4: Three-dimensional dynamic evolution at increasing $$E(\mathbf{Q}, u) = \int_\Omega \left[ f_n(\mathbf{Q}, \nabla \mathbf{Q}) + f_s(u) + f_{int}(\mathbf{Q}, u) \right] dx$$6, visualizing nested smectic layering and director organization during topological transitions.

Adaptive Time Stepping and Stability Comparison

Adaptive time-stepping profiles Figure 5 and stability contrasts versus standard ETD reveal that only the SAV-based schemes combine robust phase transition resolution with high efficiency (no step size chattering or numerical blow-up when adaptive stepping is activated).

Figure 5: Adaptive time step profiles, illustrating intelligent step-size control by sESAV versus erratic behavior in non-stabilized ETD.

Implications and Future Directions

This work advances the state of computational methods for complex multi-order parameter phase field models. In particular, it sets a new standard for numerical fidelity in smectic liquid crystal simulation by ensuring unconditional energy dissipation, maximum bound principle enforcement, and provable temporal and spatial error bounds. The R-GSAV-EI framework is highly extensible; its analytic structure and operator splitting approach are compatible with a broad family of high-order nonlinear dissipative systems, not limited to liquid crystals.

Practically, this enables high-throughput, physically reliable simulation of domain-scale systems in both 2D and 3D, which is critical for understanding topological defect evolution, grain boundary kinetics, and external field-driven pattern formation in smectics. The theoretical approach facilitates further developments in adaptive-in-time, adaptive-in-space, and multi-scale numerical algorithms for continuum complex fluids, and the SAV/GSAV methodology can be extended to stochastic and multiphase models.

Conclusion

The R-GSAV-EI method presented in this work robustly resolves the nonlinear, coupled dynamics of the mLdG SmA liquid crystal model, delivering unconditional energy stability, optimal error estimates, and strict conformance to the physical admissibility constraints throughout all computational regimes. Through comprehensive theoretical analysis and extensive computational experiments, the approach is shown to be effective for large-scale, high-fidelity dynamical simulations of smectic layers, defects, and phase transitions. These results have concrete implications for the design of next-generation structure-preserving integrators in liquid crystal modeling and other nonlinear dissipative continuum systems.