Nonlinear partial differential equations (NPDEs) are fundamental mathematical models for describing complex phenomena in various branches of science, including fluid dynamics, plasma physics, optical communications, and quantum mechanics [1]. Optical solitons are self-reinforcing solitary wave packets that maintain their shape and velocity during propagation due to a precise balance between nonlinear self-phase modulation and linear group velocity dispersion [2]. In optical fiber systems, bright solitons form when anomalous dispersion compensates Kerr nonlinearity, while dark solitons intensity dips on a continuous wave background-arise in the normal dispersion regime.

Bifurcation theory is the mathematical study of qualitative changes in the structure of the phase portrait of a dynamical system as one or more parameters vary continuously [3]. In the context of nonlinear wave equations, bifurcation analysis identifies critical parameter values at which the number or stability of equilibrium points changes, revealing qualitative transitions such as the emergence of periodic orbits or the merging of equilibria. The natural nonlinear and dispersive nature of these models enable the generation of localized wave structures, called solitons, that keep their shape and velocity during propagation and collision processes [4]. Recently, much research has been done, establishing powerful analytical and numerical techniques to discover exact solutions to NPDEs such as the inverse scattering transform [5], the bilinear method proposed by Hirota [6], Bácklund transformations [7], Darboux transformations [8], the extended tanh-function method [9], the Exp-function method [10], $(G^{^{\prime}}/G)$-expansion method [11], modified Kudryashov method [12] and Sardar sub-equation method [13]. The nonlinear Schrödinger equation (NLSE) and its variants are among the various classes of NPDEs, and they have received a lot of interest lately, on account of their central relevance for describing the evolution of a wave packet in a dispersive medium [14]. The NLSE for pulse propagation in optical fibers is determined by dispersion, where balance between dispersion and nonlinearity cause soliton formation [15]. Numerous general forms of the Schrödinger equation have been investigated, for example, cubic-quintic NLSE [16], higher-order NLSE with fourth-order dispersion [17], derivative NLSE [18], coupled NLSE systems [19], as well as space-time fractional NLSE [20], amongst others. Soliton solutions have been generated by employing various sophisticated methods from these models, such as the modified simple equation method [21], the new extended direct algebraic method [22], the unified Riccati equation expansion method [23], and the generalized exponential rational function method [24]. In the case of multi-dimensional wave propagation phenomena, integrable systems with (2+1)-dimensional structures have been of considerable importance in modeling transverse perturbations and instabilities [25]. Originally intended to characterize the nonlinear spin dynamics for ferromagnetics, the Myrzakulov family of equations has been used in both nonlinear optics as well as in plasma physics [26]. The previous studies on Myrzakulov-type equations and their solution structure have greatly improved. Kumar et al [27] studied the Myrzakulov-I equation using the modified exp-function method and identified numerous soliton families. Rezazadeh et al [28] proposed the new Kudryashov method on the Myrzakulov-II equation, and obtained hyperbolic and trigonometric function solutions. Ali and Mehanna [29] successfully used the unified method to extract optical solitons from the Myrzakulov-III equation resulting in bright, dark, and singular soliton propagation. As a newly proposed integrable (2+1)-dimensional model we extend the classical Myrzakulov framework through Akbota-Myrzakulov-Tolkynay-Zhaidary (AMTZ) equation. The AMTZ system can be illustrated with the following coupled NPDEs [30, 31]: where $\Omega(x,y,t)$ is the primary wave amplitude function, $\Pi(x,y,t)$ is the auxiliary wave function and $X(x,y,t)$ is a potential function. Our independent variables are $x$ and $y$ (space coordinates) and $t$ (temporal variable). The parameter $\beta$ is the nonlinearity coefficient related to the amplitude of self-phase modulation and cross-phase modulation effects, $b$ is a characteristic scaling factor associated with spatial dimensions and properties of materials, and $a$ is considered the intensity of higher-order dispersive effects. In equation (1), the temporal derivative $\Omega_t$ expresses the temporal evolution, the advection parameter $-X\Omega_x/b$ describes the transport influenced by the potential field, the nonlinear parameter $\beta\Omega\Omega_y/b$ indicates the transverse modulation, and $\beta\Pi_y$ denotes the coupling between fields. In equation (2), the higher-order dispersion term $\beta\Omega_{xxy}/(4ab)$ creates mixed derivative effects combining longitudinal and transverse differences. The constraint equation (3) establishes an integral connection between $X$ and $\Omega$. In spite of its mathematical framework and practical usage in electromagnetic wave propagation and magneto-acoustic wave dynamics [32], the AMTZ equation has never been much studied. Mathanaranjan and Myrzakulov [31] established its Lax pair representation to validate integrability and did not present explicit analytical solutions. Myrzakulov et al [30] have presented the AMTZ system in integrable spin models but did not examine solutions comprehensively. More advanced solution approaches have shown promise in solving similar models recently. Bin et al [33] used modified extended tanh expansion for coupled dispersive systems. Alquran [34] applied the sine-Gordon expansion method for (2+1)-dimensional equations. Raza et al [35] used unified Riccati equation technique for intricate Ginzburg-Landau equations. Unfortunately, the previous studies concentrated mainly on solution generation rather than dynamical system structure. Kamel et al have proposed an improved modified Sardar sub-equation method (IMSSEM) [36], which allows us to take on greater flexibility as the auxiliary equation $G^{^{\prime} 2}(\xi) = C G^4(\xi) + B G^2(\xi) + A$ allows us to derive larger solution families. This provides a valuable insight for assessing the stability of a network [37]. Bifurcation analysis investigates how equilibrium structures shift with parameter variations [38]. Chu et al [39] studied bifurcation phenomenon in Boussinesq equation. Jhangeer et al [40] studied chaos of perturbed NLSEs through the use of Lyapunov exponents and Poincaré sections. Nevertheless, despite these developments, no study has studied the dynamic behavior of bifurcation structures and chaotic dynamics in the AMTZ equation. The motivation to investigate the AMTZ equation derives from its application in multi-mode fiber communications, pulse dynamics in the context of dispersive optical systems, as well as wave propagation phenomena [39, 41-45]. However, a number of pressing research gaps persist: no systematic classification of exact analytical solutions has been established for the AMTZ system; relationship between model parameters and solution properties has not been established; phase plane analysis and equilibrium systems have not yet been investigated; the system's external periodic forcing response has not yet been considered; and it remains unstructured to compare the presented approaches with the existing solution methods demonstrating that the IMSSEM provides the advantages. The challenges and limitations identified in previous research could be addressed through a systematic method of constructing exact soliton solutions using the improved modified simple equation method, bifurcation analysis through the use of equilibrium point classification and phase portrait building, chaotic dynamics studies regarding periodic perturbations, validated by Lyapunov exponents, and solutions stability characterization through the analysis of temporal evolution, and physical implications for applications in nonlinear optics and plasma physics. The novelty of this work is threefold. First, this is the first study to apply the IMSSEM to the AMTZ equation, generating seventeen distinct exact solution families (hyperbolic, trigonometric, and rational) in a single unified framework a scope that surpasses all prior analytical treatments of Myrzakulov-type equations (see table 2). Second, this work presents the first systematic bifurcation analysis and phase portrait construction for the AMTZ system, classifying all equilibrium points and revealing critical parameter thresholds governing qualitative transitions. Third, this is the first study to demonstrate chaotic dynamics in the AMTZ equation via positive Lyapunov exponents and fractal Poincaré sections under external periodic perturbations, establishing that this integrable system possesses rich dynamical complexity when disturbed.

