A 100× Faster Beam Propagation Method for Nonlinear Optical Wave Propagation
Based on Discrete Exterior Calculus

Abstract

Efficient simulation of nonlinear light propagation in complex photonic structures remains a major challenge because these systems combine intricate transverse geometries with propagation over distances spanning many diffraction lengths. Existing numerical methods often require computationally intensive uniform discretizations or struggle to accurately represent complex material boundaries, limiting the practical simulation of multiscale nonlinear photonic devices. Here we introduce a computational framework for solving the nonlinear Schrödinger equation that accelerates simulations by more than two orders of magnitude (over 100×) compared with conventional approaches while maintaining spectral-level accuracy, thereby reducing the computational time required to solve large-scale and multiscale problems from days to hours. The method is based on discrete exterior calculus, enabling geometry-conforming discretization directly on unstructured meshes without the weak formulations required by conventional finite-element methods. In contrast to Fourier-based spectral solvers, it avoids global oversampling, eliminates Gibbs-type oscillations at material interfaces, naturally incorporates absorbing boundary conditions, and preserves the topological structure of the underlying differential operators. A key feature of the framework is that higher-order propagation operators are generated algorithmically from lower-order discrete operators, providing a systematic route to extending simulations beyond the standard nonlinear Schrödinger equation. Benchmarks on fundamental solitons, asymmetric beams, and optical vortices confirm spectral-level accuracy while demonstrating computational speedups exceeding 100×. By combining geometric flexibility, computational efficiency, and a systematic framework for higher-order wave propagation, our approach enables routine simulations of nonlinear dynamics in previously inaccessible multiscale photonic systems, including photonic crystal fibers, large-core multimode waveguides, and other complex integrated photonic platforms. In addition, the reported simulation speedup enables numerical optimization of nonlinear optical structures beyond what is currently feasible using simplified models or intuition alone.

I Introduction

Nonlinear dynamics has become a central framework for understanding complex phenomena across science and engineering. Within this broad domain, nonlinear wave dynamics plays a particularly important role, governing systems in which wave evolution cannot be described by linear superposition. Such waves arise in a wide range of physical contexts: in fluid dynamics they underlie turbulence and soliton formation; in biology they describe processes such as nerve pulse propagation, cardiac dynamics, and morphogenesis; in chemical systems they appear as reaction–diffusion waves; and in atmospheric physics they govern large-scale structures such as Rossby waves and shock fronts. In optics, nonlinear wave dynamics is responsible for phenomena including self-focusing, soliton propagation, supercontinuum generation, and frequency comb formation, while in quantum systems it governs the behavior of Bose–Einstein condensates and matter-wave solitons.

A particularly important class of nonlinear wave phenomena involves the evolution of a transverse field profile along a preferred propagation direction under a nonlinear evolution equation. Such systems exhibit rich dynamics arising from the interplay between dispersion or diffraction and nonlinearity, giving rise to phenomena including soliton formation, frequency comb generation, and supercontinuum broadening.

Among the models used to describe these systems, the nonlinear Schrödinger equation (NLSE) provides one of the most universal and widely applicable frameworks. It governs the evolution of slowly varying wave envelopes in nonlinear dispersive media and emerges naturally under envelope and multiple-scale approximations in the weakly nonlinear regime. In this sense, the NLSE serves as a canonical model from which a variety of other nonlinear wave equations can be systematically derived or approximated under appropriate asymptotic limits. For example, in regimes where wave dynamics is dominated by unidirectional propagation and weak nonlinearity with long-wavelength dispersion, reductions of the NLSE can lead to effective equations such as the Korteweg–de Vries (KdV) equation. More broadly, different scaling limits and perturbative expansions connect the NLSE to a hierarchy of nonlinear evolution equations, highlighting its role as a unifying framework for nonlinear wave phenomena. In nonlinear optics, the NLSE describes light propagation in nonlinear media, capturing key effects such as self-phase modulation and soliton dynamics. In quantum gases, it appears as the Gross–Pitaevskii equation, modeling the mean-field behavior of Bose–Einstein condensates. In both cases, the NLSE captures the essential interplay between longitudinal evolution and transverse structure that governs nonlinear wave propagation.

The central role of the NLSE across multiple physical systems makes its accurate and efficient numerical solution essential for modeling nonlinear wave dynamics. In this work, we focus on nonlinear optics, specifically on nonlinear wave propagation in optical fibers and waveguides, which underpins a wide range of modern technologies, including laser science, optical communications, ultrafast optics for probing chemical dynamics, and integrated photonics, just to mention a few examples. However, simulating optical wave propagation remains challenging due to the intricate interplay between nonlinearity, dispersion or diffraction, and the multiscale spatial and temporal features inherent to many practical settings. These challenges are particularly pronounced in structured optical media, where fine transverse geometric features coexist with long propagation distances, as well as in highly multimode fibers, where modal interference can generate fine spatial structures, and in strongly nonlinear regimes where effects such as beam filamentation may occur. Consequently, developing computational methods that are both efficient and capable of preserving essential physical and topological properties is of paramount importance.

To address these challenges, a variety of numerical approaches have been developed to study wave evolution under the NLSE in optical systems, including finite-difference methods (FDM), split-step methods (SSM) [1], and spectral methods (SM) [2, 3], as well as more recent formulations based on the finite element method (FEM) [4, 5]. Among these, FDM is straightforward to implement but lacks flexibility, as the use of nonuniform meshes often degrades accuracy. Spectral and split-step methods are widely used in optics due to their efficiency for homogeneous or weakly structured systems; however, they rely on uniform grids and periodic boundary conditions, making them less suitable for complex geometries and heterogeneous media. FEM, in contrast, enables flexible meshing and can naturally accommodate geometries with multiple characteristic length scales. However, it relies on weak formulations, which become increasingly cumbersome for generalized NLSE (GNLSE) models that include multiple nonlinearities and higher-order diffraction terms. Moreover, FEM performance depends sensitively on the choice of basis functions and may suffer from stability or accuracy issues if these are not carefully selected [6]. Time integration within FEM frameworks typically involves implicit schemes such as Crank–Nicolson, which can become computationally expensive in nonlinear regimes. Hybrid approaches that combine FEM for the linear step with exponential treatments of the nonlinear step alleviate some of these issues but remain computationally demanding.

In this respect, both SSM and SM have become the gold standard for simulating the GNLSE, largely thanks to the fast Fourier transform (FFT). To fully exploit the efficiency of the FFT, the transverse computational grid must be uniform, ideally comprising $2^{n}$ nodes in each direction or, as supported by modern FFT implementations such as FFTW, a number of nodes that factorizes into small prime numbers. While algorithms for nonuniform grids exist, they are significantly slower than uniform-grid implementations. This limitation becomes critical in multiscale problems, where different regions require different resolutions. For example, if a nonlinear wave experiences filamentation in localized regions, the grid may need to be refined only in those regions, and for an optical wave propagating in a photonic crystal fiber, which contains fine geometric structures (holes) across the transverse direction, it would be advantageous to have a fine mesh resolving these structures while keeping a coarser mesh in the core or uniform cladding. FFT-based methods are generally incapable of handling these scenarios efficiently. Another challenge is that FFT methods automatically impose periodic boundary conditions, which can produce incorrect results if part of the wave is scattered or deflected to the boundaries, as it would artificially re-enter from the opposite side.

Figure 1: Schematics of three representative optical waveguide platforms widely used in nonlinear optics: step-index optical fibers (left), photonic crystal fibers (center), and integrated photonic ridge waveguides (right). These examples illustrate the geometric diversity and heterogeneity of practical waveguide structures. In this work, we present a versatile computational framework based on discrete exterior calculus (DEC) that enables efficient simulation of nonlinear optical wave propagation across such geometries.

It is therefore highly beneficial to develop a computational scheme for solving the GNLSE that incorporates the following features: (1) flexible meshing to accurately handle multiscale geometries, (2) a direct discretization of the governing equation without the need to formulate a weak form, (3) seamless compatibility with absorbing boundary conditions to eliminate artifacts from scattered waves, (4) the ability to extend naturally to vectorial NLSEs, and (5) compatibility with adaptive step-size evolution schemes. Such a technique will enable fast and accurate simulation of problems involving multiple length scales, such as photonic crystal fibers or large fiber cores supporting many optical modes, both of which have recently attracted significant interest in nonlinear optics.

