1. Introduction
Nonlinear dynamical lattices are fundamental in mathematics, physics, and engineering, particularly in Hamiltonian systems. Atomic vibrations in molecular and crystalline structures are often analyzed using the harmonic approximation, which assumes small displacements and retains only quadratic terms [1–5]. However, in solid crystals with on-site potentials, discrete translational symmetry is broken, leading to the loss of total momentum conservation, an effect that becomes particularly significant in higher-dimensional atomic lattices. Consequently, on-site potentials influence anomalous thermal transport in monoatomic lattices with intersite nonlinearities, revealing insights into energy localization and wave propagation. Nonlinear lattice systems are known to sustain localized excitations, such as soliton-like behavior [6–9]. Studies on one-dimensional lattices, beginning with the Fermi–Pasta–Ulam models, have revealed key soliton dynamics. While defining solitons in two-dimensional systems is challenging, intrinsic localized modes at supersonic velocities have been observed, with theoretical predictions confirmed in fluid mechanics, optics, acoustics, and nonlinear lattices [10–15].
Among widely studied nonlinear wave models, the Davey–Stewartson (D–S) equation and its generalizations provide fundamental descriptions of (2+1)-dimensional wave packets in dispersive and nonlinear media. Most studies treat the D–S equation as a continuum model with phenomenological coefficients. In real lattice systems, however, these coefficients emerge from microscopic interatomic interactions [16, 17]. Recent work has derived effective multidimensional equations from Hamiltonian lattices [18, 19], producing physically interpretable parameters that allow tunable soliton behavior. Cubic and on-site potentials critically shape soliton dynamics in nonlinear optical lattices. On-site potentials affect the stability and mobility of localized waves, while cubic nonlinearities govern soliton formation and interactions, especially in higher-dimensional lattices [20]. Evaluating how these potentials influence soliton structure, collision outcomes, and modulational instability (MI) is crucial for designing nonlinear photonic devices, optical communication systems, and engineered lattices that enable controlled energy transport.
Motivated by these developments, we derive a generalized (2+1)-D equation from a monoatomic lattice with cubic intersite and nonlinear on-site interactions. In optical lattices with self-repulsive cubic nonlinearities, multiple soliton families, such as fundamental, dipole, and multi-component modes, can emerge across regular as well as defect-engineered configurations [21]. Overall, the study demonstrates how lattice-induced cubic and on-site potentials regulate soliton dynamics, influence interaction patterns, and enable controlled energy localization in complex nonlinear systems [22]. While numerous studies analyze the D–S equation using phenomenological parameters, very few works derive a (2+1)-dimensional system directly from a discrete Hamiltonian lattice containing both cubic intersite interactions and nonlinear on-site potentials. Unlike previous DS analyses, all coefficients in our model originate from the lattice Hamiltonian, allowing us to reveal how specific physical interactions govern multidimensional soliton dynamics. The lattice driven D–S framework reveals mechanisms that are missing in continuum models, shedding new light on multidimensional nonlinear transport in photonic and atomic lattices [23–25].
This paper examines the dynamics of a generalized (2+1)D D–S-type nonlinear evolution equation. In reference [26], the authors investigated a related nonlinear soliton equation within an analytical framework. In the present work, we employ the Hirota bilinear method to construct collisional soliton solutions of the (2+1)-dimensional D–S equation. In particular, we derive a double-soliton solution to explore the dynamical behavior of nonlinear optical lattices. The manuscript is structured as follows: section 2 develops the governing equations using a quasi-discrete multiple-scale approach, followed by section 3 , which derives one- and two-soliton solutions via the Hirota bilinear method.
2. Formulation and dynamics of (2+1)-dimensional D–S equations in multidimensional lattices
Let us consider a monoatomic lattice in $N$ dimensions, consisting of particles each having mass $M$, where each particle interacts harmonically with its nearest neighbors and is subject to both cubic and on-site potential forces. Assuming small vibrational displacements, the interaction potential can be approximated using a Taylor series. The derivation proceeds in three steps. First, we formulate the Hamiltonian, including harmonic, cubic, and on-site potentials. Second, we obtain the equations of motion and apply a quasi-discrete multi-scale perturbation to separate fast oscillations from the slow envelope. Third, we extract the envelope dynamics and arrive at an effective (2+1)-dimensional DS-type evolution equation. The structured reduction directly associates each nonlinear coefficient in equation (3 ) with its underlying lattice interaction. This leads to a system described by the following Hamiltonian [27]:
$\begin{align} H & = \sum_{n} \Bigg[ \frac{1}{2} M \dot{u}(n)^{2} + \sum_{j = 1}^{d} \Bigg\{ \frac{1}{2} k_{2j} \big(u(n+a_{j}) - u(n)\big)^{2} \nonumber\\ & \quad + \frac{1}{3} k_{3j} \big(u(n+a_{j}) - u(n)\big)^{3} + \frac{1}{4} k_{4j} \big(u(n+a_{j}) - u(n)\big)^{4} \Bigg\} \nonumber\\ & \quad + \frac{1}{2} k_{5} u(n)^{2} + \frac{1}{3} k_{6} u(n)^{3} + \frac{1}{4} k_{7} u(n)^{4} \Bigg]. \end{align}$
Define $u(n)$ as the shift of an atom from its rest position at the lattice site $\sum_{j = 1}^{d}n_{j}a_{j}$, where $n_j$ are integers and $a_j$ are the basis vectors defining a $N$-dimensional lattice structure. The interaction coefficients $k_{r,j}$ for ($r = 0, 1, 2, 3..$) describe both harmonic and anharmonic forces between atoms. The term $k_{5,j}$ corresponds to the frequency of oscillation due to the external substrate potential, while $k_{6,j}$, $k_{7,j}$ account for on-site anharmonic contributions. The model focuses on a monoatomic lattice exhibiting interactions between adjacent atoms. The following equations of motion dictate the dynamics of the lattice, $\begin{align} M\,\frac{\mathrm{d}^{2}u(n)}{\mathrm{d}t^{2}}& = \sum_{j = 1}^{d} \Big[ k_{2j}\big(u(n+a_{j}) + u(n-a_{j}) - 2u(n)\big) \nonumber\\ & \quad+ k_{3j}\Big( \big(u(n+a_{j}) - u(n)\big)^{2} - \big(u(n) - u(n-a_{j})\big)^{2} \Big) \nonumber\\ & \quad+ k_{4j}\Big( \big(u(n+a_{j}) - u(n)\big)^{3} + \big(u(n-a_{j}) - u(n)\big)^{3} \Big) \Big] \nonumber\\ &\quad - k_{5}\,u(n) - k_{6}\,u(n)^{2} - k_{7}\,u(n)^{3}. \end{align}$
To keep the notation consistent between the microscopic Hamiltonian and the reduced envelope model, we explicitly connect the lattice force constants to the effective coefficients introduced later. The coefficients $k_{2j}$, $k_{3j}$, and $k_{4j}$ represent the harmonic, cubic, and quartic intersite interaction strengths in the $j\mathrm{th}$ lattice direction, while $k_5$, $k_6$, and $k_7$ are the on-site harmonic and anharmonic potentials. After performing the quasi-discrete multiple-scale expansion and averaging over the lattice directions, the above microscopic coefficients are reduced to the effective coupling constants. $\begin{eqnarray*} J_{2}\langle \equiv k_{2j}\rangle , J_{3}\equiv \langle k_{3j}\rangle , J_{4}\equiv\langle k_{4j}\rangle , J_{6}\equiv k_{6}, J_{7}\equiv k_{7}.\end{eqnarray*}$
Therefore, the $J$-parameters in the D–S equation are the renormalized parameters that are directly obtained from the Hamiltonian coefficients. Here, $J_{2}$ governs linear dispersion, $J_{3}$, $J_{4}$ represent cubic and quartic intersite nonlinearities, while $J_{6}$, $J_{7}$ correspond to cubic and quartic on-site nonlinear potentials. These equations incorporate the average dynamics driven by the oscillating wave packet via cubic interatomic interactions, providing deeper insights into nonlinear lattice dynamics.In their study, Huang and colleagues examined the nonlinear behavior of a lattice using a quasi-discrete multi-scale perturbation method [28]. This research emphasizes the effects of cubic and on-site potentials, underscoring their substantial influence on lattice dynamics. Through this analysis, a system of (2+1)-D nonlinear wave equations is developed to precisely characterize the modulation behavior of wave packets within the optical lattice, as shown below, 3 ) are described below 3 ).
$\begin{eqnarray} \mathrm{i}B_{t}+R_{1}B_{xx}+R_{2}B_{yy}+R_{3}B_{xy}+R_{4}|B|^2B = 0.\end{eqnarray}$
Unlike the standard D–S equation for continuous media, this formulation is derived directly from a discrete Hamiltonian lattice incorporating cubic intersite and on-site potentials. This tunable D–S system has coefficients $R_1$-$R_4$ determined by the lattice couplings $J_2$, $J_3$, and $J_6$, enabling microscopic control of soliton amplitude, width, and orientation. Hence, this study extends the classical D–S framework by linking discrete lattice physics with higher-dimensional nonlinear wave propagation. Here, $B$($x$, $y$, $t$) represents the slowly varying complex amplitude of the oscillating lattice wave. The parameters used in equation ( $\begin{align} R_{1} = &\frac{1}{2\omega}[J_{2}(\lambda_{1}^{2}\cos q_{1})+J_{2}(\lambda_{2}^{2}\cos q_{2})-v_{g}^2],\nonumber\\ R_{2} = &\frac{1}{2\omega}[J_{2}(\lambda_{1}^{2}\cos q_{1})+J_{2}(\lambda_{2}^{2}\cos q_{2})],\nonumber\\ R_{3} = &\frac{J_{2}}{2\omega}\lambda_{1}\lambda_{2}[\cos q_{2}-\cos q_{1}],\nonumber\\ R_{4} = &\frac{1}{2m}[8\mathrm{i}J_{3}(i\alpha_{1}+\alpha_{2})\nonumber\\ & \times [(1-\cos q_{1})\sin q_{1}+(1-\cos q_{2}) \sin q_{2}],\nonumber\\ \alpha_{1} = &\frac{4J_{2}[\sin q_{1}(\cos q_{1}-1)+\sin q_{2}(\cos q_{2}-1)]}{4[\omega(q)]^2-[\omega(2q)]^2}, \end{align}$
$\begin{align*} \alpha_{2} = &\frac{J_{6}}{4[\omega(q)]^2-\omega[(2q)]^2} .\end{align*}$
Here, $\omega$($q$) denotes the linear dispersion relation evaluated at wave vector $q$. The coefficients $R_1$–$R_4$ are direct measures of the lattice system. Specifically, $R_1$ is responsible for longitudinal dispersion in the $x$-direction, $R_2$ is responsible for transverse dispersion in the $y$-direction, $R_3$ captures the mixed derivative coupling and lattice anisotropy, and $R_4$ is the effective nonlinear coefficient due to the joint effect of the cubic intersite interaction $J_3$ and the on-site anharmonic potential $J_6$. Therefore, $J_2$ is the dominant coefficient for dispersive spreading, whereas $J_3$ and $J_6$ regulate nonlinear self-focusing/defocusing, soliton amplitude, and collision properties. This provides a clear link between the lattice interactions at the microscopic level and the macroscopic envelope dynamics captured by equation (Additionally, $v_{g}$ represents the group velocity, $\omega$ corresponds to linear vibrations at a specific frequency, $J_{2}$ describes the harmonic potential, $J_{3}$ accounts for the cubic potential and $J_{6}$ represents the on-site potential in this study. These equations possess a physically transparent structure: In contrast, $R_4$ characterizes the effective nonlinear response that emerges from the combined effects of the cubic intersite coupling $J_3$ and the nonlinear on-site potential $J_6$. Specifically, the cubic term drives intersite energy transfer, whereas the on-site potential restricts local trapping and amplitude. Thus, tuning $J_3$ and $J_6$ enables direct microscopic control over soliton characteristics, including width, orientation, interaction type ($X$/$Y$-pattern), and the MI bandwidth. These features, absent in standard D–S equations, emerge naturally in lattice-derived systems.
3. Derivation of solitary wave solutions in (2+1) D–S equation via Hirota bilinear formalism
$~~~$Although the D–S equation is a well-known integrable system with established soliton forms, the present analysis distinguishes itself by incorporating lattice-derived coefficients formulated within a Hamiltonian framework. The use of the Hirota bilinear method has long been known to be effective. However, there is no intent to propose any methodological originality here. The technique instead provides a viable way of extracting exact one and two-soliton solutions for a D–S-type equation derived from a lattice, where the coefficients are defined microscopically. The inverse scattering transform [29], Darboux transformation [30], Backlund transformation [31] are further options for producing soliton solutions. The technical complexity of these approaches increases when working with multidimensional systems that have anisotropic, lattice-based coefficients. Also, variational and numerical methods provide some qualitative insights into soliton behavior, but they do not offer exact forms for analyzing the redistribution of collision energy [32]. Among these, the Hirota method stands out as one of the most efficient and effective tools for generating multi-soliton solutions, as well as certain special solutions, for integrable nonlinear evolution equations. To discover the soliton collision for the obtained NLSE, we present the following rational transformation: 3 ) is derived as follows: 6 ) is solvable using power series expansions in terms of $M$ and $N$: 6 ) into equation (7 ), grouping terms with the same power of $\epsilon$, and successively solving the corresponding partial differential equations.
$\begin{align}B = \frac{M\left(x,y,t\right)}{N\left(x,y,t\right)}.\end{align}$
Let $M(x,y,t)$ denote a function that takes complex values, and let $N(x,y,t)$ represent a function that takes only real values. Through a series of symbolic transformations, the bilinear representation corresponding to equation ( $\begin{align} & \qquad\qquad [R_{2}D_{y}^{2}+R_{3}D_{xy}-R_{1}D_{x}^{2}] = 0,\nonumber\\ & \qquad [\mathrm{i}D_{t}+R_{1}D_{x}^{2}+R_{2}D_{y}^2+R_{3}D_{xy}]N^{(1)} = 0,\nonumber\\ & \qquad [R_{2}D_{y}^{2}+R_{3}D_{xy}-R_{1}D_{x}^{2}][2N^{(2)}-R_{4}M^{(1)}M^{(1)^*}] = 0,\nonumber\\ & [\mathrm{i}D_{t}+R_{1}D_{x}^{2}+R_{2}D_{y}^{2}+R_{3}D_{xy}][N^{(3)}+M^{(2)}N^{(1)}] = 0,\nonumber\\ & [R_{2}D_{y}^{2}+R_{3}D_{xy}-R_{1}D_{x}^{2}] [2N^{(4)}+2N^{(2)}-M^{(1)}M^{(3)^*}\nonumber\\ & \quad +M^{(1)^*}M^{(3)}] = 0.\end{align}$
The bilinear operators $D_x$, $D_y$ and $D_t$ are defined as follows $\begin{align*} D_x^p D_y^q D_t^r (M \cdot N) & = \left( \frac{\partial}{\partial x} - \frac{\partial}{\partial x^{^{\prime}}} \right)^p \left( \frac{\partial}{\partial y} - \frac{\partial}{\partial y^{^{\prime}}} \right)^q \left( \frac{\partial}{\partial t} - \frac{\partial}{\partial t^{^{\prime}}} \right)^r\nonumber\\ &\quad M(x, y, t) N(x^{^{\prime}}, y^{^{\prime}}, t^{^{\prime}}) \Big|_{x^{^{\prime}} = x, y^{^{\prime}} = y, t^{^{\prime}} = t}.\end{align*}$
Equation ( $\begin{align}M & = \epsilon M^{(1)}+\epsilon^{3}M^{(3)}+\epsilon^{5}M^{(5)}......,\nonumber\\N &= 1+\epsilon^{2} N^{(2)}+\epsilon^{4}N^{(4)}+\epsilon^{6}N^{(6)}.......\end{align}$
By substituting the series from equation (3.1. Single-wave soliton solution
In the process of formulating the single-soliton solution to the D–S equation, we take: 8 ) into equation (6 ) yields the explicit one-soliton solution:
$\begin{align}M &= \epsilon M^{(1)},\nonumber\\N &= 1+\epsilon^2 N^{(2)}.\end{align}$
Substituting equation ( $\begin{align}B = \frac{\mathrm{e}^{\eta_{1}}}{1+N^{\left(2\right)}},\end{align}$
where, $\begin{align}N^{(2)} & = \frac{R_{4}}{2}\frac{\mathrm{e}^{(\eta_{1}+\eta_{1}^{*})}}{(k_{1}+k_{1}^{*})^{2}},~~~\eta_{1} = k_{1}x+\omega_{1}y-\mu_{1}t-\eta_{1}(0),\nonumber\\ \mu_{1} & = \mathrm{i}{R_{1}k_{1}^2+R_{2}\omega_{1}^{2}+R_{3}k_{1}\omega_{1}}.\end{align}$
Here, $\eta_{1}$ and $\omega_{1}$ are complex wave parameters, and $k_{1}$ and $\eta_{1}(0)$ are arbitrary complex constant. The magnetization dynamics of the medium are studied by analyzing the intensity profile while varying these complex numbers and other physical parameters
3.2. Double-wave soliton solution
To derive the two-soliton solution, we assume that: 6 ) and solving the PDE gives the explicit double-soliton solution.
$\begin{align}M & = \epsilon M^{(1)}+\epsilon^{3} M^{(3)},\nonumber\\N & = 1+\epsilon^2 N^{(2)}+\epsilon^{4} N^{(4)}.\end{align}$
Substituting into equation ( $\begin{align}B = \frac{\epsilon M^{\left(1\right)}+\epsilon^{3} M^{\left(3\right)}}{1+\epsilon^2 N^{\left(2\right)}+\epsilon^{4} N^{\left(4\right)}}.\end{align}$
with, $\begin{align}M^{(1)} = &\mathrm{e}^{\eta_{1}+\eta_{2}},\nonumber\\M^{(3)} = &\mathrm{e}^{\eta_{1}+\eta_{1}^{*}+\eta_{2}+a_{1}}+\mathrm{e}^{\eta_{1}+\eta_{2}+\eta_{2}^{*}+b_{1}},\nonumber\\N^{(2)} = &\frac{R_{4}}{2(R_{2}+R_{3}-R_{1})}\Bigg[\mathrm{e}^{\eta_{1}+\eta_{1}^{*}+L_{1}}+\mathrm{e}^{\eta_{1}+\eta_{2}^{*}+L_{2}}\nonumber\\ & +\mathrm{e}^{\eta_{2}+\eta_{1}^{*}+L_{3}}+\mathrm{e}^{\eta_{2}+\eta_{2}^{*}+L_{4}}\Bigg],\end{align}$
$\begin{align*}N^{(4)} = &\mathrm{e}^{\eta_{1}+\eta_{1}^{*}+\eta_{2}+\eta_{2}^{*}+\Xi},\nonumber\\ \mathrm{e}^{a_{1}} = &\frac{R_{4}}{2(R_{2}+R_{3}-R_{1})}\Bigg[\frac{(k_{2}-k_{1})^{2}(\omega_{2}-\omega_{1})^{2}}{(k_{1}+k_{1}^{*})^{2}(k_{2}+k_{1}^{*})^{2}(\omega_{2}+\omega_{1}^{*})^{2}}\Bigg],\nonumber\\ \mathrm{e}^{b_{1}} = &\frac{R_{4}}{2(R_{2}+R_{3}-R_{1})}\Bigg[\frac{(k_{2}-k_{1})^{2}(\omega_{2}-\omega_{1})^{2}}{(k_{2}+k_{2}^{*})^{2}(k_{1}+k_{2}^{*})^{2}(\omega_{2}+\omega_{2}^{*})^{2}}\Bigg],\nonumber\\ \mathrm{e}^{L_{1}} = &\frac{1}{(\omega_{1}+\omega_{1}^{*})^{2}+(k_{1}+k_{1}^{*})(\omega_{1}+\omega_{1}^{*})-(k_{1}+k_{1}^{*})^{2}},\nonumber\\ \mathrm{e}^{L_{2}} = &\frac{1}{(\omega_{1}+\omega_{2}^{*})^{2}+(k_{1}+k_{2}^{*})(\omega_{1}+\omega_{2}^{*})-(k_{1}+k_{2}^{*})^{2}},\nonumber\\ \mathrm{e}^{L_{3}} = &\frac{1}{(\omega_{2}+\omega_{1}^{*})^{2}+(k_{2}+k_{1}^{*})(\omega_{2}+\omega_{1}^{*})-(k_{2}+k_{1}^{*})^{2}},\nonumber\\ \mathrm{e}^{L_{4}} = &\frac{1}{(\omega_{2}+\omega_{2}^{*})^{2}+(k_{2}+k_{2}^{*})(\omega_{2}+\omega_{2}^{*})-(k_{2}+k_{2}^{*})^{2}},\end{align*}$
$\begin{align} \mathrm{e}^{\Xi} &\frac{R_{4}}{4(R_{2}+R_{3}-R_{1})}\Bigg[\frac{(k_{2}^{*}-k_{1}^{*})^{2}(\omega_{2}^{*}-\omega_{1}^{*})^{2}}{(k_{2}^{*}+k_{2})^{2}(k_{1}^{*}+k_{2})^{2}(\omega_{1}^{*}+\omega_{2})^{2}}+\frac{(k_{2}^{*}-k_{1}^{*})^{2}(\omega_{2}^{*}-\omega_{1}^{*})^{2}}{(k_{1}^{*}+k_{1})^{2}(k_{2}^{*}+k_{1})^{2}(\omega_{2}^{*}+\omega_{1})^{2}}\nonumber\\ & +\frac{(k_{2}-k_{1})^{2}(\omega_{2}-\omega_{1})^{2}}{(k_{2}+k_{2}^{*})^{2}(k_{1}+k_{2}^{*})^{2}(\omega_{1}^{*}+\omega_{2})^{2}}+\frac{(k_{2}-k_{1})^{2}(\omega_{2}-\omega_{1})^{2}}{(k_{1}+k_{1}^{*})^{2}(k_{2}+k_{1}^{*})^{2}(\omega_{2}+\omega_{1}^{*})^{2}}\nonumber\\ & - \frac{(k_{2}+k_{2}^{*})^{2}+(k_{1}+k_{1}^{*})^{2}-2(k_{2}+k_{2}^{*})(k_{1}+k_{1}^{*})(\omega_{2}+\omega_{1}^{*})+(\omega_{1}+\omega_{1}^{*})^2}{(k_{1}+k_{1}^{*})^{2}(k_{2}+k_{2}^{*})^{2}(\omega_{1}+\omega_{1}^{*})^{2}(\omega_{2}+\omega_{2}^{*})^{2}}\nonumber\\ & - \frac{(k_{1}+k_{2}^{*})^{2}+(k_{2}+k_{1}^{*})^{2}-2(k_{2}+k_{2}^{*})(k_{1}+k_{1}^{*})(\omega_{1}+\omega_{2}^{*})+(\omega_{2}+\omega_{2}^{*})^2}{(k_{1}+k_{1}^{*})^{2}(k_{2}+k_{1}^{*})^{2}(\omega_{1}+\omega_{1}^{*})^{2}(\omega_{2}+\omega_{1}^{*})^{2}}\Bigg].\end{align}$
4. A study of Modulational Instability in lattice systems
Linear stability analysis is a fundamental tool widely applicable across nonlinear wave systems in nature. This analysis reveals how small-amplitude perturbations often arising from background noise can grow rapidly due to the combined effects of nonlinearity and diffraction [33]. The connection between stability and solitons becomes evident when one observes that the filaments formed through this instability process often evolve into trains of nearly ideal solitons.
According to the linear stability analysis [34], the stationary solution of the form: 15 ) into the considered equation, we get 17 ) into considered equation and linearizing, yields 18 ) into equation (19 ), we obtain following equations for $q_1$ and $q_2$, 20 ) has nonzero solutions if and only if
$\begin{align} B = b_{0} \mathrm{e}^{\mathrm{i} \left(\kappa x+\nu y-\omega t\right)},\end{align}$
where $b_0$ is the initial incidence power while $ \kappa, \nu $ and $\omega$ are the wave numbers and perturbation frequency. Substituting equation ( $\begin{equation} \omega = \kappa ^2 R_1+\nu ^2 R_2+\kappa \nu R_3- b_0^2 R_4.\end{equation}$
Assume that the perturbed solution of the form: $\begin{equation} B = \left[b_0+\theta \phi \left(x,y,t\right)\right]~\mathrm{e}^{\mathrm{i}\left( \kappa x+ \nu y-\omega t\right)},\end{equation}$
where $\theta \lt \lt 1$ is a perturbation parameter and $\phi $ can be expressed as $\begin{equation} \phi \left(x,y,t\right) = q_1~\mathrm{e}^{\mathrm{i}\left( \hat{\kappa} x+\hat{\nu } y-\hat{\omega } t\right)}+q_2~\mathrm{e}^{-\mathrm{i}\left( \hat{\kappa} x+\hat{\nu } y-\hat{\omega } t\right)},\end{equation}$
where $q_1$ and $q_2$ are the constants. Inserting equation ( $\begin{align} & R_4 b_0^2 \left(\phi +\phi ^*\right)+\mathrm{i} \phi _t+\mathrm{i} \left(2\nu R_2+\kappa R_3 \right)\phi _y\nonumber\\ & \quad +R_2 \phi _{\text{yy}}+\mathrm{i}\left( 2\kappa R_1+\nu R_3 \right)\phi _x+R_3 \phi _{\text{xy}}\nonumber\\ & \quad +R_1 \phi _{\text{xx}} = 0.\end{align}$
Substituting equation ($\begin{align*} & \psi _{11} q_1+\psi _{12}q_2 = 0,\end{align*}$
$\begin{align} & \psi _{21} q_1+\psi _{22} q_2 = 0,\end{align}$
where $\begin{align*} \psi _{11} &= R_4 b_0^2, \nonumber\\ \psi _{12} &= -\hat{\omega }+\left(2 \kappa -\hat{\kappa }\right) \hat{\kappa } R_1+\left(2 \nu -\hat{\nu }\right) \nonumber\\ & \quad \times \hat{\nu } R_2+\left(\hat{\kappa } \left(\nu -\hat{\nu }\right)+\kappa \hat{\nu }\right) R_3+R_4 b_0^2, \nonumber\\ \psi _{21} &= \hat{\omega }-\hat{\kappa } \left(2 \kappa +\hat{\kappa }\right) R_1-\hat{\nu } \left(2 \nu +\hat{\nu }\right)\nonumber\\ & \quad R_2-\left(\kappa \hat{\nu }+\hat{\kappa } \left(\nu +\hat{\nu }\right)\right) R_3+R_4 b_0^2, \nonumber\\ \psi _{22}& = R_4 b_0^2 .\end{align*}$
Equation ( $\begin{equation} \begin{vmatrix} \psi _{11} & \psi _{12} \\ \psi _{21} & \psi _{22} \end{vmatrix} = 0.\end{equation}$
Solving the above determinant, the following dispersion relation can be obtained $\begin{align} \hat{\omega } = 2 \kappa \hat{\kappa } R_1+2 \nu \hat{\nu } R_2+\left(\nu \hat{\kappa }+\kappa \hat{\nu }\right) R_3\pm \sqrt{\triangle },\end{align}$
where $\Delta = (\hat{\kappa }^2 R_1+\hat{\nu }^2 R_2+\hat{\kappa } \hat{\nu } R_3) (\hat{\kappa }^2 R_1+\hat{\nu }^2 R_2+\hat{\kappa } \hat{\nu } R_3$ $-2 R_4 b_0^2)$.If $\Delta \unicode{x2A7E} 0$, then $\hat{\omega }$ is real and hence the state for the considered equation is stable. Contrarily, if $\Delta \lt 0$, then $\hat{\omega }$ is imaginary and the state becomes unstable. In contrast to classical D–S systems, the modulation instability spectrum in this model is fully tunable through the microscopic lattice potentials. The deformation of the MI surface, from an isotropic bowl to a directional instability tongue, directly results from the lattice-induced anisotropy encoded in $R_1$, $R_2$, and $R_3$. The modulation instability gain spectrum $G(\kappa)$ can be derived as follows:
$\begin{align} G\left(\kappa\right)& = 2 {\text{Im}}\left(\hat{\omega }\right) \nonumber\\ &= 2 {\text{Im}}\left[2 \kappa \hat{\kappa } R_1+2 \nu \hat{\nu } R_2+\left(\nu \hat{\kappa }+\kappa \hat{\nu }\right) R_3\pm \sqrt{\triangle }\right].\end{align}$
The spectrum of MI gain is influenced by the values for the anisotropic dispersion coefficient $R_1$, $R_2$, and the mixed term $R_3$, which together shape the MI gain spectrum. The $R_1$ and $R_2$ provide control of the instabilities occurring along the principal lattice axes, and $R_3$ couples the longitudinal and transverse perturbations, thereby rotating and distorting the MI gain surface in the ($\kappa$, $\nu$)-plane. Thus, this procedure anisotropically selects MI instability tongues and directionally selects MI determined by the lattice anisotropy.5. Results and discussion
We used two different analytical approaches to investigate the nonlinear dynamic behavior of the (2+1) dimensional D–S system that arise from a lattice. The Hirota bilinear method produces exact localized wave solutions while MI can provide information about the stability of extended wave backgrounds and which parameter ranges that are most conducive for soliton formation. It is found that the two approaches can be integrated to obtain a comprehensive understanding of nonlinear energy localization within multi-dimensional lattices.
5.1. Soliton dynamics and collision characteristics
It is noteworthy that the double soliton solution is defined by arbitrary complex parameters $k_{j}$ and initial phases $\eta_{j}(0)$, $j$ = $1$, $2$, $3$ and it displays variations in amplitude intensity across individual solitons. This study analyzes soliton collision dynamics in multidimensional lattices, emphasizing the impact of the cubic potential $J_{3}$ and the on-site potential $J_6$ on lattice behavior. Evolution plots based on equations (9 ) and (12 ) illustrate key aspects of nonlinear energy transport. These visuals illustrate soliton propagation and interaction in nonlinear optical lattices, where unique behaviors arise from the interplay of harmonic coupling and external potentials. The soliton solutions from equations (9 ) and (12 ) give a clear idea of how they behave, as shown in figures 1–4. Since the solution, equation (9 ) produces soliton profiles similar to earlier ones. Figures 1(a)–(c) show how the soliton behaves in three different planes based on the parameter values $J_{3}$ = $0.1$, $J_{6}$ = $0.1$, $J_{2}$ = $1$, $q_{2}$ = $0.1$, $q_{1}$ = $0.4$, $\lambda_{1}$ = $0.55$, $\lambda_{2}$ = $0.06$, $n$ = $0.1$, $V_{g}$ = $0.6$, $m$ = $0.1$, $\omega_{1}$ = $0.15$, $\omega_{1}^{*}$ = $0.15$, $k_{1}$ = $1.2$, $k_{1}^{*}$ = $1.2$. In the $x-t$ plane (figure 1(a)), the soliton’s core features are clear. In the $x-y$ plane (figure 1(b)), it shows periodic amplitude and orientation changes, while in the $y-t$ plane (figure 1(c)), it narrows without altering its peak. Figure 2 shows that increasing $J_3$ rotates and compresses the soliton, whereas varying $J_6$ mainly shifts it rightward. These changes highlight soliton sensitivity to parameter variations, crucial in optical applications. In fiber communications and photonic crystals, solitons maintain pulse shape over long distances, while external potentials enable control for improved signal transmission, nonlinear waveguides, and all-optical switching.
Figure 1. Propagation of the one-soliton solution in the (a) $x-t$ plane at $y$=$0.1$, (b) $x-y$ plane at $t$=$0.1$, and (c) $y-t$ plane at $x$=$0.1$, illustrating its spatial-temporal evolution. |
Figure 2. One-soliton profiles in the $y-x$ plane for varying cubic ($J_3$) and on-site ($J_6$) potentials, showing width narrowing, orientation change, and positional shifts. |
Figure 3. $X$-type two-soliton collision profiles for increasing cubic potential $J_3$ = $0.1$, $0.3$, and $0.6$, with 3D intensity surfaces, peak-amplitude evolution, and corresponding 2D distributions. |
Figure 4. $Y$-type soliton interactions for increasing on-site potential $J_6$ = $0.1$, $0.3$, and $0.6$, showing amplitude suppression, positional shifts, and corresponding 2D distributions. |
The collision behaviors discussed, specifically the controllable transition between $X$-type and $Y$-type interactions, originate solely from the nonlinear coefficients induced by the lattice. Classical D–S equations with constant phenomenological parameters cannot reproduce these parameter-controlled collision transitions. We next examine the dynamics of two-soliton structures under cubic and on-site nonlinearities. To obtain the two-soliton solutions of equation (12 ), we approximate equation (3 ) and derive the field $B(x,y,t)$. The intensity variations observed during soliton collisions stem from the combined effects of the cubic interaction term $J_{3}$, the on-site potential $J_{6}$, and the intrinsic system parameters.
We first examine how the cubic potential $J_{3}$ governs the pairwise interactions, primarily controlling the energy exchange between solitons. This results in energy-sharing collisions characterized by noticeable intensity redistribution and amplitude-dependent position shifts, making the interaction inherently inelastic. This behavior is illustrated in figures 3(a)–(c) for the parameter values, $J_{3}$ = $0.1$, $0.3$ and $0.6$, $J_{6}$ = $0.5$, $m$ = $0.1$, $q_{1}$ = $0.1$, $q_{2}$ = $0.4$, $\lambda_{1}$ = $0.5$, $\lambda_{2}$ = $1$, $n$ = $5$, $V_{g}$ = $1$, $k_{1}$ = $9+1i$, $k_{1}^{*}$ = $4.1$, $k_{2}$ = $2.5-1i$, $k_{2}^{*}$ = $2+1i$, $\omega_{1}$ = $1.1+1i$, $\omega_{1}^{*}$ = $0.8-1i$, $\omega_{2}$ = $5.1$, $\omega_{2}^{*}$ = $4.4$, $\eta_{1}(0)$ = $7-1i$, $\eta_{1}(0)^{*}$ = $2.5+1i$, $\eta_{2}(0)$ = $-2.1+1i$ and $\eta_{2}(0)^{*}$ = $2+1i$. The contour and 2D profiles in figure 3 illustrate how the soliton pair redistributes energy during collision, showing a clear exchange and localization of energy as the parameter $J_{3}$ increases. The emergence of collisional solitons adds richer dynamical behavior beyond standard interactions, with $J_{3}$ playing a central role in shaping their collision patterns and energy flow. By tuning $J_{3}$, the strength of nonlinear coupling and the degree of soliton trapping can be precisely controlled. In multidimensional optical lattices, this cubic nonlinear potential stabilizes localized beams, governs their transverse spreading, and enables fine control of $X$-type soliton steering, confinement, and interaction strength features relevant for all optical switching and energy routing applications in nonlinear photonic lattices.
Secondly, the $Y$-shaped collision dynamics in figures 4(a)–(c) show clear intensity redistribution and amplitude-dependent position shifts between interacting solitons. As the on-site potential parameter $J_6$ increases from $0.1$ to $0.3$ and $0.6$, the overall soliton amplitudes gradually diminish, giving rise to pronounced $Y$-type interaction patterns. For $J_6$, one soliton’s amplitude grows from about $7$ to $19$ units, while the other increases from roughly $13$ to $18$ units, reflecting strong energy exchange during the interaction. The analysis is conducted using parameter values, $J_{6}$ = $0.1$, $0.3$ and $0.6$, $J_{3}$ = $0.2$, $J_{2}$ = $0.05$, $m$ = $0.1$, $q_{1}$ = $0.1$, $q_{2}$ = $0.4$, $\lambda_{1}$ = $0.5$, $\lambda_{2}$ = $1.1$, $n$ = $5$, $V_{g}$ = $1.5$, $k_{1}$ = $8.1+5i$, $k_{1}^{*}$ = $2-1.5i$, $k_{2}$ = $4.3-1i$, $k_{2}^{*}$ = $2+1.1i$, $\omega_{1}$ = $2.45+1i$, $\omega_{1}^{*}$ = $0.5-1i$, $\omega_{2}$ = $3$, $\omega_{2}^{*}$ = $1.1$, $\eta_{1}(0)$ = $4.5+1i$, $\eta_{1}(0)^{*}$ = $9.5-1i$, $\eta_{2}(0)$ = $-1.9+1i$, $\eta_{2}(0)^{*}$ = $4.5+1i$. The interactions exhibit clear shape-changing dynamics, most prominently inelastic $Y$-shaped solitary-wave collisions. As the on-site potential parameter $J_6$ increases, the soliton amplitudes gradually diminish and their positions shift. At $J_6$ = $0.3$, the solitons stay well apart, whereas increasing $J_6$ to $0.6$ brings them noticeably closer, as shown in the 2D profiles.
These snapshots, taken before and after the collision, show how stronger on-site potentials markedly alter both the spacing and amplitude of the interacting waves. Such tunable lattice–soliton behavior is crucial for applications ranging from magnetic data storage to optical communication systems. Specifically, $Y$-shaped solitons in multidimensional lattices with tunable on-site potentials enable multidirectional signal routing and controlled energy redistribution in photonic networks. By manipulating soliton shapes and trajectories through on-site potential engineering, such control allows the implementation of robust soliton-based logic gates and multidimensional information processing. Table 1 summarizes the novelty of our work, highlighting how embedding the (2+1)-dimensional D–S model in a lattice uncovers new soliton dynamics and how the combined cubic and on-site potentials produce unique $X$- and $Y$-type collisions. The transition from $X$-type to $Y$-type interactions reflects a competition between energy transfer (controlled by $J_3$) and local confinement (controlled by $J_6$).
Table 1. Comparison of parameter effects between previous studies and current findings. |
| Figure | What the plot shows (this work) | Comparison with earlier literature | Physical / technical impact of parameters |
|---|---|---|---|
| Figure 1 | Single-soliton evolution in $x$–$t$, $x$–$y$, and $y$–$t$ planes showing periodic amplitude changes, orientation shifts, and $y$-direction narrowing for baseline values ($J_3 = J_6 = 0.1$). | Similar propagation patterns reported in DS/NLS studies, but earlier papers use continuum parameters. Present work emphasizes lattice-derived $J$-coefficients. | Small $J_3$ and $J_6$ preserve amplitude and shape, indicating weak-lattice limit behaves like continuum DS. Serves as baseline for collision analysis. |
| | |||
| Figure 2 | Single-soliton profile variations in the $y$–$x$ plane. Increasing $J_3$ narrows and tilts the soliton; increasing $J_6$ shifts it rightward. | Previous gap-soliton and DS studies show width/orientation changes via dispersion–nonlinearity interplay. Here changes are explicitly linked to $J_3$ and $J_6$. | $J_3 \uparrow$: narrower, tilted soliton. $J_6 \uparrow$: positional drift due to on-site potential. Enables experimental steering ($J_6$) and confinement ($J_3$). |
| | |||
| Figure 3 | Bright–dark soliton interaction for varying $J_3$ ($0.1 \rightarrow 0.6$). Higher $J_3$ produces inelastic collisions with energy transfer and position shifts. | Earlier Hirota/DS works show elastic or mildly inelastic collisions. Here, strong inelasticity arises solely from lattice-induced $J_3$. | $J_3 \uparrow$: stronger energy exchange, amplitude/trajectory shift. Cubic inter-site coupling enables energy routing in lattices. |
| | |||
| Figure 4 | $Y$-shaped soliton interactions vs $J_6$ ($0.1 \rightarrow 0.6$). Solitons move closer and amplitudes decrease, producing $Y$-type fission/fusion. | $Y/X$ patterns reported in multi-component DS models, usually requiring phase tuning. Present results show on-site potential $J_6$ alone produces $Y$-structure. | $J_6 \uparrow$: amplitude suppression, stronger interactions, $Y$-type merging/fission. On-site potential acts as trap/damping for soliton morphology control. |
When the cubic intersite coupling dominates, solitons efficiently exchange energy, leading to amplitude-dependent trajectory shifts and $X$-type branching. When the on-site potential becomes comparable or stronger, solitons experience localized trapping that suppresses amplitude growth and forces a merging–splitting behavior characteristic of $Y$-type interactions. These transitions have no analogue in continuum D–S equations, demonstrating that lattice potentials create new regimes of multidimensional soliton physics. Through experimental control of lattice geometries and trapping potentials, it is possible to control the cubic intersite and on-site nonlinearities, thus offering a route to test the predicted $X$ and $Y$-type soliton interactions and control energy localization in real systems.
This is an important time to explain clearly the purpose of using the analytic method to derive two-soliton solutions to the D–S equation. The Hirota bilinear method provides explicit two-soliton solutions of a (2+1)D D–S equation and closed-form solutions for multidimensional parameters with lattice properties.
This means that by having an explicit closed-form representation for the D–S equation in terms of the microscopic parameters $J_{3}$ and $J_{6}$ the exact two-soliton solutions can directly provide collision geometries ($X$ and $Y$-type). These solutions also show how much the amplitude redistributes and how the solitons deviate from their original paths of travel through the lattice system. In contrast, different analytic methods in solving (2+1) systems with anisotropic and lattice-dependent coefficients become significantly more complicated to use for multi-dimensional lattice systems. For these reasons, the Hirota bilinear method provides the best opportunity to determine the dynamics of controlled soliton interactions utilized throughout this paper.
5.2. Modulational Instability and stability landscape
To investigate the MI of the multidimensional lattice system, we evaluated equation (23 ) for various parameter sets, as shown in figures 5(a)–(c). These figures illustrate how the MI growth rate $G$($\kappa$, $\nu$) varies with the longitudinal ($\kappa$) and transverse ($\nu$) perturbation wavenumbers. The results indicate that energy localization, the cubic potential, and the on-site potentials collectively support the formation of pulse-like solitary excitations. The optical solitons supported by these lattices exhibit well-defined stable and unstable regions. Initially, tuning the cubic potential parameter $J_3$ reshapes the MI surface, revealing the impact of nonlinear coupling on multidimensional instability.
Figure 5. MI modulational instability surfaces and gain cross-sections for $J_3$ = $0.1$, $0.3$, and $0.6$, illustrating the transition from a broad instability peak to a tilted, narrowed, anisotropic band. |
In figure 5(a) (weak cubic regime), the MI surface shows a broad but shallow peak with an almost symmetric instability region, indicating moderate perturbation growth across a wide wavenumber range. This corresponds to the parameter values, $J_3$ = $0.1$, $0.3$, $0.6$, $J_6$ = $1.1$, $J_2$ = $0.1$, $m$ = $15.1$, $q_1$ = $0.1$, $q_2$ = $0.4$, $\lambda_1$ = $0.1$, $\lambda_2$ = $0.1$, $n$ = $0.5$, and $V_g$ = $1.5$. This behavior indicates weak nonlinearity, where the cubic term induces only minor deviations from the quadratic lattice MI profile. The 2D plots show that the gain starts from negative values and rises gradually to a modest peak near $\kappa \approx 0$. For $J_3$ = $0.3$ (figure 5(b)), the MI peak sharpens and becomes more localized, reflecting a higher growth rate. The asymmetric surface suggests that the cubic potential enhances nonlinear energy transfer along specific directions, producing anisotropic instability. Overall, a stronger cubic term amplifies modulation growth and shifts the most unstable mode, yielding a steeper, more pronounced instability profile.
When $J_3$ = 0.6 (figure 5(c)), the MI surface exhibits a noticeable tilt, and the instability peak shifts toward a specific region in the ($\kappa$, $\nu$) plane. This behavior indicates that stronger cubic interactions narrow the unstable band while enhancing the maximum growth rate, resulting in highly localized and directionally biased instability. The corresponding 2D MI curves become strongly asymmetric, displaying a steep rise and an extended unstable region near positive $\kappa$. Overall, increasing the cubic interaction intensifies MI by
| • Raising the peak growth rate, | |
| • Narrowing and shifting the instability band, | |
| • Introducing anisotropic mode growth, | |
| • Selectively stabilizing/destabilizing modes. |
Consequently, the MI spectrum becomes highly sensitive to the lattice’s nonlinear structure, significantly affecting energy localization Secondly, fine-tuning the on-site potential $J_6$ produces a smooth, bowl-shaped 3D MI gain surface $G$($\kappa$, $\nu$) with gentle curvature, indicating only mild variations in instability across the ($\kappa$, $\nu$) plane. The analysis is carried out under the parameter sets, $J_6$ = $0.1$, $0.3$, $0.6$, $J_3$ = $1.2$, $J_2$ = $0.1$, $m$ = $4.1$, $q_1$ = $0.1$, $q_2$ = $0.4$, $\lambda_1$ = $0.5$, $\lambda_2$ = $1$, $n$ = $5$ and $V_g$ = $1$. The minimum gain occurs near $\kappa$ = $0$, suggesting that long-wavelength perturbations experience only weak growth. The corresponding 2D gain curve $G$($\kappa$) exhibits a symmetric $U$-shaped profile, confirming that MI persists over a broad range of wavenumbers but with relatively low growth rates as depicted in figure 6(a). Because a small $J_6$ provides only a weak restoring force, MI can emerge, but its overall strength remains limited. As the on-site potential increases to $J_6$ = $0.3$, the MI gain surface becomes more irregular and asymmetric, as illustrated in figure 6(b). The central valley deepens, and sharper gradients emerge along the $\kappa$-direction, indicating that certain wavenumber combinations amplify much faster than others. The 2D cut clearly reflects this, showing a rapid growth for negative $\kappa$ and much weaker amplification for positive $\kappa$. A moderate $J_6$, therefore, enhances nonlinear coupling and produces a direction-dependent MI response typical of multidimensional lattices. In this regime, the system becomes more unstable overall, but the instability spectrum loses its symmetry.
Figure 6. MI surfaces and gain curves for $J_6$ = $0.1$, $0.3$, and $0.6$, showing increasing asymmetry, central suppression, and confinement of instability to higher wavenumbers. |
In the strong on-site potential regime, the MI gain surface develops a broad, flat region around $\kappa$ = $0$, indicating complete suppression of MI at small wavenumbers, as shown in figure 6(c). The surface rises only near the edges of the ($\kappa$, $\nu$) domain, meaning that instability survives only for sufficiently large $\kappa$. The corresponding 2D gain curve $G$($\kappa$) reinforces this trend, displaying a zero-gain plateau in the central region and nonzero growth only near the boundaries. These features show that a strong on-site potential stabilizes the system by suppressing low-frequency perturbations and confining MI to high-frequency modes. Collectively, the three figure sets clearly demonstrate that the on-site potential $J_6$ is a key stabilizing factor, directly shaping the strength, symmetry, and bandwidth of MI. We have also demonstrated the 3D MI growth-rate surfaces in figures 7(a), (b) and 8(a), which show that increasing the cubic potential progressively deforms the instability landscape in the ($\kappa$, $\nu$) plane. For small cubic potential values, the surface forms a closed, symmetric bowl centered at the origin, matching the elliptical instability region in figure 7(a). In this weakly nonlinear regime, MI is essentially isotropic, allowing a broad range of ($\kappa$, $\nu$) perturbations to grow. As the cubic potential increases, the 3D surface stretches and tilts, opening the unstable region primarily along $\kappa \gt 0$, as seen in figure 7(b). This deformation indicates the emergence of anisotropic MI, with instability concentrated in a narrower band of modes. Thus, adjusting the cubic potential reshapes the MI surface from a symmetric pocket to a directional instability tongue, showing how nonlinear lattice potentials control MI strength and directionality. When the on-site potential is varied, the 3D MI growth-rate surface changes markedly, revealing how the local restoring force governs the stability of long-wavelength perturbations. For a weak on-site potential (figure 8(a)), the surface becomes broad and relatively flat, creating a large, nearly uniform unstable region in the ($\kappa$, $\nu$) plane. As the on-site potential increases (figure 8(b)), the 3D surface contracts and develops a sharp central peak. The instability is then confined to a small neighbourhood around ($\kappa$, $\nu$) = ($0$, $0$), with the outer region becoming rapidly stable. Physically, increasing the on-site potential strengthens local binding, suppressing transverse distortions and limiting MI to perturbations in the near-zero-wavenumber regime. Thus, adjusting the on-site potential narrows and sharpens the MI spectrum, reducing the instability window and stabilizing most high-wavenumber modes. These results show that tuning the cubic nonlinearity and on-site lattice potential provides precise control over MI bandwidth, gain, and instability onset. Such tunability enables stable soliton channels, suppresses beam breakup, and supports robust spatiotemporal confinement. The nonlinear lattice potentials offer a powerful route for stabilizing and shaping multidimensional solitons in modern photonic systems.
Figure 7. 3D MI growth-rate snapshots for weak and strong cubic potentials ($J_3$ = $0.1$ and $J_3$ = $0.16$), highlighting the deformation from an isotropic pocket to an anisotropic instability tongue. |
Figure 8. 3D MI growth-rate snapshots for weak and strong on-site potentials ($J_6$ = $0.1$ and $0.6$), showing the transition from a broad unstable region to a sharp, localized instability peak. |
6. Conclusions
We have established a (2+1)-dimensional D–S-type equation derived from the lattice, whose coefficients reflect the cubic intersite and nonlinear on-site potentials of the underlying Hamiltonian. This microscopic derivation reveals aspects of the system’s physics that remain hidden in conventional D–S models with phenomenological parameters.
| • The cubic intersite potential $J_3$ serves as a nonlinear energy-transfer channel, regulating soliton width, orientation, and the inelasticity of two-soliton interactions. Higher $J_3$ values enhance amplitude redistribution, yielding $X$-type collision patterns with tunable energy exchange, an effect not previously reported in DS studies. | |
| • The on-site nonlinear potential $J_6$ introduces localized trapping and amplitude suppression, promoting the emergence of $Y$-type collision geometries. This mechanism offers a new interpretation of shape-changing collisions based on microscopic lattice physics rather than external phase engineering. | |
| • MI analysis shows that $J_3$ and $J_6$ distinctly shape the instability landscape: $J_3$ enhances directional gain, forming anisotropic tongues, whereas $J_6$ suppresses long-wavelength perturbations, narrowing the unstable band. This highlights the direct role of lattice potential strengths in governing multidimensional soliton stability. |
Overall, the study demonstrates that microscopic nonlinearities in discrete lattices provide powerful control over multidimensional soliton propagation, collision morphology, and MI. The lattice-derived D–S framework enables tunable nonlinear waveguides, photonic lattices, and energy-transport systems, with soliton behavior governed by intrinsic lattice interactions.