The IMSSEM offers three key advantages over standard methods such as the extended tanh method [46], the $(G^{^{\prime}}/G)$-expansion method [47], and the modified Kudryashov method [48]. (i) Solution diversity: the auxiliary equation $G^{^{\prime} 2} = CG^4 + BG^2 + A$ encodes three independent parameters $(A, B, C)$, allowing simultaneous derivation of rational, hyperbolic, and trigonometric families in a single algebraic system, whereas methods such as $(G^{^{\prime}}/G)$-expansion produce only one family type per run. (ii) Generality: the IMSSEM subsumes several sub-equation methods as special cases (e.g. setting $A = 0$ recovers Sardar sub-equation results), making it strictly more general. (iii) Computational efficiency: the polynomial structure of the method reduces the solution search to a single algebraic system solvable by Maple, avoiding the iterative steps required in numerical or transformation-based methods. The AMTZ scheme can thus accommodate a variety of analytical families via IMSSEM, showing computational efficiency gains compared to traditional approaches. Bifurcation analysis discloses critical transitions of individual parameters and chaos analysis demonstrates that the integrable system is sensitive to external perturbations that can be crucial in applications with real-world effects.

The remainder of the paper is organized as follows. Section 2 introduces the traveling wave transformation that reduces the AMTZ system to an ODE. Section 3 applies the IMSSEM to derive explicit solution families. Section 4 performs bifurcation analysis. Section 5 investigates chaotic dynamics. Section 6 discusses results and physical implications. Section 7 presents conclusions, and section 8 outlines key findings, limitations, and future directions.