In this work, we develop a computational scheme based on discrete exterior calculus (DEC) that achieves all of the aforementioned objectives. By leveraging the inherent structure-preserving properties of DEC, our approach enables flexible meshing capable of resolving multiscale geometries, direct discretization of the governing nonlinear evolution equations without the need for a weak form, and seamless integration with absorbing boundary conditions to prevent spurious reflections. While in this work, we focus on nonlinear optical systems, the computational scheme introduced here can be applied equally well to the study of nonlinear matter waves in Bose–Einstein condensates.

II Results

In this section, we present the conceptual and mathematical framework of the discrete exterior calculus beam propagation method (DEC-BPM), while deferring detailed derivations to the Supplementary Information. We then demonstrate its performance through representative examples. Before proceeding to these details, it is instructive to first provide a high-level comparison between DEC-BPM, standard spectral BPMs, and FEM-based BPMs, as summarized in Table 1.

Table 1: Comparison of numerical methods for beam propagation and waveguide analysis, highlighting the structural and mathematical distinctions of the proposed DEC-BPM framework.

Feature	Spectral Method	Finite Element Method	Discrete Exterior Calculus
Mesh/Grid	Uniform Cartesian: Costly for multiscale features.	Unstructured Primal: Adaptable to complex geometries.	Primal-Dual Complex: Barycentric/Voronoi duals.
Formulation	Utilizing global basis functions, FFT-based.	Weak Form via shape functions.	Direct discretization of differential forms. Preserves exact metric-free topology ($\mathbf{d}^{2}=0$).
Propagation	Typically explicit; it may use Runge-Kutta propagation, though implicit schemes are permissible.	Typically implicit (e.g., Crank-Nicolson via a split-step approximation). Explicit schemes require mass matrix inversion or artificial mass lumping.	Naturally explicit: Geometrically constructs a diagonal Hodge star $\fourierbb{H}_{0}$, enabling explicit ODE integration (e.g., via RK45) without mass matrix inversion.
Laplacian ($\nabla^{2}$)	$-k^{2}$ (Spectral derivative)	Stiffness matrix via weak form: $\int\nabla u\cdot\nabla v\,d\Omega$	Combinatorial incidence & Hodge stars matrices: $-\fourierbb{H}_{0}^{-1}\fourierbb{D}_{0}^{\mathsf{T}}\fourierbb{H}_{1}\fourierbb{D}_{0}$
Higher-order ($\nabla^{2n}$)	$(-k^{2})^{n}$ (Spectral derivative)	Requires $\boldsymbol{H}^{n}$-conforming elements or recursive mixed formulation splitting.	Cascaded matrix multiplication: $(-\fourierbb{H}_{0}^{-1}\fourierbb{D}_{0}^{\mathsf{T}}\fourierbb{H}_{1}\fourierbb{D}_{0})^{n}$
Accuracy	Exhibits Gibbs phenomenon at sharp index steps. Less optimal for complex beam profiles.	Accurate when explicitly aligned with material interfaces, but computationally heavier for propagation.	Captures sharp contrasts accurately (step-index). Highly convergent for non-standard inputs (Elliptic Gaussian, vortex beams).

A key advantage of DEC-BPM is its natural compatibility with unstructured meshes, which enables accurate representation of complex geometries and sharp material boundaries. This flexibility also allows the computational domain to conform to the physical structure (e.g., circular or irregular cross-sections), avoiding the artificial rectangular truncation required by Cartesian grids in spectral methods. Such truncation leads to unnecessary domain padding, increased memory usage, and higher computational cost. Moreover, the DEC framework readily accommodates absorbing boundary conditions, which are essential for modeling unguided radiation and leakage in open systems. In contrast, spectral methods inherently assume periodic boundary conditions, making them ill-suited for such scenarios.

Compared to FEM-based approaches, DEC offers a significant conceptual and practical simplification: the discretization proceeds directly from the governing equations without requiring a weak formulation. This facilitates the treatment of generalized NLSE models with multiple nonlinear terms. In addition, higher-order diffraction operators can be constructed systematically through repeated application of discrete operators associated with lower-order terms, avoiding the stringent basis function requirements typically encountered in FEM. Taken together, these features enable the use of explicit time-stepping schemes, such as fourth-order Runge–Kutta (RK4), for propagation. The compatibility of such schemes with adaptive step sizing further enhances computational efficiency, allowing faster simulations while maintaining high accuracy.

Conceptual Framework:— As discussed above, we focus on modeling nonlinear optical wave propagation in guiding structures using the nonlinear Schrödinger equation (NLSE). In normalized form, this equation can be written as:

$$i\frac{\partial\psi}{\partial z}=-\sum_{n}p_{n}\nabla_{\perp}^{2n}\psi-q_{L}V(x,y)\,\psi-q_{NL}\mathcal{N}(|\psi|^{2})\,\psi,$$ (1)

where $\psi(x,y,z)$ is the normalized field envelope, $z$ is the dimensionless propagation coordinate, and $\nabla_{\perp}^{2}=\partial_{x}^{2}+\partial_{y}^{2}$ denotes the transverse Laplacian. The terms $\nabla_{\perp}^{2n}$, with $n$ an integer, account for higher-order diffraction effects (typically, only the leading term is retained for weakly diffracting beams). The linear potential $V(x,y)$ represents the transverse refractive index profile of the waveguide, while $\mathcal{N}(|\psi|^{2})$ describes the nonlinear response. For Kerr media, the nonlinearity takes the form $\mathcal{N}(|\psi|^{2})=\chi(x,y)\,|\psi|^{2}$, where $\chi(x,y)$ is the Kerr coefficient, which may vary across the transverse plane. Finally, the constants $p_{n}$, $q_{L}$, and $q_{NL}$ ensure consistency of units across all terms. In practice, Eq. 1 is often expressed in normalized units, which we adopt after developing the formalism. Although the derivation and normalization of Eq. 1 are standard, they are included for completeness in Supplementary Note 1.

The first step in using DEC with any differential equation is to recast the equation in the language of exterior calculus (EC) [7, 8, 9]. A key idea in this formulation is that different physical quantities naturally live on different geometric objects. For example, a scalar field (such as refractive index or intensity) is associated with points in space, while differences of a field occur along edges, and fluxes naturally pass through areas (faces). This is not merely a mathematical construction—it reflects how these quantities are defined physically. This organization is formalized using differential forms of different orders. In an $n$-dimensional space, a $k$-form represents a quantity that is naturally associated with $k$-dimensional objects: $0$-forms correspond to pointwise values (scalars), $1$-forms to quantities along edges, $2$-forms to quantities over surfaces, and so on up to $n$-forms, which are associated with volumes. The order of the form therefore reflects both the type of physical quantity and the dimensionality of the domain it interacts with. This viewpoint provides a unified way of describing different physical quantities based on how they interact with the geometry of the domain. The above discussion illustrates, at a conceptual level, why DEC in two dimensions employs two complementary meshes (see table 1. Each mesh represents different types of quantities that naturally reside on distinct geometric elements, such as edges or faces.

Next, the various operations—whether differential operators or linear and nonlinear mappings—must be expressed in terms of a set of standard exterior calculus (EC) operators that relate differential forms of different orders. These operators are defined in a way that separates the topology of the problem from its metric properties. One of these fundamental operators is the exterior derivative, $\mathbf{d}$, which unifies the notions of gradient (when acting on $0$-forms), curl (when acting on $1$-forms), and divergence (when acting on $2$-forms) into a single, coordinate-independent operation. Importantly, this operator is purely topological, meaning it depends only on how points are connected, rather than on distances or angles. Figure 2 illustrates this concept from both continuous and discrete perspectives. Panel (a) shows the scalar function $f(x,y)$ in the original coordinates $(x,y)$, while panel (b) depicts the same function in the deformed coordinate system $(X,Y)$. Although the scalar field itself is unchanged, its gradient and directional derivatives generally differ in the two coordinate systems. Nevertheless, the integral quantity $\oint df$ around any closed loop always vanishes. In numerical analysis, one is often interested in approximating quantities such as $\nabla f$. However, many discretization schemes do not guarantee that the identity $\nabla\times\nabla f=0$ is exactly preserved, even though it should hold independently of the chosen coordinate system or deformation. In contrast, discrete exterior calculus (DEC), through the definition of the exterior derivative operator $\mathbf{d}$, preserves this topological structure by construction, regardless of the coordinate representation. Panel (c) further illustrates this idea using a discrete network representation, where vertices are connected by weighted edges encoded by their colors. Conceptually, the information contained in the network can be separated into two parts: the connectivity of the graph (topology) and the edge weights (metric information). Within the framework of DEC, the operator $\mathbf{d}$ isolates precisely this topological information. To recover physically meaningful quantities, such as gradients per unit length or fluxes per unit area, the geometric information of the domain must be incorporated. This is achieved through the Hodge star operator, $\star$, which is metric-dependent and accounts for lengths, areas, and volumes in the system. The Hodge star operator plays a complementary role by connecting these different types of quantities through geometric information such as lengths, areas, and volumes. In particular, it maps a $k$-form into an $(n-k)$-form, allowing conversion between quantities defined on complementary geometric objects (for example, between edge-based and face-based quantities).

Figure 2: (a) A scalar function $f(x,y)$ and (b) its mapped representation under the coordinate deformation $x=X^{2}$ and $y=Y^{2}$. Although the scalar field itself remains unchanged, its gradient (or directional derivative along any given direction) depends on the local coordinate system. Nevertheless, the integral quantity $\oint df$ around any closed loop, such as those shown in the figures, always vanishes. This reflects a topological property of the connected vertices (white dots) that is independent of the local coordinate deformation. (c) Illustration of the same concept in discrete networks. The vertices are connected through weighted edges, with the weights encoded by the edge colors. Equivalently, the network can be represented by an unweighted connectivity graph (white edges), which encodes the topology, together with a separate list of edge weights. One of the main advantages of DEC over many other numerical schemes is that it naturally separates topology from geometry. As a result, identities and conservation laws of the original continuous problem, such as $\mathrm{curl}\,\mathrm{grad}\,f=0$, remain exactly preserved after discretization.

In DEC these ideas are implemented directly on a computational mesh. The exterior derivative, $\mathbf{d}$, is represented exactly using matrices that describe how mesh elements are connected—for example, which edges connect to which nodes, or which faces are bounded by which edges. These matrices act as discrete versions of derivatives, capturing how quantities change across the mesh. On the other hand, the Hodge star operator, which depends on geometry, is approximated using matrices that encode the size and shape of the mesh elements. These matrices enable consistent conversion between quantities defined on different parts of the mesh, ensuring that both the geometry and the underlying physical structure of the problem are properly captured. It is important to emphasize that, in the context of DEC, the separation between topological and metric operators serves a purpose far beyond mere mathematical elegance. As discussed earlier, this separation guarantees that key structural properties of the underlying continuous problem—such as the identity $\nabla\times\nabla f=0$—are preserved exactly at the discrete level. More generally, DEC ensures that the fundamental relation $d^{2}f=0$, which encompasses identities of this kind, is always satisfied on the discrete mesh of the computational domain. This stands in contrast to standard discretization schemes, such as finite differences (FD), finite elements (FEM), or spectral methods, where such identities are typically satisfied only approximately.

Mathematical Formulation of the EC-NLSE:— Having outlined the conceptual framework of EC, we now present the basic mathematical formulation for representing the NLSE in the EC language. However, since EC and DEC, along with their associated concepts such as forms and differential forms, are not widely familiar in the context of nonlinear optics, and to preserve the flow of the manuscript, we restrict ourselves here to introducing only the essential concepts, notations, and final formulations. A more detailed mathematical treatment is deferred to Supplementary Note 2. Before proceeding, we make three important remarks. First, we emphasize that EC is a framework rather than a rigid representation of the underlying PDE (the NLSE in our case). Consequently, there is more than one way to define the variables and to cast the PDE in the language of exterior calculus. Second, once a particular choice is made, it is essential to ensure that all terms in the resulting formulation correspond to the same $k$-form. In this sense, the requirement is analogous to enforcing consistency of physical units across all terms. Third, while the operators $\mathbf{d}$ and $\star$ are defined abstractly for differential forms of any degree, their explicit action depends on both the degree $k$ of the form and the dimensionality of the space, $n$. For example, applying $\mathbf{d}$ to a $0$-form (a scalar function defined at points) produces a $1$-form (quantities associated with edges), corresponding to the gradient operator. When applied to a $1$-form, $\mathbf{d}$ yields a $2$-form (quantities associated with oriented areas), corresponding to the curl operator. In general, the exterior derivative maps a $k$-form to a $(k+1)$-form. On the other hand, the Hodge star operator $\star$ maps a $k$-form to an $(n-k)$-form.

In this work, we begin by noting that the optical field is described by a scalar function $\psi(x,y,z)$, which we naturally treat as a $0$-form. Similarly, the scalar potential is also represented as a $0$-form. Other scalar quantities, such as $\chi(x,y)$ and $|\psi|^{2}$, are likewise treated as $0$-forms. Within the DEC framework, these quantities are therefore defined on the vertices of the primal mesh. Their multiplication is understood as pointwise multiplication, which, in the language of EC, corresponds to the wedge product $\wedge$ between $0$-forms. Next, as shown in Supplementary Note 2.1, the Laplacian operator can be expressed as $\star\mathbf{d}\star\mathbf{d}$. Taken together, these considerations lead to the following exterior calculus formulation of the NLSE:

$$i\frac{\partial\psi}{\partial z}=-\star\mathbf{d}\star\mathbf{d}\psi-V\wedge\,\psi-\,\chi\wedge|\psi|^{2}\wedge\,\psi.$$ (2)

In the above equation, we retain only the leading-order diffraction term and consider Kerr nonlinearity, both of which are relevant to many practical applications. The more general form of the equation, including higher-order terms, is presented in Supplementary Note 2.2. For clarity, we also set $p_{1}=q_{L}=q_{NL}=1$, although other normalizations can be adopted straightforwardly.

Note that every term in this equation is consistently a $0$-form, as expected. While this formulation is mathematically equivalent to the original equation, it is not ideally suited for numerical implementation. The primary reason is that the potential function may be discontinuous for sharply varying profiles, making its value less well-defined at boundaries. The exterior calculus formalism provides a systematic framework for addressing this problem by representing the linear potential and nonlinear response using pulse basis functions defined on the triangles, rather than in terms of delta functions at the nodes [10, 11]. In other words, the potentials are defined as a constant value per triangle. In particular, by applying the $\star$ operator to Eq. 2 from the left, and using the property that applying $\star$ twice yields the identity, we obtain:

$$i\star\frac{\partial\psi}{\partial z}=-\mathbf{d}\star\mathbf{d}\psi-\star V\wedge\,\psi-\star\,\chi\wedge|\psi|^{2}\wedge\,\psi.$$ (3)

In this form, every term in the equation is a $2$-form, which can be interpreted as representing an averaged quantity over the area associated with each vertex (see Supplementary Note 2.3 for details). This interpretation naturally introduces a form of local averaging, improving the numerical robustness of the formulation, particularly in the presence of sharp spatial variations. Next, since the quantities $V$, $\psi$, $\chi$, and $|\psi|^{2}$ are all $0$-forms, the above equation can be recast as (see Supplementary Note 2.4 for more details):

$$i\star\frac{\partial\psi}{\partial z}=-\mathbf{d}\star\mathbf{d}\psi-\star_{V}\psi-\star_{\mathcal{N}}\psi,$$ (4)

where $\star_{V}$ and $\star_{\mathcal{N}}$ are linear operators associated with the potential $V$ and the nonlinear interaction $\chi|\psi|^{2}$, respectively. The role of these operators will become clearer when we discuss the discrete formulation of this equation below.

DEC computational framework:— Next, we present the DEC computational framework employed in this work. As discussed earlier, different physical quantities are represented on different elements of the DEC computational mesh. To achieve this, DEC utilizes complementary primal and dual meshes [12, 13, 14, 15], as illustrated schematically in Fig. 3(a) for a two-dimensional computational domain.

Figure 3: (a) Geometric representation of a two-dimensional primal simplicial complex and its barycentric dual. The top panel shows a local mesh patch surrounding a central primal vertex, $v_{0}$. The shaded red region denotes the corresponding dual $2$-cell, constructed by connecting the barycenters of the neighboring primal $2$-simplices (e.g., $b_{012}$) through the barycenters (midpoints) of the shared primal $1$-simplices. The lower panel summarizes the duality mapping induced by the discrete Hodge star operator, $\ast$, which maps each primal $k$-simplex, $\sigma^{(k)}$, to a dual $(2-k)$-cell, $\ast\sigma^{(k)}$. In this two-dimensional complex, the $2$-simplices are assigned a consistent counterclockwise orientation determined by the ordered vertices $[v_{i},v_{j},v_{k}]$, whereas the orientation of the $1$-simplices (edges) is chosen arbitrarily. (b) Discretization of the potential and nonlinear terms on the primal simplicial complex. Each triangle $t_{k}=[v_{k_{1}},v_{k_{2}},v_{k_{3}}]$ carries a piecewise-constant refractive-index potential $V_{k}$ and a nonlinear contribution $\chi_{k}\,\left|\tilde{\psi}_{k}\right|^{2}$, where $\left|\tilde{\psi}_{k}\right|^{2}=\bigl(|\psi_{{k_{1}}}|^{2}+|\psi_{{k_{2}}}|^{2}+|\psi_{{k_{3}}}|^{2}\bigr)/3$ is the averaged intensity over the three vertices $[v_{k_{1}},v_{k_{2}},v_{k_{3}}]$ of $t_{k}$.

The blue mesh represents the primal mesh, whose vertices, edges, and areas correspond to 0-, 1-, and 2-simplices, respectively. These simplices host the primal (i.e., independent and not derived from one another) 0-, 1-, and 2-forms in the same order. The red mesh, on the other hand, represents the dual mesh, which similarly consists of its own 0-, 1-, and 2-cells. The connection between the primal and dual meshes is established through the Hodge star operator. For example, applying the Hodge star operator to a primal 0-form (i.e., a quantity defined on the blue vertices) produces a dual 2-form that resides on the corresponding dual 2-cell, represented by the red-shaded hexagonal region. This primal–dual pairing separates topology, encoded combinatorially by the exterior derivative on the primal mesh, from geometry, encoded metrically by the Hodge star through ratios of primal and dual cell volumes or via the Galerkin process.

More rigorously, the simplicial complex, denoted by $\mathcal{K}$, contains a natural hierarchy of simplices: 2-simplex faces (triangles), their 1-simplex edges, and their 0-simplex nodes. Each $k$-simplex $\sigma^{k}\in\mathcal{K}$, for $k=0,\ldots,2$, is defined by its $k+1$ vertices as $\sigma^{k}=[v_{1},\ldots,v_{k+1}]$, where the subscripts denote node indices [12]. The 0-simplices are the vertices $\{v_{1},\ldots,v_{N_{0}}\}$, where $N_{0}$ is the number of nodes. The 1-simplices are edges $\{e_{1},\ldots,e_{N_{1}}\}$, with each edge $e_{k}\equiv\sigma^{1}_{k}=[v_{k_{1}},v_{k_{2}}]$ oriented according to the ordering of its endpoints, and $N_{1}$ is the number of edges. The 2-simplices are triangles $\{t_{1},\ldots,t_{N_{2}}\}$, with each triangle $t_{k}\equiv\sigma^{2}_{k}=[v_{k_{1}},v_{k_{2}},v_{k_{3}}]$ defined by three nodes and bounded by three 1-simplex edges. The orientation of the 2-simplices is assumed to be consistent throughout the mesh (e.g., counterclockwise in Fig. 3) (a), a condition typically enforced by mesh generation tools. For lower-dimensional simplices ($k<1$), we adopt a canonical orientation: each edge $[v_{k_{1}},v_{k_{2}}]$ is oriented such that $k_{1}<k_{2}$. Additionally, the orientation of each primal k-simplex induces a compatible orientation on the corresponding dual $(n-k)$-cell [13].

Discrete exterior derivative:- Since the simulation domain is represented by a simplicial complex $\mathcal{K}$, the smooth differential forms must also be discretized. In DEC, this is achieved by integrating each smooth differential $k$-form over the corresponding $k$-simplices of the mesh. The resulting discrete spaces are denoted by $\mathcal{C}^{k}(\mathcal{K})$, which represent discrete $k$-forms defined on $\mathcal{K}$.

The discretization is carried out through the de Rham map $\mathcal{R}$. For example, the primal $0$-form $\psi$ and the $1$-form $\mathbf{d}\psi$ are discretized as:

$$\mathcal{R}(\psi)=\boldsymbol{\Psi}=\left[\psi(v_{1}),\ldots,\psi(v_{N_{0}})\right]^{\mathsf{T}}\in\mathcal{C}^{0}(\mathcal{K}),$$ (5)

and

$$\mathcal{R}(\mathbf{d}\psi)=\left[\int_{e_{1}}\mathbf{d}\psi,\ldots,\int_{e_{N_{1}}}\mathbf{d}\psi\right]^{\mathsf{T}}\in\mathcal{C}^{1}(\mathcal{K}),$$ (6)

where $\mathsf{T}$ denotes the transpose. While Eq. (5) simply evaluates $\psi$ at the mesh vertices, Eq. (6) requires the discrete exterior derivative, denoted by $\fourierbb{D}_{k}$, which maps discrete $k$-forms to discrete $(k+1)$-forms.

We note that the operator $\fourierbb{D}_{k}$ is purely combinatorial and depends only on the topology of the mesh. Moreover, it is represented by a sparse $(N_{k+1}\times N_{k})$ incidence matrix encoding the signed relations between $k$-simplices and $(k+1)$-simplices. Its construction follows directly from the generalized Stokes’ theorem,

$$\int_{\partial\mathcal{M}}\boldsymbol{\alpha}=\int_{\mathcal{M}}\mathbf{d}\boldsymbol{\alpha},$$ (7)

where $\boldsymbol{\alpha}$ is a $k$-form and $\mathcal{M}$ is an oriented $(k+1)$-dimensional manifold with boundary $\partial\mathcal{M}$.

For the case of the smooth $0$-form $\psi$, the exterior derivative $\mathbf{d}\psi$ is a $1$-form integrated along each oriented edge $e_{k}=[v_{k_{1}},v_{k_{2}}]$. Applying Stokes’ theorem gives

$$[\mathcal{R}(\mathbf{d}\psi)]_{k}=\int_{e_{k}}\mathbf{d}\psi=\psi(v_{k_{2}})-\psi(v_{k_{1}})=\sum_{\ell}[\fourierbb{D}_{0}]_{k\ell}\psi_{\ell},$$

where $\psi_{\ell}=\psi(v_{\ell})$ and $k=1,\ldots,N_{1}$. The matrix $\fourierbb{D}_{0}\in\mathbb{R}^{N_{1}\times N_{0}}$ therefore maps vertex-based quantities to edge-based quantities, with entries

$$[\fourierbb{D}_{0}]_{k\ell}=\begin{cases}+1,&\text{if node $\ell$ is the head of edge $k$},\\ -1,&\text{if node $\ell$ is the tail of edge $k$},\\ 0,&\text{otherwise}.\end{cases}$$

Higher-order discrete exterior derivatives, such as $\fourierbb{D}_{1}$ (and $\fourierbb{D}_{2}$ in three dimensions), are constructed analogously.

Discrete Hodge star:- Certain operations—notably the Hodge star, which maps $k$-forms to $(n-k)$-forms, where $n$ is the dimension of $\mathcal{K}$—are often conceptually associated with a dual complex $\ast\mathcal{K}$. In the specific case of Delaunay meshes [16], the dual complex is constructed using Voronoi cells [17], ensuring that primal and dual elements are mutually orthogonal. While the discrete exterior derivative $\fourierbb{D}_{k}$ is purely topological, the discrete Hodge star $\fourierbb{H}_{k}$ encodes the metric structure of the domain. In our 2D problem, the Hodge star operator maps primal $k$-forms to dual $(2-k)$-forms by relating the measure of a primal $k$-simplex to that of its dual $(2-k)$-cell.

For general unstructured meshes, the primal-dual mesh orthogonality is not preserved. To address this, we adopt the Galerkin formulation for the Hodge star using Whitney forms [18]. The mapping is defined via the variational formulation [11], which ensures numerical stability and consistency on general meshes without requiring the explicit geometric construction of an orthogonal dual mesh.

The Galerkin Hodge star operator is denoted here by $\fourierbb{H}_{k}$. In general, we only need the Galerkin approach for $\fourierbb{H}_{1}$, leading to a sparse non-diagonal matrix. The operators $\fourierbb{H}_{0}$ and $\fourierbb{H}_{2}$ can be obtained using the barycentric dual mesh, and they are diagonal matrices with positive entries [19]. This enables efficient explicit time integration or propagation, since the inverse of $\fourierbb{H}_{0}$ is trivial. To better explain how this procedure is applied to our problem, let us recall Eq. (4):

$$i\star\frac{\partial\psi}{\partial z}=-\mathbf{d}\star\mathbf{d}\psi-\star_{V}\psi-\star_{\mathcal{N}}\psi.$$

To discretize the above equation, four distinct discrete Hodge star operators need to be defined. They correspond to the $\star$ operator on the left-hand side of the equation, $\star$ operator on the right-hand side (which are different because they act on 0-forms and 1-forms, respectively), as well as the operators $\star_{V}$ and $\star_{\mathcal{N}}$. Their discrete representations are denoted by $\fourierbb{H}_{0}$, $\fourierbb{H}_{1}$, $\fourierbb{H}_{0}(V)$, and $\fourierbb{H}_{0}(\mathcal{N})$, respectively. The operators $\fourierbb{H}_{0}$ and $\fourierbb{H}_{1}$ are purely geometrical, i.e., only constructed via the mesh metrics (triangle areas or edge lengths). On the other hand, operators $\fourierbb{H}_{0}(V)$ and $\fourierbb{H}_{0}(\mathcal{N})$ encode the information about the linear and nonlinear potential, respectively. For the linear potential $V$, we assign to each triangle $t_{k}$ a value $V_{k}$ of the potential evaluated at the centroid of triangle $t_{k}$. The nonlinear response $\mathcal{N}$ is also defined as $\chi_{k}\,\left|\tilde{\psi}_{k}\right|^{2}$, where $\left|\tilde{\psi}_{k}\right|^{2}:=\bigl(|\psi_{{k_{1}}}|^{2}+|\psi_{{k_{2}}}|^{2}+|\psi_{{k_{3}}}|^{2}\bigr)/3$ is the averaged intensity over the three vertices $[v_{k_{1}},v_{k_{2}},v_{k_{3}}]$ of triangle $t_{k}$. These details are also presented in Fig 3(b) and its caption. The entries of the diagonal matrices $\fourierbb{H}_{0}$ are then given by [19]

$$\left[\fourierbb{H}_{0}\right]_{p,p}=\frac{1}{3}\sum_{l}\left|t_{p_{l}}\right|,\quad 1\leq p\leq N_{0},$$

where $\set{t_{p_{l}}}_{l=1}^{n_{p}}$ is the set of triangles $n_{p}$ connected to the node $p$ and $\left|t_{p_{l}}\right|$ is the area of the triangle $t_{p_{l}}$. The factor $1/3$ accounts for the portion of the triangle area that is added to the total effective dual area for node $p$. Similarly, the entries of the diagonal matrices $\fourierbb{H}_{0}(V)$ and $\fourierbb{H}_{0}(\chi\left|\psi\right|^{2})$ are given by

$$\left[\fourierbb{H}_{0}(V)\right]_{p,p}=\frac{1}{3}\sum_{l}V_{p_{l}}\left|t_{p_{l}}\right|,\quad 1\leq p\leq N_{0},$$

and

$$\left[\fourierbb{H}_{0}(\chi\left|\psi\right|^{2})\right]_{p,p}=\frac{1}{3}\sum_{l}\chi_{p_{l}}\left|\tilde{\psi}_{p_{l}}\right|^{2}\left|t_{p_{l}}\right|,\quad 1\leq p\leq N_{0},$$

Finally, the Hodge operator $\fourierbb{H}_{1}$ is constructed via Whitney 1-forms. Given a generic triangle $t$ with nodes denoted locally as $[v_{1},v_{2},v_{3}]$, a Whitney 0-form $\mathcal{W}^{(0)}_{\ell}$ is associated with each node $v_{\ell},\,\ell=1,\ldots,3$, and is given by

$$\mathcal{W}^{(0)}_{v_{\ell}}=\lambda_{\ell}(\boldsymbol{r}),$$

where $\lambda(\boldsymbol{r})$ is the barycentric coordinates for the triangle $t$. Since 0-forms are scalars, there is no distinction between the Whitney 0-forms and their proxy fields [20]. The Whitney 1-forms are associated with edges. The locally indexed edges of the triangle $t$ are $e_{\ell},\,\ell=1,\ldots,3$, where the edge $e_{\ell}$ is defined by the nodes forming it as $e_{l}=[v_{l_{1}},v_{l_{2}}]$. The Whitney 1-form associated with the edge $e_{\ell}$ is

$$\mathcal{W}^{(1)}_{e_{\ell}}=\lambda_{l_{1}}(\boldsymbol{r})\mathbf{d}\lambda_{l_{2}}(\boldsymbol{r})-\lambda_{l_{2}}(\boldsymbol{r})\mathbf{d}\lambda_{l_{1}}(\boldsymbol{r}),$$

while its proxy field $\mathbf{W}^{(1)}_{e_{\ell}}$ is

$$\mathbf{W}^{(1)}_{e_{\ell}}=\lambda_{l_{1}}(\boldsymbol{r})\nabla\lambda_{l_{2}}(\boldsymbol{r})-\lambda_{l_{2}}(\boldsymbol{r})\nabla\lambda_{l_{1}}(\boldsymbol{r}).$$ (8)

The entries of the Galerkin Hodge star $\fourierbb{H}_{1}$ are then given by [11]

$$\left[\fourierbb{H}_{1}\right]_{p,q}=\int_{\Omega}\xi\,\mathbf{W}^{(1)}_{e_{p}}\cdot\mathbf{W}^{(1)}_{e_{q}}d\Omega,\quad 1\leq p,q\leq N_{1}$$ (9)

DEC beam propagation:— Having outlined how to construct the key discrete operators needed for our work, we now write the DEC form of the NLES:

$$i\frac{\partial\mathbf{\Psi}}{\partial z}=\boldsymbol{\mathcal{L}}\mathbf{\Psi}+\boldsymbol{\mathcal{N}}(\left|\psi\right|^{2})\mathbf{\Psi},$$ (10)

where the operators $\boldsymbol{L}$ and $\boldsymbol{\mathcal{N}}(\left|\psi\right|^{2})$ are given, respectively, by

$$\boldsymbol{\mathcal{L}}=-\fourierbb{H}_{0}^{-1}\left[\fourierbb{D}_{0}^{\mathsf{T}}\fourierbb{H}_{1}\fourierbb{D}_{0}+\fourierbb{H}_{0}(V)\right],$$ (11)

and

$$\boldsymbol{\mathcal{N}}(\left|\psi\right|^{2})=-\fourierbb{H}_{0}^{-1}\fourierbb{H}_{0}(\chi|\mathbf{\psi}|^{2}).$$ (12)

A key computational advantage of the DEC formulation is that $\fourierbb{H}_{0}$ is diagonal, making the computation of $\fourierbb{H}_{0}^{-1}$ particularly inexpensive. Consequently, the resulting numerical scheme naturally lends itself to explicit time-integration and beam-propagation methods. An important implementation detail is that the nonlinear term requires only one sparse matrix-vector product per stage. First, the values of the $\left|\tilde{\psi_{i}}\right|^{2}$ at triangles $i=1..N_{2}$ are found at each step using a pre-computed projection operator $\mathbf{Q}$ that maps vertex-based values to element-based values by averaging.

$$\mathbf{Q}=\frac{1}{6}|\fourierbb{D}_{1}||\fourierbb{D}_{0}|.$$ (13)

Second, notice that the transpose $\mathbf{Q}^{\mathsf{T}}$ is an $(N_{0}\times N_{2})$ matrix that maps element-based values back to vertex-based values by accumulating contributions from all adjacent triangles, which is indeed the operation we need to accumulate dual area portions and element-wise potential values of the triangles that share a given node. So, in matrix notation, the discrete nonlinear operator in Eq. (12) is given by

$$\boldsymbol{\mathcal{N}}(\left|\psi\right|^{2})=-\fourierbb{H}_{0}^{-1}\text{diag}\big(\mathbf{Q}^{\mathsf{T}}\mathbf{D}\mathbf{Q}|\mathbf{\Psi}|^{2}\big)\mathbf{\Psi},$$ (14)

where $\mathbf{D}=\text{diag}(\chi_{i}|t_{i}|)$ is an $(N_{2}\times N_{2})$ diagonal matrix with entries $(\chi_{i}|t_{i}|)$ for each triangle $t_{i}$. Notice that the combined action of calculating the nonlinear potential and the construction of the Hodge star operator $\fourierbb{H}_{0}(\chi\left|\psi\right|^{2})$ is efficiently computed via the pre-computed sparse matrix $\mathbf{Q}^{\mathsf{T}}\mathbf{D}\mathbf{Q}$, given in Eq. (14).

For beam propagation of both methods, we employ the adaptive-step Runge–Kutta–Fehlberg (RK45) scheme [26], as implemented in scipy.integrate.solve_ivp. The RK45 method is a fifth-order accurate explicit integrator with an embedded fourth-order error estimator, allowing for automatic step-size control based on a user-specified tolerance. In all of our simulations, and for both the DEC and spectral methods, we chose a fixed relative tolerance of rtol = $10^{-9}$ and an absolute tolerance of atol = $10^{-11}$ for the RK45 solver. For well-conditioned triangular meshes, both $\boldsymbol{\mathcal{L}}$ and $\boldsymbol{\mathcal{N}}(\left|\psi\right|^{2})$ are sparse matrices with total non-zero entries $n_{nz}$, yielding a computational complexity of $\mathcal{O}(n_{nz})$ per stage. This linear scaling, combined with the adaptive step-size control, enables highly efficient propagation over long distances without sacrificing accuracy. This is significantly more efficient than split-step and spectral methods, which typically require multiple FFT operations per step (each scaling as $\mathcal{O}(n\log n)$ for $n$ grid points), or implicit FEM schemes, which necessitate solving a nonlinear system at each step via iterative methods such as Newton–Raphson, incurring substantial additional overhead.

Numerical Benchmarking and Performance:— We now evaluate the performance of the proposed DEC-BPM framework in terms of accuracy, convergence, and computational efficiency by benchmarking it against the gold-standard FFT-based algorithm. In particular, we compare it with the Spectral Method in the Interaction Picture (SIP), as detailed in Supplementary Note 3, which is widely regarded as one of the most computationally efficient methods for solving the nonlinear Schrödinger equation.

The geometry considered throughout this work is a standard step-index fiber. This choice is motivated both by its role as the canonical waveguide geometry in nonlinear optics applications and by the stringent numerical challenge posed by its discontinuous refractive-index profile, which is difficult to represent accurately using many conventional methods. The physical and computational parameters used in the simulations are summarized in Table 2. For computational efficiency, all simulations are performed using the standard normalized form of the NLSE (see Supplementary Note 1 for details).

To provide a comprehensive assessment across a broad range of conditions, we consider several benchmark scenarios of progressively increasing geometric and topological complexity:

1. Fundamental soliton:- A stationary nonlinear eigenmode is propagated over long distances to evaluate physical fidelity and the preservation of conserved quantities.

2. Off-axis elliptic Gaussian beam:- A non-eigenmode excitation with broken cylindrical symmetry is used to examine geometric robustness in the presence of a sharp core–cladding interface.

3. Elliptic optical vortex:- A beam carrying orbital angular momentum ($m=1$), centered on the fiber axis, is considered to assess the method’s ability to resolve higher-order modal distributions that interact strongly with the core–cladding interface.

To ensure a fair comparison, all simulations employ the same adaptive RK45 integrator (from SciPy) to integrate the equation of motion along the propagation direction. Under this common integration protocol, the adaptive step size is governed by the spectral radius of each transverse discretization (the largest eigenvalue magnitude of the matrix representing the linear differential operator) rather than by the propagation physics alone (i.e. geometry and initial condition). The reported wall-clock times therefore reflect both the accuracy achieved per degree of freedom and the stiffness that each spatial discretization presents to an explicit integrator, both of which are intrinsic properties of the underlying numerical scheme.

Table 2: Physical and computational parameters used in all test cases.

Parameter	Symbol	Value
Fiber Properties
Operating wavelength	$\lambda$	$1.55\,\mu\text{m}$
Core radius	$R_{0}$	$25.0\,\mu\text{m}$
Normalization constant along $z$	$L_{d}$	$7347.27\mu\text{m}$
Core linear index	$n_{\text{core}}$	$1.450$
Cladding linear index	$n_{\text{clad}}$	$1.449$
Core nonlinear index	$n_{2,\text{core}}$	$2.8\times 10^{-20}\,\text{m}^{2}/\text{W}$
Cladding nonlinear index	$n_{2,\text{clad}}$	$2.2\times 10^{-20}\,\text{m}^{2}/\text{W}$
Computational Domain (normalized to $R_{0}$)
SIP	$[-4,4]\times[-4,4]$
DEC	$\{(x,y)\in\mathbb{R}^{2}\mid\sqrt{x^{2}+y^{2}}\leq 4\}$

Crucially, to validate and benchmark the proposed method, reliable reference solutions must first be established. For stable soliton propagation, the analytical soliton solution provides a natural benchmark. In scenarios where analytical solutions are unavailable, reference solutions are instead generated numerically using highly resolved simulations with the SIP method, the gold-standard FFT-based solver for the nonlinear Schrödinger equation. Specifically, the reference solutions are constructed as follows:

• Soliton Case: We compute the fundamental soliton eigenmode on a dense ($1024\times 1024$) grid using a self-consistent nonlinear eigenmode solver, as described in the Supplementary Note 4. Because this is a stationary self-consistent solution with an invariant intensity profile, the reference intensity is simply the input beam profile.

• Elliptic Gaussian and Vortex Cases: For these two examples, the nonlinear Schrödinger equation is propagated using the SIP method on a dense ($512\times 512)$ grid, and the resulting field distribution at the target propagation distance is taken as the reference solution. To keep the computational cost of generating these reference solutions manageable, and because of the SIP method’s inability to accurately handle radiated fields over longer propagation distances, we restrict the benchmark to a maximum normalized propagation distance of ($z=0.1$).

We will refer to the reference solutions discussed above as $\psi_{\mathrm{ref}}(\mathbf{r},z)$ and we will use the relative $L_{2}$-norm as a measure of the error:

$$E_{\mathrm{rel}}(z)=\frac{\lVert|\psi(\mathbf{r},z)|^{2}-|\psi_{\mathrm{ref}}(\mathbf{r},z)|^{2}\rVert_{2}}{\lVert|\psi_{\mathrm{ref}}(\mathbf{r},z)|^{2}\rVert_{2}},$$ (15)

Importantly, we note that all reported wall-clock times, including mesh generation and initialization overheads, were measured on a 32-core AMD Ryzen Threadripper Pro workstation equipped with 256 GB of RAM.

Fundamental soliton:- As mentioned earlier, our first example considers the propagation of a soliton in a standard optical fiber geometry (see Table 2). Since the soliton is a stationary solution of the NLSE, any distortion of its shape during propagation provides a direct measure of numerical error. This example therefore serves as a stringent test of the DEC method.

Figure 4: Initial intensity profile of soliton solutions obtained via the Spectral and DEC methods.

The input soliton profiles used in the SIP and DEC simulations are shown in Fig. 4. These profiles were obtained numerically using the self-consistent method described in Supplementary Note 4. Each soliton is computed self-consistently on the same grid or mesh subsequently used for its propagation, so that the launch field is a discrete eigenmode of the operator that evolves it. Figure S1 in the Supplementary Note 5 compares the transverse propagation of the soliton simulated via the SIP and DEC methods across varying spatial resolutions. Figure 5(a) presents the convergence of the computed soliton eigenvalue $\mu$ as a function of the number of degrees of freedom $N$ for both methods.

As a reference, we compute the soliton eigenvalue $\mu_{\mathrm{ref}}=-3.5537$ using the SIP solver on a $1024\times 1024$ reference grid. The two solvers are mutually consistent at the $10^{-3}$ level, with the DEC solver reaching this agreement using approximately an order of magnitude fewer unknowns ($\mu_{\mathrm{DEC}}=-3.5540$ at $N=25{,}217$ versus $\mu_{\mathrm{SIP}}=-3.5547$ at $N=262{,}144$). This discrepancy in performance is primarily attributed to the sharp material discontinuity introduced by the step-index fiber profile, which induces Gibbs oscillations in Fourier-based spectral methods and consequently slows their convergence. By contrast, the second-order DEC formulation naturally accommodates local discontinuities without introducing global oscillatory artifacts, enabling it to accurately capture the step-index profile while achieving substantially higher computational efficiency.

Figure 5: Fundamental soliton propagation in a step-index fiber. (a) Convergence of the soliton eigenvalue $\mu$ versus degrees of freedom $N$. The DEC method (blue circles) reaches agreement with the reference value using $\approx 10\times$ fewer unknowns than the SIP method (red squares), whose convergence is slowed by the staircase representation of the circular core–cladding interface on the uniform Cartesian grid. (b) Long-term shape-preservation error $E_{\text{rel}}$ over a propagation distance $z=25$. SIP fluctuates in the $10^{-3}$–$10^{-2}$, whereas DEC maintains fidelity near the numerical noise floor ($\sim 10^{-6}$).

To study the propagation dynamics, we consider a normalized propagation distance of $z=25$, corresponding to a physical propagation length of $25L_{d}=18.37$ cm. The long-term stability advantage of the DEC approach is illustrated in Fig. 5(b), which shows the evolution of the relative error $E_{\text{rel}}$ as given by Eq. (15) over the entire propagation distance. Throughout the simulation, the error associated with the SIP approach oscillates between $10^{-3}$ and $10^{-2}$. In contrast, the DEC method preserves the soliton profile with errors remaining near the numerical noise floor ($\sim 10^{-6}$) over the entire propagation range, demonstrating superior long-term stability and shape-preserving fidelity.

Off-axis elliptic Gaussian beam:- The second test case is specifically designed to stress-test the ability of each numerical method to accurately handle sharp material discontinuities in the presence of a dynamically evolving and spatially asymmetric optical field. To this end, we launch an off-axis elliptic Gaussian beam. Since this input profile is not an eigenmode of the step-index waveguide, it undergoes pronounced beam breathing and self-phase modulation during propagation, causing the optical energy to repeatedly interact with the core–cladding interface. The initial field envelope is given by:

$$\psi_{0}(U,V)=A\exp\!\left(-\frac{U^{2}}{w_{x}^{2}}-\frac{V^{2}}{w_{y}^{2}}\right),$$ (16)

where $(U,V)$ represent the rotated transverse coordinates centered at the beam’s offset position $(x_{0},y_{0})$ and inclined by an azimuthal angle $\theta$. These are defined through the spatial transformation:

$$U =(x-x_{0})\cos\theta+(y-y_{0})\sin\theta,$$ (17)

$$V =-(x-x_{0})\sin\theta+(y-y_{0})\cos\theta.$$

In our simulations, we used the following parameters: amplitude $A=1.0$, asymmetric beam waists $w_{x}=1.2$ and $w_{y}=0.6$, tilt angle $\theta=\pi/4$, and an off-axis shift of $(x_{0},y_{0})=(0.1,-0.1)$. The resulting transverse intensity profiles, $|\psi(x,y,z=0)|^{2}$, initialized on both the Cartesian spectral grid and the unstructured DEC mesh, are shown in Fig. 6. By intentionally breaking the cylindrical symmetry of the waveguide, this configuration provides a stringent test of the unstructured-mesh capabilities of the DEC solver.

Figure 6: Initial transverse intensity profiles ($|\psi|^{2}$) of the off-axis elliptic Gaussian beam in normalized units. (a) Initialization on the Cartesian grid for the SIP method. (b) Initialization on the unstructured circular domain for the DEC method. The beam is shifted off-center to $(x_{0},y_{0})=(0.1,-0.1)$, rotated by an angle of $\theta=\pi/4$, and exhibits an elliptical aspect ratio of $w_{x}/w_{y}=2$. The dashed white circle designates the core–cladding interface ($r=R_{0}$) of the step-index fiber.

Figure 7: Performance benchmark for the off-axis elliptic Gaussian beam in a step-index fiber. (a) Relative $L_{2}$ error versus degrees of freedom $N$. The DEC method (blue circles) maintains second-order convergence $\mathcal{O}(h^{2})$, corresponding to the observed $\mathcal{O}(N^{-1})$ scaling. Conversely, the SIP method (red squares) is limited to first-order $\mathcal{O}(h^{1})$ convergence (observed here as $\mathcal{O}(N^{-0.5})$) by the Gibbs phenomenon at the refractive-index discontinuity. (b) Error versus compute time. At the representative accuracy threshold $E_{\text{rel}}=5\times 10^{-3}$ (dashed magenta line), DEC requires $\approx 1.19$ s versus $\approx 87.15$ s for SIP, a speedup of 73x.

Figure 8: Visual comparison of the computed transverse intensity profiles ($|\psi|^{2}$) for the off-axis elliptic Gaussian beam at the dimensionless propagation distance of $z=0.1$. (Top row, a–d) Solutions obtained using the SIP method on a square Cartesian grid, illustrating the grid density required to control Gibbs ringing at the refractive-index discontinuity. (Bottom row, e–h) Corresponding solutions using the DEC method on an unstructured circular domain. Each subplot is annotated with the total number of degrees of freedom ($N$) and the required compute time ($t$) to reach $z=0.1$. The dashed white circle denotes the core-cladding interface of the step-index fiber. The DEC method achieves a comparable high-resolution field representation to the spectral approach but with substantially fewer unknowns, highlighting its superior computational efficiency for discontinuous refractive-index profiles.

To compare the performance of the SIP and DEC methods, we plot the convergence behavior and computational cost of both approaches for this input beam in Fig. 7(a)–(b), respectively. The DEC formulation consistently outperforms the SIP solver, requiring a factor of 6.2 fewer degrees of freedom to achieve comparable accuracy and yielding a computational speedup of approximately 73 times. More specifically, the DEC method exhibits clear second-order convergence, $\mathcal{O}(h^{2})$, corresponding to an $\mathcal{O}(N^{-1})$ scaling with respect to the number of degrees of freedom $N$. This behavior is confirmed numerically as the DEC mesh is refined from $N=409$ to $25{,}217$ nodes, for which the measured convergence order steadily improves from 2.04 to 2.15, while the relative error decreases to $1.09\times 10^{-3}$. In contrast, the SIP method suffers from Gibbs oscillations induced by the sharp refractive-index discontinuity at the core–cladding interface, which degrades its theoretical exponential convergence to an effective first-order behavior, $\mathcal{O}(h^{1})$  (observed numerically as $\mathcal{O}(N^{-0.5})$), consistent with the $L_{2}$ representation error of a discontinuous coefficient sampled on a uniform grid. As a result, the spectral solver remains severely bottlenecked across the tested resolutions, causing its computational cost to increase rapidly with accuracy demands. For instance, achieving a representative relative error threshold of $5\times 10^{-3}$ requires approximately 87.15 seconds using the SIP solver, whereas the DEC method—which is naturally insensitive to local discontinuities—reaches the same accuracy in only 1.19 seconds. This corresponds to a $73$-fold speedup under the common-integrator protocol, highlighting both the efficiency and robustness of the DEC framework for modeling optical fibers with step-index profiles.

The qualitative behavior underlying the convergence results is illustrated in Fig. 8, which shows the transverse intensity profiles, $|\psi|^{2}$, of the input beam after propagation to the normalized distance $z=0.1$ for different numbers of degrees of freedom (grid points for the SIP method and mesh elements for the DEC method). At this stage of propagation, the beam has already begun to interact strongly with the step-index interface while undergoing nonlinear reshaping. The top row corresponds to the SIP solution on a Cartesian grid, whereas the lower row shows the DEC solution on an unstructured circular mesh conforming to the fiber geometry. Figure S2 in Supplementary Note 5 compares the transverse propagation dynamics of the elliptic Gaussian beam simulated via the SIP and DEC methods across increasing spatial resolutions.

Although both methods capture the overall confinement of the optical field within the fiber core (marked by the dashed circle), the SIP solution exhibits pronounced Gibbs ringing near the sharp refractive-index discontinuity unless extremely fine grids are employed. Consequently, achieving acceptable accuracy with the SIP solver requires a substantial increase in computational resolution and runtime, even for relatively short propagation distances. In contrast, the DEC discretization naturally conforms to the material interface and therefore produces smooth, artifact-free field distributions using only a small fraction of the degrees of freedom. This geometric adaptability enables the DEC framework to achieve high-fidelity solutions with significantly reduced computational cost.

Elliptic optical vortex:- Next, we consider the third example: an elliptic optical vortex beam carrying orbital angular momentum (OAM) with topological charge $m=1$. To construct the elliptic vortex beam, we again employ the rotated coordinate system defined in Eq. (17), now centered on the fiber axis, $(x_{0},y_{0})=(0,0)$, with rotation angle $\theta=\pi/4$. The input field profile is given by:

$$\psi_{0}(U,V)=A\rho(U,V)e^{i\Phi(U,V)}\exp\!\left(-\frac{U^{2}}{w_{x}^{2}}-\frac{V^{2}}{w_{y}^{2}}\right),$$ (18)

where the complex prefactor has been expressed in polar form with an elliptical amplitude $\rho(U,V)=\sqrt{(U/w_{x})^{2}+(V/w_{y})^{2}}$ and a spatially varying phase $\Phi(U,V)=\arg\left(\frac{U}{w_{x}}+i\frac{V}{w_{y}}\right)$. This representation explicitly highlights the characteristic phase winding associated with the vortex structure. Here, we set $A=1.0$, $w_{x}=1.0$ and $w_{y}=0.8$. Figures 9(a) and (b) illustrate the transverse intensity profile of this beam on both the Cartesian spectral grid and the unstructured DEC mesh prior to propagation.

Figure 9: Initial transverse intensity profiles ($|\psi|^{2}$) of the elliptic optical vortex beam in normalized units. (a) Initialization on the Cartesian grid for the SIP method. (b) Initialization on the unstructured circular domain for the DEC method. The beam is centered at the origin $(0,0)$, rotated by $\theta=\pi/4$, and carries an orbital angular momentum (OAM) with topological charge $m=1$.

Having defined the input profile, we next propagate the beam in both solvers over a normalized distance of $z=0.1$. Figure 10(a) and (b) compare the convergence behavior and computational efficiency of the two numerical schemes. From the plots, it is evident that the DEC method maintains a highly stable second-order convergence rate, $\mathcal{O}(h^{2})$, corresponding to an $\mathcal{O}(N^{-1})$ scaling with respect to the number of degrees of freedom. In particular, as the DEC mesh is refined from $N=409$ to $25{,}217$ nodes, the measured convergence order remains consistently between 2.02 and 2.17, while the relative error decreases to $1.56\times 10^{-3}$ in only 18.28 seconds. Conversely, the SIP method initially exhibits a noticeably shallow convergence trajectory at lower resolutions. Although its convergence rate improves as the spectral grid is further refined, the method remains significantly less efficient overall. At its highest resolution of $N=65{,}536$, the spectral solver takes nearly 269 seconds only to reach a relative error of $7.06\times 10^{-3}$. Figure 10(b) further highlights this disparity at the accuracy threshold $E_{\text{rel}}=8\times 10^{-3}$: the DEC method reaches this target in $\approx 1.2$ seconds, whereas interpolation along the SIP curve—whose finest point, $7.06\times 10^{-3}$ in $269$ seconds, already exceeds this accuracy—gives $\approx 181$ seconds at the threshold, a speedup of $\approx 150\times$ under the common-integrator protocol.

Figure 10: Performance benchmark for the elliptic optical vortex ($m=1$). (a) Relative $L_{2}$ error versus degrees of freedom $N$. DEC (blue circles) maintains second-order convergence $\mathcal{O}(h^{2})$ that corresponds to the observed $\mathcal{O}(N^{-1})$ scaling, while the SIP error (red squares) decreases non-monotonically, consistent with the resolution-dependent staircase bias of the discrete interface. (b) Error versus compute time. At the accuracy threshold $E_{\text{rel}}=8\times 10^{-3}$ (dashed line), DEC requires $\approx 1.2$ s; interpolation of the measured SIP curve gives a speedup of $\approx 150\times$ under the common RK45 protocol.

The widening of the performance gap relative to the Gaussian case calls for care in attribution. The vortex field itself is smooth: $\psi_{0}\propto(U/w_{x}+iV/w_{y})\exp(-U^{2}/w_{x}^{2}-V^{2}/w_{y}^{2})$ is an entire function of the transverse coordinates, and the on-axis phase singularity coincides with an intensity null at which the field is locally linear—neither discretization has intrinsic difficulty representing it. The enhanced gap instead originates at the material interface. First, the $m=1$ launch projects predominantly onto odd-azimuthal guided modes, including the near-cutoff LP₃₁ group (cutoff $V=5.14$ against the fiber’s $V=5.46$), whose fields concentrate at and beyond the core–cladding boundary; the geometric error of the staircased interface is therefore sampled precisely where the field resides. Second, the vortex launch has poorer overlap with the guided-mode set and sheds a larger radiated fraction than the Gaussian; on the uniform spectral grid, this halo must be represented at full resolution over the entire rectangular domain, whereas the graded DEC mesh resolves the core finely and the outer region coarsely. Third, the $\mathcal{O}(h)$ variation of the staircased core radius produces the non-monotone SIP convergence noted above.

Figure 11: Visual comparison of the computed transverse intensity profiles ($|\psi|^{2}$) for the elliptic optical vortex ($m=1$) at the dimensionless propagation distance of $z=0.1$. (Top row, a–d) Solutions obtained using the SIP method on a square Cartesian grid, illustrating the grid density required to control Gibbs errors at the step-index interface. (Bottom row, e–h) Corresponding solutions using the DEC method on an unstructured circular domain. Each subplot is annotated with the total number of degrees of freedom ($N$) and the required compute time ($t$) to reach $z=0.1$. The DEC method effectively preserves the characteristic doughnut-shaped intensity profile with substantially fewer unknowns.

Figure 11 presents the output intensity profiles at the normalized propagation distance $z=0.1$ for different grid and mesh resolutions. The top row corresponds to the SIP method on a Cartesian grid, whereas the bottom row shows the DEC solution on the unstructured circular mesh. Evidently, the spectral approach requires fine grids and long computation times to suppress Gibbs errors at the step-index interface, which for this launch is strongly illuminated by near-cutoff odd-azimuthal mode content. As a consequence, preserving the characteristic doughnut-shaped intensity profile of the $m=1$ vortex becomes computationally demanding for the SIP method. By contrast, the DEC formulation naturally conforms to the material boundary and maintains a high-fidelity representation of the vortex ring using only a small fraction of the degrees of freedom and computational cost.

As a final remark, we emphasize that the propagation distances considered in this work are not an inherent limitation of the DEC-BPM itself, but rather of the current implementation, which employs Neumann boundary conditions. These distances can be extended straightforwardly by incorporating absorbing boundary conditions. To illustrate this point, we performed simulations over substantially longer propagation distances by enlarging the computational domain. For centrally launched fields, this incurs only a modest additional cost. In Supplementary Note 6, we extend the Gaussian beam propagation to five times the original distance using a computational window only 2.5 times as large, while maintaining essentially the same accuracy and relative performance reported in the main text. In contrast, for annular fields launched near the core–cladding interface, such as vortex beams, radiated energy reaches the domain boundary much sooner. Consequently, maintaining the same level of accuracy requires a much larger computational window, making long-distance simulations increasingly uneconomical without the use of absorbing boundary conditions.

III Conclusion

In this work, we introduced a topology-preserving computational framework based on DEC for the efficient simulation of nonlinear wave propagation in complex media. By replacing the continuous domain with a simplicial complex and representing field quantities as discrete differential forms, the governing equations emerge directly from the discrete exterior derivative and Hodge star, without requiring the variational formulation and element-level matrix assembly characteristic of conventional finite-element methods, while at the same time avoiding the geometric restrictions of FFT-based spectral solvers. Unlike Cartesian discretizations, the proposed framework conforms naturally to complex geometries and material interfaces while preserving the topological structure of the underlying differential operators.

Our results demonstrate that this structure-preserving formulation enables accurate and computationally efficient simulations of multiscale nonlinear photonic systems containing sharp refractive-index discontinuities. In these challenging regimes, DEC-BPM resolves complex geometries with substantially fewer degrees of freedom, achieving spectral-level accuracy while delivering computational speedups exceeding two orders of magnitude over conventional spectral methods. Beyond improving computational performance, these gains make routine simulations of previously inaccessible large-scale nonlinear photonic systems computationally practical, reducing runtimes from days to hours and enabling numerical optimization beyond what is currently feasible using simplified models or physical intuition alone.

Another important advantage of the DEC-BPM scheme presented here is its natural extensibility to rigorous absorbing boundary conditions, such as perfectly matched layers (PMLs). In particular, the geometric structure of the discrete operators in the DEC formulation naturally accommodates PMLs through complex stretching of the discrete Hodge stars. Since this stretching modifies the discrete metric rather than any individual operator, the same construction carries over directly to full-vector formulations, in the spirit of exterior complex scaling [21]. Incorporating such absorbers is therefore a natural next step toward accurate long-distance propagation in these more general settings. By contrast, FFT-based algorithms do not enjoy this flexibility. Their efficiency relies on the Fourier basis diagonalizing the constant-coefficient Laplacian, a property that is destroyed by the spatially varying coefficients introduced by a PML transformation. Although absorbing boundary conditions can be incorporated into FFT-based methods through approaches such as complex absorbing potentials [22, 23, 24, 25], these techniques are generally more ad hoc and lack the same geometric and mathematical foundation as PMLs. Consequently, they do not provide comparable guarantees of convergence or a systematic path toward extension to more complex settings, such as waveguides with arbitrary geometries, asymmetric refractive-index profiles, or fully vectorial electromagnetic fields.

More broadly, the present framework establishes discrete exterior calculus as a general computational paradigm for nonlinear wave propagation. In addition to its demonstrated computational efficiency, the structure-preserving formulation naturally accommodates geometry-conforming discretizations, and the systematic construction of higher-order propagation operators from lower-order discrete operators. These capabilities provide a natural pathway toward more general models of nonlinear wave propagation, including ultrashort pulse propagation, vectorial and polarization-dependent formulations, and higher-order nonlinear effects in complex photonic structures. Moreover, because the discrete operators developed here constitute the same algebraic building blocks required for more general electromagnetic models, the framework can be extended systematically while retaining its geometric and topological structure. We anticipate that these capabilities will enable the simulation, optimization, and inverse design of increasingly sophisticated nonlinear photonic devices, including photonic crystal fibres, multicore fibres, multimode waveguides, and integrated nonlinear photonic platforms.