The traveling wave transformation method must be used in order to solve equation (1). Here, $U(\xi)$, $V(\xi)$, and $R(\xi)$ are real-valued profile functions of the traveling wave coordinate $\xi = px + qy - vt$, where $p$ and $q$ are the wave numbers in the $x$- and $y$-directions, respectively, and $v$ is the wave propagation speed. The functions $\Omega$, $\Pi$, and $X$ depend only on $\xi$ under this reduction, converting the PDE system into a system of ODEs.

The following are obtained by substituting the transformations from equation (4) into equations (1)-(3): Integrating equation (7) with respect to $\xi$ and applying the localized wave boundary conditions $U \to 0$ and $V \to 0$ as $|\xi| \to \infty$ (i.e. requiring solutions to vanish at infinity) , the integration constant is set to zero [49-54], yielding substitute equation (8) into equation (5), we obtain Integrating equation (9) and setting the integration constant to zero under the same localized boundary conditions, we obtain Substitute equations (8) and (10) into equation (6), we obtain Integrating equation (11) and setting the integration constant to zero under the same localized boundary conditions, we obtain We find that the balancing parameter is $N = 1$ by homogeneous balance analysis of equation (12), taking into account the link between the nonlinear term $V^{3}$ and the highest-order derivative term $V_{\xi\xi}$.

In this section, we apply the IMSSEM to construct novel optical soliton solutions of the governing model. We assume that the solution of equation (12) can be expressed in the finite series form where $\delta_{i}$ $(i = 0,1,2,\ldots,N)$ are unknown constants to be determined and $N$ is a positive integer obtained via the homogeneous balance principle. By balancing the highest-order derivative term with the nonlinear term in equation (12), we obtain $N = 1$. Consequently, equation (13) reduces to where the auxiliary function $G(\xi)$ satisfies the nonlinear ordinary differential equation with $A$, $B$, and $C$ being real constants.

Applying the Galilean transformation for equation (12), we obtain the following dynamical system: To find the equilibrium points of system (51), we set $V_\xi = 0$ and $\Phi_\xi = 0$. From the first equation, $\Phi = 0$. Substituting into the second equation gives $4ab^2v^2 V + 12abpv V^2 + 8ap^2 V^3 = 0$, which factors as $4aV(bv + 2pV)(bv + pV) = 0$. This yields three solutions: $V = 0$, $V = -bv/(2p)$, and $V = -bv/p$, giving the equilibrium points where $S_1 = -\frac{b v}{p}$, $S_2 = -\frac{b v}{2 p}$.

Moreover, this is also the first integral of the system equation (51) which we have used along with the points for phase portrait

Therefore we have the following property:

Then there are the following different dynamical system cases:

Chaotic behavior in nonlinear dynamical systems is a fundamental phenomenon showing extreme sensitivity to initial conditions [40]. The perturbed dynamical system studied is equation (55), where the unperturbed Hamiltonian system (51) is forced by the periodic term $\epsilon\cos(\theta\xi)$. This perturbation breaks the Hamiltonian structure and can induce chaotic motion via homoclinic/heteroclinic orbit intersections near saddle equilibria (Melnikov theory). Numerical integration is performed using the Runge-Kutta 4th-order (RK4) method with step size $h = 0.001$ and absolute/relative tolerance $10^{-8}$ over the interval $\xi \in [0, 300]$. The Lyapunov exponents are computed using the standard Gram-Schmidt orthonormalization algorithm [55]: two tangent vectors are evolved alongside the trajectory, and at each step their lengths are recorded and the vectors re-orthonormalized; the Lyapunov exponents are the time-averaged logarithms of the successive expansion factors. A positive largest exponent $\lambda_1 \gt 0$ confirms chaos. Initial conditions and forcing parameters for each case are specified in the captions of figures 6-8.

It is important to distinguish between the unperturbed ($\epsilon = 0$) and perturbed ($\epsilon \neq 0$) cases. When $\epsilon = 0$, system (51) is a Hamiltonian (integrable) system with conserved energy $H(V,\Phi)$ given by equation (53), and all trajectories are either periodic orbits (around center points) or heteroclinic/homoclinic connections (near saddle points). No chaos is possible in the unperturbed case, as the phase space is entirely organized by these invariant structures. When $\epsilon \neq 0$, the periodic external forcing $\epsilon\cos(\theta\xi)$ breaks the Hamiltonian structure and can cause homoclinic orbit splitting, leading to the formation of chaotic layers around the formerly ordered saddle connections. This is the mechanism by which disorder (chaos) arises: the addition of even a small external periodic disturbance can destroy the integrability and produce exponential sensitivity to initial conditions, as confirmed by the positive Lyapunov exponents in figure 8. To investigate the emergence of chaos in the AMTZ equation, we introduce an external periodic forcing term into the dynamical system derived in section 4. The perturbed system is expressed as where $\epsilon$ in turn is the amplitude of the external perturbation and $\theta$ is the forcing frequency. This perturbation mimics realistic environmental disturbances in experiments [39]. Lyapunov exponents, as quantitative measures of chaotic dynamics, measure the average rate of divergence or convergence in phase space of neighboring trajectories [37]. For example, a positive largest Lyapunov exponent $\lambda_1 \gt 0$ indicates that the behavior of the dynamics may be chaotic and highly sensitive to the initial condition, $\lambda_1 = 0$ is indicative of periodic or quasi-periodic motion, and $\lambda_1 \lt 0$ denoting asymptotic stability. Numerical integration techniques are used to compute the time-dependent Lyapunov exponents. The Poincaré section presents a geometric visualization method for studying dynamic systems, as it involves sampling trajectories at regular intervals in the phase space [38]. For periodic orbits, the Poincaré section is comprised of discrete points, while chaotic attractors are fractal structures with complex patterns. Multi-stability analysis demonstrates how multiple attractors maintain the same initial characteristics under different, though identical, parameters, showing the sensitivity of the system to initial conditions. The chaotic behavior, phase portraits, Poincaré sections, Lyapunov exponents, and multi-stability plots of the systems for distinct parameter sets is shown in figures 6-8, showing that the AMTZ system shows strong chaotic dynamics during external periodic perturbations.

After implementing the IMSSEM, it has shown that it is the most satisfactory method to obtain seventeen true analytical solution families from the AMTZ equation. And it also found a complete and diverse data of seventeen exact analytical product families across hyperbolic, trigonometric, rational function types. Direct back-substitution of the reduced ODE equation (12) (substitution into both PDE (1)-(3)) results in a pure 0 difference returned by mathematical computation software Maple 2024 in every case. These solutions capture specific nonlinear wave phenomena in multi-dimensional dispersive systems, which are highly relevant to the future field of optical fiber communication, plasma physics, and fluid dynamics. Figure 1 gives rise to 3D surface and 2D contour plots of the soliton solutions calculated in equation (35) that is a member of the hyperbolic cosh/sinh family. Figure 1(a) is the result of a bright soliton structure with focused amplitude peak and is propagated through the spatial domain. This corresponds to, in a physical meaning, the optical energy level in a fiber waveguide with the Kerr nonlinearity effectively offsetting the group-velocity dispersion and generating a self-sustaining pulse. The bell-shaped symmetric profile is a trademark of anomalous-dispersion soliton propagation, and it is typical of long-haul optical communication systems in which pulse distortion has to be minimized. Figure 1(b) shows a kink-type soliton of the potential function $X(x,y,t)$, representing a smooth monotonic transition between two modes. This arrangement is comparable to domain walls of ferromagnetic spin systems, in which magnetization spirals continuously between two stable orientations, and also to shock-like forms of magneto-acoustic plasma waves. The sign change in amplitude verifies that the system supports topological excitations (largely robust to small perturbations). Figure 1(c) reveals a dark soliton of the auxiliary field $\Pi(x,y,t)$, appearing as a localized dip in intensity over a finite continuous wave background. In nonlinear optics, dark solitons represent the phase singularities in the normal dispersion regime, which are well established in Bose-Einstein condensates and in photonic lattices. In the $(x,t)$ plane, the corresponding two-dimensional contour charts, figures 1(d)-(f), corroborate the diagonal propagation pathways, confirming that each region is conveying at a constant velocity without deformation and, in turn, is the necessary pre-condition for a correct traveling wave solution. Figure 2 studies the temporal evolution of the solution under equation (35) with five constant time conditions $ t = 1,3,6,8,10 $. Figure 2(a) illustrates that the primary wave amplitude $\Omega(x,1,t)$ consistently undergoes rightward translation, such that the pulse peak is always the same height and width at any time, indicating that dispersion and nonlinearity are in perfect balance in the AMTZ system. Figure 2(b) also indicates that $X(x,1,t)$ propagates at the same speed as $\Omega$, retaining its kink-like nature without any amplitude decay, since the system obeys general principles of conservation. Figure 2(c) indicates $\Pi(x,1,t)$ remains a localized dark-type region over the time sequence as we can show, which reveals that the coupled field dynamics of the AMTZ are completely coherent. That time-evolution invariance physically is the defining property of a soliton: the wave packet is not an imperfect, transient perturbation but a nonlinear, persistent state that persists indefinitely, in the absence of noise and external inputs, suitable for the encoding of information over fiber optic transmission lines. Figure 3 displays the three-dimensional surface and contour plots of the sech-type soliton solutions as predicted from equation (37). Figure 3(a) shows a single bright soliton of $\Omega(x,y,t)$, with pronounced amplitude peak and the fast exponential decay in both directions, which demonstrates the practical appearance of an ultra-short optical pulse in the mode-locked laser cavity with temporal confinement in the femtosecond region. The steepness of gradient on the face of the peak corresponds to very strong nonlinear self-trapping, where amplitude depends very closely on space to narrow a region in an area at the expense of dispersion. Figure 3(b) shows potential function $X(x,y,t)$ with anti-soliton with inverted polarity, which is in the perspective of, say, a hole-type excitation in plasma physics or a defocusing nonlinear medium in optics, with the wave amplitude being a trough not a crest. Figure 3(c) shows a composite wave structure of $\Pi(x,y,t)$ with both a localized positive peak and an extended asymptotic behavior that implies higher-order nonlinear coupling of transverse and longitudinal wave characteristics, based on the $(2+1)$-dimensional nature of the AMTZ system. The two-dimensional contour plots figures 3(d)-(f) confirm that each of these distinct wave structures propagate along well-defined linear features in the $(x,t)$ plane, validating the traveling wave ansatz and the agreement of the solutions with the original PDE system.

Equation (37) at $t = 1,3,6,8,10$. We display in figure 4(a) that $\Omega(x,1,t)$ is characterized simply by a single soliton profile that maps exactly to the scale of the whole cosine spectrum, without changing the amplitude peak and width, and this offers adequate evidence that although the solution looks sharply singular, it is indeed solitonic. Figure 4(b) consistently propagates $X(x,1,t)$, a bell-shaped positive shape that never changes. Figure 4(c) shows that $\Pi(x,1,t)$ maintains its characteristic positive-amplitude localized structure over the total time order. Physically, the thin but time invariant profile of this solution corresponds well with the behavior of ultra-short pulses in high-energy laser systems, where a strong nonlinearity stifles dispersive broadening even for very steep waveforms. The soliton solutions of equation figure 5 give a tight temporal stability confirmation for equations (35) and (37) at three well separated instances within time. Panels figures 5(a)-(c) are for $\Omega$, $X$ and $\Pi$ in equation (35) at $t = 20,50,80$ and panels figures 5(d)-(f) provide the corresponding fields from equation (37) both at the same time instances. The patterns of the wave profiles are visually identical across the three time snapshots: the amplitude, width, and spatial placement of each feature transition at a constant velocity with no deformation or amplitude decay in any panel. Note that this long-range shape preservation in the stretched time period $t\in[20,80]$ is direct evidence that these solutions are stable traveling waves in the formal sense. To a certain extent, this stability derives from a compromise between the non-linear self-phase modulation conditions and the higher-order dispersive conditions of the AMTZ regime, and it appears to show that these methods can be employed as strong carrier signals in optical fiber networks for long periods of propagation. Figure 6 presents the three-dimensional scale lines and two-dimensional contour maps of the tanh-type soliton solutions obtained from equation (39). Figure 6(a) depicts a kink-antikink primary wave $\Omega(x,y,t)$ which is smooth between two asymptotic values, which physically simulates the propagation of a domain-wall excitation in a ferromagnetic medium or a rarefaction wave in a dispersive plasma. The tanh profile is the canonical solution of the nonlinear-wave equations having competing cubic nonlinearity and linear dispersion, and it can be seen that the AMTZ system has the qualitative nonlinear balance phenomena to exhibit. Figure 6(b) shows the potential field $X(x,y,t)$ propagate as a smooth kink with complementary sign, and expresses the boundary between $X$ and $\Omega$ based on the constraint in equation (8). Figure 6(c) reveals the auxiliary field $\Pi(x,y,t)$ in a sech-squared-shaped shape obtained from tailoring the tanh profile, which is the localized energy density pulse accompanying along the kink. The contour plots 6(d)-(f) show the uniform behavior in the $(x,t)$ plane of all three coupled fields on straight plane paths in figures 6(d)-(f). Figure 7 provides temporal stability validation for equation (39) at $t = 50,100,150$ in panels figures 7(a)-(c), and equation (48) at $t = 20,50,80$ in panels figures 7(d)-(f). For equation (39), the tanh-type kink structure of $\Omega$ is rigidly mapped across the whole domain, and the potential and auxiliary fields all pass at the same rate, which further verified full system coherence on a long temporal domain of $\Delta t = 100$. For equation (48), the compound trigonometric soliton $\Pi$ has a periodic-singular shape at $t = 20,50,80$. The amplitude is indeed strikingly the same. In terms of spatial frequency, even the family members of the trigonometric solution retain their structural identity over large propagation distances. Physically, this robustness indicates that these solutions are not just marginally stable. They lie on the part of completely persistent invariant manifolds of the dynamical system, which makes them ideal candidates for use in dispersive wave transmission. Figure 8 uses colors which depict the time evolution at $t = 1,3,6,8,10$. Figure 8(a) shows the tanh-type dark soliton of $\Omega(x,1,t)$ evenly to the right as a function of time and the transition plateau of the two asymptotic states is perfectly horizontal and with the transition width set. This is the signature of a kink soliton with velocity controlled completely by the wave number parameters $p$, $q$, and $v$ across the propagation coordinate $\xi = px + qy - vt$. Figure 8(b) verifies that the potential $X(x,1,t)$ propagates uniformly, whereas figure 8(c) indicates that the localized sigmoidal shape of $\Pi(x,1,t)$ never abates. The steady rightward drift indicates that the group velocity of the nonlinear wave packet is positive and constant, as can be expected considering the sub-luminal propagation of the dispersive optical media.

Contex soliton solutions presented in equation (48), which involves the combination Figure 9(a) shows that $\Omega(x,y,t)$ is a spatially periodic wave of characteristic asymmetry that comes from the interaction of the tangent component and secant components. With its physical shape, this structure reflects a nonlinear periodic wave train of a class common in plasma lattices and photonic crystal waveguides. As a result of this configuration time is imposed by the medium rather than initial circumstances. Figure 9(b) shows that $X(x,y,t)$ is in fact an osmotically kink type spatial modulation system simultaneously with the oscillatory behavior of $\Omega$, showing clearly the rich multi-component dynamics of an affordable, multi-dimensional AMTZ system. Figure 9(c) depicts $\Pi(x,y,t)$, containing very distinct peaks that directly emerge from a unique $\sec^2$ term at its output, and physically correspond to intense-burst peaks, such as high-intensity burst propagation events in a plasma resonance configuration or multi-photon ion formation in a laser-matter interaction. The contour plots figures 9(d)-(f) also show the regular diagonal propagation of the three fields, and this confirms the regular traveling wave nature of this trigonometric solution family. It displays the time evolution of equation (48) with $t = 1,3,6,8,10$. Figure 10(a) shows that the $\Omega$ component moves uniformly with its asymmetric tan$+$sec profile consistent, indicating that the combined trigonometric nonlinearity supports a stable propagating structure. Figure 10(b) shows that $X(x,1,t)$ has a smooth kink evolution which occurs and the transition location will change at constant speed since this is in agreement with the current velocity of the passing wave. Figure 10(c) shows $\Pi(x,1,t)$ in the form of sharp spikes at points across space, which are always shifted uniformly over time as a function of the poles of $\csc$ and $\sec$ functions, whose coordinates rigidly translate to the current wave. Physically, the periodic singularity structure underlying this solution has the property that the wave phenomena in dissipationless plasmas can be modeled at resonance points at very large amplitudes in the wave field which could have a considerable impact on the mechanism of wave to particle interaction and energy transfer in space plasmas. Figure 11 shows the phase portrait analysis derived from the Galilean transformation (51) for the reduced dynamical system. Figure 11(a) is corresponding to Case 1 ($S_1 \lt 0$, $S_2 \lt 0$), where the equilibrium points $(0,0)$ and $(-1,0)$ are saddle points with real eigenvalues $\lambda_{1,2} = \pm 2$ and $(-\frac{1}{2},0)$ is a center point with purely imaginary eigenvalues $\lambda_{1,2} = \pm\sqrt{-2}$. Phases portrait introduces heteroclinic distances between the pair of saddle points in the center region, with the closed orbits around the center corresponding to periodic wave solutions of AMTZ equation. From a physical point of view, the saddle equilibria correspond to unstable wave propagation modes in which any minor perturbation will cause the trajectory to diverge while the center equilibrium represents stable oscillatory wave packets that move endlessly in phase space without energy loss. Figure 11(b) is Case 2 ($S_1 \gt 0$, $S_2 \gt 0$), and the equilibria configuration is reversed with respect to Case 1: the point $(1,0)$ is a center and has purely imaginary eigenvalues $\lambda_{1,2} = \pm 2$i, and $(\frac{1}{2},0)$ is a saddle with real eigenvalues $\lambda_{1,2} = \pm\sqrt{2}$. Based on this complementary structure, we find that the qualitative dynamics depend on the sign of the product $bv$, which determines the direction the non-linear wave energy concentrates on positive or negative amplitude excitations. Figure 11(c) is Case 3 ($v = 0$) where all equilibrium points coalesce at the origin $(0,0)$ into a cuspidal degenerate equilibrium. This case is related to the static limit of the wave equation where the propagation speed is no longer active, and the phase image becomes cusp-like separatrix. In essence, through the lens of quantum mechanics itself, this critical case represents where qualitatively separated wave propagation regimes cross the threshold for an associated bifurcation of the AMTZ system, that is, from the existence of previously isolated equilibria to a degenerate state at a higher order critical point, and this is what determines how signal propagation systems (nonlinear optical fibers) may undergo a wave amplitude collapse.

Figure 12 shows the chaotic dynamics of the perturbed system (55) under the parameter set $a = 1$, $b = -0.1$, $\beta = 1$, $p = 1$, $q = 1$, $v = 1$, $\theta = 3$, $\epsilon = 3.6$, for initial condition $(0.91,\,0.99)$. We see the two-dimensional chaotic phase portrait in $(V,\Phi)$ plane with a characteristic attractor with an inscrutable dense, non-repeating trajectory that runs through a bounded region without connecting to a stationary point or a limit cycle figure 12(a). The nonlinear, crossing pattern is a signature of deterministic chaos: the system progresses according to deterministic equations, but its trajectories are exponentially sensitive to the initial state, which makes long-term prediction in practice impossible. Physically, this indicates that the smallest variation on the initial amplitude or phase of the optical pulse will cause quite different wave propagation effects, after sufficient propagation distances, an interesting problem in the design of excellent optical communication systems. Figure 12(b) gives a 3D phase view at $(V,\Phi,W)$ space with $W$ as the time-extended coordinate revealing toroidal geometry of the strange attractor and demonstrating that chaotic dynamics are not a two-dimensional character but rather a real 3D phenomenon with fractal dimension greater than 2. Figure 12(c) is the Poincaré section which resulted through sampling the trajectory over regular phase intervals: non-periodic and scattered distribution of points around the section confirms an absence of regular orbit and geometrically verify their fractal structure. On the other hand, a periodic orbit would only form a finite set of isolated points on the Poincaré section. Figure 12(d) is the three-dimensional Poincaré map with phase-angle colored phases, and we are seeing that the color gradient from blue (low phase) to red (high phase) results in the complex multi-layer structure of the attractor in the extended phase space and confirms that the chaotic motion has a well-defined geometric structure, consistent with a strange attractor, not random noise. Figure 13 has a second chaotic characterization under another parameter regime: $a = 1$, $b = -0.5$, $\beta = 1$, $p = 1$, $q = 1$, $v = 2$, $\theta = -2$, $\epsilon = 0.16$, with a condition of $(0.7,\,0.6)$. Whilst the perturbation amplitude $\epsilon = 0.16$ is much smaller than $\epsilon = 3.6$ of the previous case, the AMTZ system has stable chaotic behavior, indicating that the chaos in this system is not a high-perturbation artifact, but a structurally persistent phenomenon. Figure 13(a) depicts a figure-eight-shaped exotic attractor in $(V,\Phi)$ plane, a geometry that can be physically associated with an unstable equilibrium neighborhood, where a steady state pattern in phase space oscillates between the two unstable neighborhoods in the oscillating wave trajectory, and which corresponds with the homoclinic orbit splitting projected by Melnikov theory near the saddle equilibria described in section 4. Figure 13(b) shows that the three-dimensional size of the attractor is a twisted toroidal structure. The fractal attractor is found to survive the reduced amplitude forcing. Figure 13(c) shows a Poincaré section with widely scattered points, establishing a clear geometric evidence of chaos: the irregular distribution tells us there are no periodic orbit under these forcing conditions. The 3D Poincaré map is illustrated in figure 13(d) with phase coloring, the random arrangement of colored spots in the section confirms that the system is ergodically sampling all the phases, which is a key characteristic of chaotic attractors. Quantitative confirmation that the system is chaotic is provided by Lyapunov exponent analysis and multi-stability plots (figure 14). Figure 14(a) presents time-step transition of two great Lyapunov exponents $\lambda_1$ and $\lambda_2$ in figure 12: after a transient $\lambda_1$ reaches a positive value (approximately $\lambda_1 \approx 0.5 \gt 0$), but $\lambda_2$ shows a negative value implying good quantitative confirmation of chaos. In phase space, the largest Lyapunov exponent in positive direction implies that trajectories near each one of the stations separate exponentially faster with the difference e$ ^{\lambda_1 t}$, and therefore the system is unpredictable to a time of Lyapunov $\tau_L = 1/\lambda_1$. Figure 14(b) depicts the Lyapunov exponents for figure 13 in which, while the forcing amplitude $\epsilon = 0.16$ is smaller, the highest Lyapunov exponent is positive confirming that no matter how small the perturbation, the chaos can be persistent, which is directly relevant to real-world optical systems that always experience small environmental perturbations. Figure 14(c) shows a multi-stability diagram for figure 12, where two trajectories originating from slightly different initial conditions $V(t)$, $\Phi(t)$ gradually diverge into qualitatively different attractors at $t \in [0,200]$, representing visually with a dramatic representation the delicacy of dependence on initial conditions. Figure 14(d) displays multi-stability for figure 13 parameters and the resulting divergence is similar over the time period $t \in [0,300]$, indicating that the intricate basin structure and co-existing attractors are a robust part of the AMTZ system in response to periodic external force. The IMSSEM is much more effective than classic ones including the extended tanh method and the $(G ^{^{\prime}}/G)$ expansion method as it achieves production of all 17 solution families simultaneously from the same algebraic system in this case rather than each function type having to run in a separate manner. All solutions obtained satisfy the original PDE system (1)-(3) and the reduced ODE equation (12) as such, and Maple 2024 back substitution confirmed the zero residuals in every case that was obtained, adding considerable mathematical validation. The comparative table 2 confirms that neither the Myrzakulov-I nor the Myrzakulov-II equation has been investigated via bifurcation or chaos analysis in the existing literature, thus providing the first step in the comprehensive dynamical characterization for the whole family of equations by the Myrzakulov family.

This study presented the first comprehensive analytical and dynamical characterization of the AMTZ equation. Using the IMSSEM with the auxiliary ODE $G^{^{\prime} 2} = CG^4 + BG^2 + A$, seventeen families of exact solutions were derived, including hyperbolic (sech, csch, tanh, coth), trigonometric (tan, cot, sec, csc), and rational function types. Bifurcation analysis of the reduced dynamical system identified three equilibrium configurations-saddle, center, and cuspidal points-whose classification is governed by the sign of the Jacobian determinant $J(V,\Phi) = -4b\beta^2p^2q^2a(6V^2p^2 + 6Vbpv + b^2v^2)$. Temporal stability analysis confirmed shape-preserving soliton propagation over long distances. Chaos analysis under external periodic forcing $\epsilon\cos(\theta\xi)$ revealed positive Lyapunov exponents and fractal Poincaré sections, confirming that the integrable AMTZ system exhibits strong sensitivity to perturbations. The key findings, limitations, and future directions are detailed in section 7.