EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6723 ISSN 1307-5543 – ejpam.com Published by New York Business Global Frequency Locking and Bifurcation Analysis in Asymmetrically Forced Van der Pol-Duffing Model of Glacial Cycles Ibrahim Alraddadi 1,∗, Saad M. Almuaddi 2,3,∗ 1 Department of Mathematics, Faculty of Science, Islamic University of Madinah, Madinah, Saudi Arabia 2 Basic Applied Scientific Research Center, Imam Abdulrahman Bin Faisal University, P.O. Box 1982, Dammam 31441, Saudi Arabia 3 Mathematics Department, College of Science, Imam Abdulrahman Bin Faisal University, Dammam 31441, Saudi Arabia Abstract. This paper examines the application of self-sustained oscillator systems, particularly the van der Pol-Duffing oscillator, to understand the complex dynamics of Pleistocene glacial cy- cles. We investigate how asymmetric forcing (β ̸= 0) influences frequency locking and bifurcation behavior through geometric singular perturbation theory (GSPT) analysis and Poincaré return map construction. Our numerical results demonstrate that the van der Pol-Duffing oscillator pos- sesses substantially larger regions of stable periodic behavior in parameter space compared to standard van der Pol oscillators. As asymmetry increases from β = 0.25 to β = 1.2, we observed progressive narrowing of Arnold tongue structures with most frequency locking regions requiring stronger forcing amplitudes (a ≥ 1.5) to initiate synchronization. However, remarkably resilient 2:1 frequency locking regions persist across all asymmetry levels. This provides a mathematical frame- work for explaining dominant frequency transitions observed in paleoclimate records, particularly the Mid-Pleistocene Transition from 41 kyr to 100 kyr glacial cycles (4-significance) 2020 Mathematics Subject Classifications: Primary: 4C23, 37C25, 37-04, 37C35, 34C25, 34C26 Key Words and Phrases: Dynamical system, Periodic solutions, bifurcations 1. Introduction The Pleistocene glacial record presents a fundamental puzzle in climate dynamics that has captivated researchers for decades [1]. While astronomical forcing is dominated by spectral components at periods of 41 kyr (obliquity) and 19-23 kyr (precession), the climate ∗Corresponding author. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6723 Email addresses: ialraddadi@iu.edu.sa (Ibrahim Alraddadi), smuaddi@iau.edu.sa (Saad M. Almuaddi) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 2 of 20 response during the late Pleistocene exhibits dominant periodicity at approximately 100 kyr. Furthermore, the system underwent a significant transition approximately 1 Myr ago—the Mid-Pleistocene Transition (MPT)—shifting from lower-amplitude 41 kyr cycles to higher-amplitude 100 kyr cycles [2]. This mismatch between forcing and response frequencies, combined with the MPT bifurcation, indicates that Earth’s climate system responds nonlinearly to astronomical forcing. Linear models are insufficient to explain these phenomena, necessitating exami- nation of nonlinear oscillator systems with rich dynamical behavior [1, 3]. The van der Pol oscillator, originally introduced by van der Pol (1926) [4] to model electrical circuits, has emerged as a powerful framework for understanding self-sustained oscillations in climate systems. Van der Pol and van der Mark (1927, 1928) [5] first ob- served frequency locking phenomena—now known as synchronization—when examining the response of relaxation oscillators to periodic forcing. This phenomenon, where the os- cillator’s frequency becomes locked to rational multiples of the forcing frequency, provides a crucial mechanism for understanding glacial cycle dynamics [1, 3]. Modern applications of van der Pol oscillators to climate science began with Crucifix (2012) [3], who used forced van der Pol models to examine astronomical forcing and asym- metry between ice-forming and melting phases during the late Pleistocene. Subsequent work by de Saedeleer et al. (2013) [1] and Ditlevsen and Ashwin (2018) [2] has extended these models to explain the MPT as a synchronization transition rather than a change in external forcing while contemporary studies have explored Arnold tongue structures in various nonlinear systems [6] Advanced bifurcation analysis techniques have been applied to oscillatory systems [7], with computational methods for nonlinear dynamics expanding rapidly [8]. Recent advances in geometric singular perturbation theory (GSPT) have provided rig- orous mathematical foundations for analyzing these slow-fast climate systems [9]. The work of Guckenheimer et al. (2003) [9] established the theoretical framework for under- standing folded singularities and canard trajectories in forced van der Pol systems, while Ashwin et al. (2018) [10] demonstrated that van der Pol-Duffing variants exhibit enhanced chaotic behavior crucial for climate modeling. 2. Nonlinear Oscillator Models for Glacial Cycles 2.1. Theoretical Framework Autonomous nonlinear oscillator systems exhibit several mathematical properties that make them particularly relevant for understanding glacial cycle dynamics [3, 10]. These systems naturally produce self-sustained oscillations even in the absence of external forc- ing, creating limit cycles that can serve as baseline climate states. When subjected to external periodic forcing, they demonstrate n : m phase-locking behavior, where the ratio of response frequency ω to forcing frequency ωF satisfies ω/ωF = n/m with n,m ∈ Z+. This frequency locking mechanism provides a mathematical foundation for explaining how climate systems can respond at frequencies different from those of astronomical forcing. I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 3 of 20 The temporal asymmetry observed in glacial cycles—with slow ice accumulation fol- lowed by rapid deglaciation—emerges naturally from these oscillators through their gen- eration of non-sinusoidal waveforms. Unlike linear systems that preserve the spectral characteristics of their forcing, nonlinear oscillators can transform symmetric forcing into asymmetric responses that mirror the observed glacial-interglacial patterns. Parameter- dependent bifurcations in these systems correspond to observed climate transitions, pro- viding a mechanism for understanding sudden changes in system behavior such as the Mid-Pleistocene Transition. Perhaps most importantly for climate applications, these oscillators can exhibit chaotic dynamics characterized by positive Lyapunov exponents under quasi-periodic forcing [10] demonstrated through numerical computation that several low-order nonlinear models, particularly those with nonlinear restoring forces, exhibit positive Lyapunov exponents under quasi-periodic forcing across significant regions of parameter space. This result implies fundamental limits to the predictability of trajectories in such systems despite deterministic forcing, which aligns with the irregular timing variations observed in actual glacial records. The van der Pol-Duffing oscillator provides a particularly powerful framework because it incorporates nonlinear damping from the van der Pol equation, nonlinear restoring forces from the Duffing equation, and parametric forcing terms [11–13]. This combination creates a system capable of exhibiting the complex dynamics observed in paleoclimate records while maintaining sufficient mathematical tractability for detailed analysis. 2.2. The Modified Van der Pol Oscillator We consider the modified van der Pol oscillator proposed as a low-order model for ice-age cycles [10] : τ2κ2 d2y dt2 − ατκ(1− y2) dy dt + g(y)− γF (t) + β = 0 (1) where the nonlinear restoring force in the van der Pol-Duffing variant is: g(y) = −y + 1.2y3 (2) This differs significantly from the standard van der Pol oscillator where g(y) = y. The cubic nonlinearity introduces multiple equilibria and enhanced shear in phase space, essential ingredients for robust chaotic behavior [10]. Parameter Interpretation: • α = 11.11: Controls fast/slow timescale separation • β: Symmetry-breaking parameter (β = 0 for symmetric case) • γ: Effective forcing amplitude • κ = 52.0: Sets unforced oscillation timescale I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 4 of 20 • τ = 1: Time scaling factor (τ = 1 gives ∼100 kyr unforced period) Using the Liénard transformation x = y − y3/3 − y/(τκα), we obtain the first-order system [9]: x′ = 1 τκ (γk sin(2πωt)− β − g(y)) (3) y′ = α τκ (y − y3/3 + x) (4) (a) (b) Figure 1: (a) Phase portrait of system (3-4) demonstrating different dynamical regimes with parameter values given above. The trajectories show the relaxation oscillation pattern with distinct slow and fast phases of motion. (b) The trajectories show the relaxation oscillation pattern with distinct slow and fast phases of motion Figure 1 illustrates the rich dynamical behavior of system (3-4). The phase portraits reveal the characteristic features of van der Pol duffing oscillators including the S-shaped nullclines and the alternating slow-fast motion that creates relaxation oscillations. Setting τκ = 1 and defining ε = 1/α, a = γ, θ = ωt, we get the non-autonomous slow-fast system: x′ = a sin(2πθ)− g(y)− β (5) εy′ = y − y3/3 + x (6) θ′ = ω (7) I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 5 of 20 3. Geometric Singular Perturbation Analysis 3.1. Critical Manifold Structure When 0 < ε ≪ 1, system (5)–(7) exhibits slow-fast dynamics with x and θ as slow variables and y as the fast variable. Following the GSPT framework of Guckenheimer et al [9], we analyze the singular limit ε = 0. The critical manifold is defined as: S = {(x, y, θ) ∈ R3 | y − y3/3 + x = 0} (8) This gives x = y3/3 − y, defining an S-shaped curve in the (x, y) plane. The critical manifold has: • Repelling sheet: Sr = S ∩ {−1 < y < 1} where ∂/∂y(y − y3/3) = 1− y2 > 0 • Attracting sheets: Sa = S ∩ {y < −1} ∪ S ∩ {y > 1} where ∂/∂y(y − y3/3) = 1− y2 < 0 Fold lines occur at y = ±1: L− = {(−2/3,−1, θ) : θ ∈ (0, 2πn)} (9) L+ = {(2/3, 1, θ) : θ ∈ (0, 2πn)} (10) 3.2. Slow and Fast Subsystems Slow subsystem (ε = 0): x′ = a sin(2πθ) + y − 1.2y3 − β (11) 0 = y − y3/3 + x (12) θ′ = ω (13) Fast subsystem (layer problem): dx dT = 0 (14) dy dT = y − y3/3 + x (15) dθ dT = 0 (16) where T is fast time (t = εT ). I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 6 of 20 3.3. Desingularized Reduced System The duffing nonlinearity term g(y) = −y + 1.2y3 significantly affects global dynamics compared to the standard van der Pol oscillator [10]. Differentiating the critical manifold condition and applying the chain rule: ∂f ∂y y′ + ∂f ∂x x′ + ∂f ∂θ θ′ = 0 (17) With ∂f/∂y = 1− y2, ∂f/∂x = 1, ∂f/∂θ = 0: (1− y2)y′ = −(a sin(2πθ) + y − 1.2y3 − β) (18) Rescaling time by t = (y2 − 1)s yields the desingularized reduced system: dθ ds = ω(y2 − 1) (19) dy ds = a sin(2πθ) + y − 1.2y3 − β (20) (a) symmetric case (b) asymmetric case Figure 2: (a) Trajectories of the desingularized reduced system (19-20) for symmetric case β = 0 with parameters a = 0.55 and ω = 1.5. The phase portrait shows four folded singularities on the fold lines y = ±1, with blue trajectories on the upper stable branch and red trajectories on the lower stable branch. (b) Trajectories of the desingularized reduced system (19-20) for asymmetric case β = 0.5 with same parameters a = 0.55 and ω = 1.5. The phase portrait shows two folded singularities on the fold lines y = −1, with blue trajectories on the upper stable branch and red trajectories on the lower stable branch. I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 7 of 20 (a) β = 0.25 (b) β = 0.5 Figure 3: (a) Trajectories of the desingularized reduced system (19-20) for β = 1.2 with parameters a = 3.5 and ω = 1.5. The phase portrait shows four folded singularities on the fold lines y = ±1, with blue trajectories on the upper stable branch and red trajectories on the lower stable branch. (b) Trajectories of the desingularized reduced system (19-20) for β = 1.2 with parameters a = 3.5 and ω = 1.5. The phase portrait shows four folded singularities on the fold lines y = ±1, with blue trajectories on the upper stable branch and red trajectories on the lower stable branch. Notes some blue trajectories attract to stable folded singularities equilibrium point on lower stable branch I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 8 of 20 (a) β = 0.8 (b) β = 1.2 Figure 4: (a)Trajectories of the desingularized reduced system (19-20) for β = 0.8 with parameters a = 3.5 and ω = 1.5. The system exhibits similar folded singularity structure but with altered positions due to the different symmetry-breaking parameter value. (b) Trajectories of the desingularized reduced system (19-20) for β = 1.2 with parameters a = 3.5 and ω = 1.5. The phase portrait shows four folded singularities on the fold lines y = ±1, with blue trajectories on the upper stable branch and red trajectories on the lower stable branch. Figures 2 to 4 illustrate the evolution of system (19-20) as the symmetry-breaking parameter β varies from 0 to 1.2. The phase portraits reveal how β systematically af- fects the position and stability of folded singularities on the fold lines y = ±1. As β decreases, the system transitions from highly asymmetric behavior toward the symmetric case, demonstrating the crucial role of this parameter in controlling the global dynamics of the desingularized reduced system. 3.4. Folded Singularities Folded singularities occur where dy/ds = 0 and y = ±1. For y = 1: sin(2πθ) = β − 0.8 a (21) For y = −1: sin(2πθ) = β − 0.2 a (22) The eigenvalues at folded singularities determine their type: λ1,2 = 1− 3.6y2 ± √ (1− 3.6y2)2 + 16ωπa cos(2πθ) 2 (23) I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 9 of 20 4. Poincaré Return Map Analysis 4.1. Construction of Return Map Following Guckenheimer et al. (2003) [9], we construct a Poincaré return map F : S2 → S2 from the section y = 2 to itself. The map decomposes as: F = J− ◦ P− ◦ J+ ◦ P+ (24) where: • P+ : S2 → S1 (flow on stable manifold from y = 2 to y = 1) • J+ : S1 → S−2 (jump map from fold line L+ to y = −2) • P− : S−2 → S−1 (flow on stable manifold from y = −2 to y = −1) • J− : S−1 → S2 (jump map from fold line L− to y = 2) The jump maps are explicit: J+(θ, 1) = (θ,−2) and J−(θ,−1) = (θ, 2), while P+ and P− are computed by numerical integration of the desingularized system (19)–(20) (MATLAB code is shown in [14]) 4.2. Poincaré Return Map Analysis The construction of Poincaré return maps reveals fundamental differences between symmetric and asymmetric systems. In the symmetric Van der Pol oscillator, the return map F : S2 → S2 exhibits classical circle map behavior for small forcing amplitudes, with smooth, invertible dynamics. As forcing amplitude increases beyond a = 1, the map becomes non-invertible with gap regions corresponding to folded equilibria, but these gaps appear symmetrically. (author?) [15] showed that the asymmetric Van der Pol oscillator generates return maps with fundamentally different structures. The asymmetry parameter β shifts the positions of discontinuities and creates asymmetric gap regions in the return map. These gaps occur when trajectories cannot escape from regions between folded saddle equilibria and their unstable manifolds, but their positions and sizes depend on β in a non-trivial manner. 5. Numerical Results and Discussion 5.1. Return Map Structure Poincaré return maps constructed from the section y = 2 exhibit markedly different characteristics depending on the forcing amplitude a (Figure 5). For a < 1, the maps are invertible circle maps with smooth, monotonic structure that preserve topological properties under iteration. These smooth maps typically exhibit quasiperiodic dynamics or simple periodic orbits. I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 10 of 20 When a > 1, the return maps develop discontinuities and gap regions where trajectories cannot return to the Poincaré section (Figure 6). These gaps arise from folded saddle equilibria that trap trajectories in certain phase space regions, preventing their escape to complete full orbits. The resulting non-invertible maps can exhibit complex bifurcation cascades and chaotic dynamics. (a) symmetric case (b) asymmetric case Figure 5: First return map of the desingularized reduced system (19-20) for the first parameter set. The red curve represents the return map F (θ) and the black diagonal line shows θ. (a) First return map showing symmetric and discontinuity when β = 0. (b)The intersection points indicate fixed points of the return map, demonstrating asymmetric case and continuity appear here. I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 11 of 20 (a) (b) (c) Figure 6: (a)First return map showing saddle-node bifurcation behavior for the first pa- rameter set. The red curve represents the return map F (θ) and the black diagonal line shows θ. The intersection points indicate fixed points of the return map, demonstrating the pre-bifurcation state with stable periodic solutions. (b) First return map at the saddle- node bifurcation point. The return map curve F (θ) becomes tangent to the diagonal line, indicating the critical parameter value where the saddle-node bifurcation occurs and peri- odic solutions are created or destroyed. (c) First return map in the post-bifurcation regime showing the disappearance of fixed points. The return map curve no longer intersects the diagonal line at the previous locations, demonstrating the completion of the saddle-node bifurcation process. 5.2. Bifurcation Analysis The saddle-node bifurcations that define Arnold tongue boundaries occur at different parameter values. For the symmetric case, bifurcations occur when the return map be- comes tangent to the diagonal line θn+1 = θn at symmetric positions. In the asymmetric I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 12 of 20 case, these tangency conditions shift according to the folded singularity positions, creat- ing asymmetric bifurcation curves in parameter space. Saddle-node bifurcations marking Arnold tongue boundaries satisfy the analytical condition given by equations (19)-(20), where the parameter β depends on which fold line contains the bifurcating periodic orbit. The asymmetric parameter β shifts these bifurcation boundaries according to the folded singularity conditions in equations (19) and (20). This systematic displacement of Arnold tongues provides the mechanism for frequency transitions observed during the Mid-Pleistocene climate transition, where gradually changing system parameters can move the climate system between different synchronization regimes [16]. Figures 6(a) through 6(c) illustrate the saddle-node bifurcation of the return map as system parameters are varied. These figures show the classical saddle-node bifurcation sequence where a pair of fixed points (one stable, one unstable) are created and subse- quently annihilated as the parameter crosses the bifurcation threshold. The saddle-node bifurcation of the return map corresponds to the boundaries of Arnold tongues in param- eter space, marking the transition between frequency-locked and non-frequency-locked behavior in the original dynamical system. 5.3. Arnold Tongues and Frequency Locking Properties The most significant difference between symmetric and asymmetric systems lies in their Arnold tongue structures and frequency locking properties. In the symmetric Van der Pol oscillator, Arnold tongues exhibit classical symmetric structure with well-defined boundaries determined by saddle-node bifurcations. The tongue widths scale with forc- ing amplitude, and the system exhibits clear 1:1, 2:1, and 3:1 frequency locking regions that explain the transition from 41kyr to 100kyr glacial cycles. [15] demonstrated that the asymmetric forced Van der Pol oscillator (β ̸= 0) exhibits fundamentally different Arnold tongue behavior. As β increases from 0 to 1.2, the frequency locking regions be- come progressively narrower, with the 1:1 tongue showing the most dramatic reduction in width. This narrowing occurs because the asymmetry parameter shifts the effective resonance conditions, making synchronization more difficult to maintain across parameter variations. The parameter space analysis reveals that for β = 0.25, large Arnold tongues persist for periods 2, 3, 4, and 8, maintaining robust frequency locking similar to the sym- metric case. However, as β increases to 0.8 and beyond, the tongues become significantly narrower, with periods ranging from 1 to 25 but occupying smaller parameter regions. This progressive narrowing reflects the increasing difficulty of maintaining synchronization as system asymmetry grows. The Van der Pol-Duffing system shows even more dramatic changes in Arnold tongue structure due to the combined effects of the cubic nonlinear- ity and asymmetry parameter [17] found that while individual tongues become narrower, the overall parameter space covered by chaotic regions increases substantially, providing a more realistic representation of the irregular timing variations observed in paleoclimate records. I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 13 of 20 5.4. Validation and Cross-Verification Cross-validation with symmetric case (β = 0): Our results reproduce classical Van der Pol-Duffing behavior when β = 0, matching Guckenheimer et al. (2003) canard solutions and exhibiting expected 1:1, 2:1, and 3:1 Arnold tongue structures. Physical parameter mapping: Forcing frequency ω = 1 corresponds to 100 kyr eccentricity cycles, ω = 2.44 to 41 kyr obliquity cycles, and amplitude γ = 0.1-2.0 spans realistic eccentricity variations (0.005-0.07). Paleoclimate proxy comparison: The progressive Arnold tongue narrowing with increasing β qualitatively matches the observed reduction in glacial cycle regularity during the Mid-Pleistocene Transition (1.2-0.8 Ma), supporting the model’s physical relevance. 5.5. Frequency Locking Regions Analysis The (a, ω) parameter space reveals distinct dynamical regimes as the symmetry-breaking parameter β varies (Figure 7). For the symmetric case (β = 0), we observe well-defined Arnold tongue structures with large 1:1 frequency locking regions and significant bista- bility zones between different periodic windows. The tongue boundaries correspond to saddle-node bifurcations of periodic orbits, consistent with classical synchronization the- ory. I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 14 of 20 beta=0.25 0 1 2 3 4 a 0 1 2 3 4 0 10 20 30 2:1 4:1 3:1 (a) beta=0.5 0 1 2 3 4 a 0 1 2 3 4 0 10 20 30 2:1 3:1 4:1 (b) beta=0.8 0 1 2 3 4 a 0 1 2 3 4 0 10 20 30 2:1 8:1 3:1 4:1 (c) beta=1.2 0 1 2 3 4 a 0 1 2 3 4 0 5 10 15 20 25 30 2:1 3:1 8:1 (d) Figure 7: Arnold tongues in the (a, ω) parameter plane for different β. The colorbar indicates the rotation period N of the return map. (a) β = 0.25 shows wide frequency locking regions around ω = 1 and ω = 4 with clear period-2, period-3, period-4 and period-8 structures. (b) β = 0.5 displays one dominant wide locking region around ω = 1, particularly prominent for small forcing amplitudes a. (c) β = 0.8 exhibits narrowed Arnold tongues with most frequency locking regions starting from a = 1, except one resilient region 2:1 from a = 0. (d) β = 1.2 demonstrates severely reduced tongues where most locking regions require a ≥ 1.5 to initiate, while only one wide frequency locking region 2:1 persists from a = 0 around ω ≈ 2 (MATLAB code is shown in [14]) I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 15 of 20 These results have significant implications for understanding climate system dynam- ics and the mechanisms underlying major climate transitions such as the Mid-Pleistocene Transition (MPT) [16]. Our parameter space analysis reveals a clear hierarchy of syn- chronization robustness as the asymmetry parameter β increases from 0.25 to 1.2. In the near-symmetric regime (β = 0.25), Figure 7(a) demonstrates that the system maintains wide frequency locking regions around ω = 1 and ω = 4, supporting multiple periodic structures including period-2, period-3, and period-4 oscillations. This robust synchro- nization capability suggests that climate systems with minimal asymmetry can maintain stable frequency relationships with astronomical forcing across broad parameter ranges, consistent with the relatively stable 41 kyr obliquity-paced cycles observed during the early Pleistocene [16] As asymmetry increases to β = 0.5, Figure 7(b) shows a fundamental reorganization where the system consolidates into a single dominant wide locking region around ω = 1, particularly prominent for small forcing amplitudes a. This consolidation represents a critical transition point where the system begins to lose its capacity for diverse frequency locking modes while maintaining strong synchronization in specific parameter windows. The preference for lower forcing amplitudes in this regime suggests that moderately asym- metric climate systems may be more sensitive to weak astronomical forcing components [18] The transition to strongly asymmetric behavior becomes evident at β = 0.8, where Figure 7(c) reveals significantly narrowed Arnold tongues with most frequency locking regions requiring a ≥ 1 to initiate. Notably, one resilient region persists from a = 0, indicating that even in highly asymmetric systems, certain synchronization modes remain accessible under weak forcing conditions. This threshold effect has important implications for understanding how ice sheet dynamics and climate feedbacks may influence the system’s response to orbital forcing. The most dramatic changes occur in the strongly asymmetric regime at β = 1.2, where Figure 7(d) demonstrates severe reduction in frequency locking capabilities. Most locking regions now require a ≥ 1.5 to initiate, representing a substantial increase in the forcing threshold needed to achieve synchronization. Remarkably, one wide frequency locking region 2:1 around ω ≈ 2 remains resilient from a = 0, suggesting this represents the most robust synchronization mode in highly asymmetric systems. 5.6. Physical Interpretation and Climate Implications Mechanism for Arnold tongue narrowing As asymmetry β increases, folded saddle equilibria shift systematically (Equations 19- 20), displacing bifurcation boundaries. This creates an asymmetric response to orbital forcing where glacial inception becomes increasingly sensitive to initial conditions. Mid- Pleistocene Transition mechanism: Our results suggest the 41 kyr → 100 kyr transition resulted from gradual increase in climate system asymmetry (increasing β) rather than a threshold effect. As β approached 0.8-1.0, the 2:1 frequency locking region (41 kyr response to 20.5 kyr precession) destabilized, allowing 1:1 eccentricity locking (100 kyr) I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 16 of 20 to dominate. 6. Conclusion This study has revealed fundamental new insights into the dynamical behavior of the van der Pol-Duffing oscillator through comprehensive numerical analysis of Poincaré return maps and Arnold tongue structures. Our numerical results demonstrate that the forcing amplitude parameter a serves as a critical bifurcation parameter, producing a dramatic structural transition in return map behavior: for a < 1, we obtained smooth, invertible circle maps with monotonic structure that preserve topological properties and exhibit quasiperiodic dynamics or simple periodic orbits, while for a > 1, the return maps develop discontinuities and gap regions arising from folded saddle equilibria that trap trajectories and prevent complete orbit closure, leading to non-invertible maps capable of complex bi- furcation cascades and chaotic dynamics. Most significantly, our systematic investigation of the asymmetry parameter β revealed a progressive deterioration of frequency locking capabilities as β increases from 0 to 1.2. For the near-symmetric case (β = 0.25), we found robust wide frequency locking regions around ω = 1 and ω = 4 supporting period-2, period-3, period-4, and period-8 structures across broad parameter ranges. At moderate asymmetry (β = 0.5), the system consolidated into a single dominant wide locking region around ω = 1, particularly prominent for small forcing amplitudes. The critical transi- tion occurred at β = 0.8, where Arnold tongues became significantly narrower with most frequency locking regions requiring a ≥ 1 to initiate, though one resilient 2 : 1 region persisted from a = 0. In the strongly asymmetric regime (β = 1.2), we observed severe reduction in frequency locking capabilities with most locking regions requiring a ≥ 1.5 to initiate, while remarkably, one wide 2 : 1 frequency locking region around ω ≈ 2 re- mained resilient from a = 0, representing the most robust synchronization mode in highly asymmetric systems. Our saddle-node bifurcation analysis revealed that these structural changes occur through systematic displacement of Arnold tongue boundaries according to folded singularity conditions, providing the mechanism for frequency transitions observed during major climate transitions. The numerical evidence shows that while individual Arnold tongues become progressively narrower with increasing asymmetry, the overall pa- rameter space covered by chaotic regions increases substantially, offering a more realistic representation of the irregular timing variations observed in paleoclimate records and es- tablishing enhanced chaotic behavior as a distinguishing feature of asymmetric van der Pol-Duffing systems compared to their symmetric counterparts. Future research directions should focus on several key areas to extend the asymmet- ric Van der Pol-Duffing framework. First, incorporating stochastic perturbations repre- senting high-frequency climate variability could provide more realistic ensemble forecasts and better uncertainty quantification. Second, coupling multiple oscillators with different timescales could capture the interaction between orbital forcing, ice sheet dynamics, and atmospheric composition changes. Third, systematic parameter estimation from high- resolution paleoclimate records would improve model calibration and validation. Fourth, I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 17 of 20 developing hybrid approaches that combine low-dimensional conceptual models with de- tailed Earth system components could bridge the gap between mathematical tractability and physical realism. Finally, extending the framework to other Earth system oscillations such as monsoon variability and ocean circulation patterns could demonstrate the broader applicability of these nonlinear dynamics principles. Our two-variable oscillator cannot capture ice sheet-atmosphere-ocean interactions ex- plicitly, requiring parameterization through effective forcing terms. also, Results are sen- sitive to α and β values variations can shift Arnold tongue boundaries Author Contributions I. Alraddadi: Conceptualization, methodology, investigation, data analysis, writing—original draft, supervision. S.M. Almuaddi: Formal analysis, validation, software, resources, data curation, writing—review and editing. All authors have read and agreed to the published version of the manuscript. References [1] Bernard de Saedeleer, Michel Crucifix, and Sebastian Wieczorek. Is the astronomical forcing a reliable and unique pacemaker for climate? A conceptual model study. Climate Dynamics, 2013. [2] Peter Ditlevsen and Peter Ashwin. Complex climate response to astronomical forcing: The middle-Pleistocene transition in glacial cycles and changes in frequency locking. pages 1–13, 2018. [3] Michel Crucifix. Oscillators and relaxation phenomena in Pleistocene climate theory. Trans. R. Soc. A, 370:1140–1165, 2012. [4] Balth Van der Pol. On “relaxation-oscillations”. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11):978–992, 1926. [5] Balth Van der Pol and Jan Van Der Mark. Frequency demultiplication. Nature, 120(3019):363–364, 1927. [6] Nikolay Kyurkchiev, Tsvetelin Zaevski, Maria Vasileva, Vesselin Kyurkchiev, Anton Iliev, and Asen Rahnev. Dynamics of a class of extended duffing–van der pol os- cillators: Melnikov’s approach, simulations, control over oscillations. Mathematics, 13(14), 2025. [7] Wieslaw Marszalek and Maciej Walczak. Bifurcation diagrams of nonlinear oscillatory dynamical systems: A brief review in 1d, 2d and 3d. Entropy, 26(9), 2024. [8] O. Cornejo-Pérez, P. Albares, and J. Negro. Solutions of an extended duffing–van der pol equation with variable coefficients. Physica D: Nonlinear Phenomena, 476:134675, 2025. [9] John Guckenheimer, Kathleen Hoffman, and Warren Weckesser. The Forced van der Pol Equation I: The Slow Flow and Its Bifurcations. SIAM J. Applied Dynamical Systems, 2(1):1–35, 2003. I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 18 of 20 [10] P. Ashwin, C. D. Camp, and A. S. von der Heydt. Chaotic and nonchaotic response to quasiperiodic forcing: limits to predictability of ice ages paced by Milankovitch forcing. Dynamics and Statistics of the Climate System, 3(1):dzy002, 2018. [11] K. Bold, C. Edwards, J. Guckenheimer, S. Guharay, K. Hoffman, J. Hubbard, R. Oliva, and W. Weckesser. The forced van der Pol equation II: Canards in the reduced system. SIAM Journal on Applied Dynamical Systems, 2(4):570–608, 2003. [12] G. M. Moatimid and T. S. Amer. Dynamical system of a time-delayed 6-Van der Pol oscillator: a non-perturbative approach. Scientific Reports, 13:11942, 2023. [13] M. A. Elfouly and M. A. Sohaly. Van der Pol model in two-delay differential equation representation. Scientific Reports, 12:2925, 2022. [14] Ibrahim Alraddadi. MATLAB code for computing the first return map and fre- quency locking. https://doi.org/10.5281/zenodo.17081498, sep 2025. Zenodo, DOI: 10.5281/zenodo.17081498. [15] Alraddadi. The asymmetric periodically forced van der pol oscillator. European Journal of Pure and Applied Mathematics, 18(1):5787, 2025. [16] Karl HM Nyman and Peter D Ditlevsen. The middle Pleistocene transition by fre- quency locking and slow ramping of internal period. Climate Dynamics, 53(5):3023– 3038, 2019. [17] Peter Ashwin, Charles David Camp, and Anna S von der Heydt. Chaotic and non- chaotic response to quasiperiodic forcing: limits to predictability of ice ages paced by Milankovitch forcing. Dynamics and Statistics of the Climate System, 3(1), 2018. [18] Peter Ditlevsen and Peter Ashwin. Complex climate response to astronomical forc- ing: The mpt as a change in frequency locking. In Geophysical Research Abstracts, volume 21, 2019. [19] Stephen Lynch. Dynamical Systems with Applications using MATLAB®. Springer International Publishing, 2014. Appendix Nomenclature and Parameter Definitions To ensure clarity and consistency throughout this paper, we provide complete def- initions of all mathematical symbols and parameters used in our Van der Pol-Duffing oscillator analysis. Table 1 summarizes the key variables, parameters, and their physical interpretations in the context of glacial cycle dynamics. I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 19 of 20 Table 1: Nomenclature and parameter definitions for the asymmetrically forced Van der Pol-Duffing oscillator model Symbol Definition Physical Interpreta- tion Typical Values Primary Variables x Fast variable (Liénard coordi- nate) Ice volume anomaly [−3, 3] y Slow variable Rate of ice volume change [−2, 2] θ Phase variable Orbital phase [0, 1] System Parameters α Timescale separa- tion parameter Controls relaxation os- cillation sharpness 11.11 β Asymmetry pa- rameter Climate system asym- metries/thresholds [0, 1.2] γ Effective forcing amplitude Astronomical forcing strength [0.1, 3.5] κ Unforced oscilla- tion timescale Sets natural climate os- cillation period 52.0 τ Time scaling fac- tor Normalizes to glacial timescales 1.0 ω Forcing frequency Orbital frequency ratios [0.5, 5.0] GSPT Analysis S0 Critical manifold Equilibrium ice configu- rations Surface Sa Attracting sheet Stable ice states Subset of S0 Sr Repelling sheet Unstable ice states Subset of S0 θs Folded saddle lo- cation Critical orbital phases Function of a, β θn Folded node loca- tion Turning point phases Function of a, β Physical Timescales T0 Unforced period Natural glacial cycle length ∼ 100 kyr Tf Forcing period Orbital period (eccen- tricity) ∼ 100 kyr T41 Obliquity period Earth’s obliquity cycle 41 kyr T23 Precession period Precession cycle 23 kyr I alraddadi al / Eur. J. Pure Appl. Math, 18 (4) (2025), 6723 20 of 20 Computational Implementation of the First Return Map and frequency locking The first return map F is calculated using a MATLAB-based numerical integration approach following the methodology described in [19]. In this computational framework, we evaluate the component maps P+ and P− independently, whereas the jump maps J+ and J− are analytically defined as specified in equation. To approximate maps P+ and P−, we solve the dynamical system by initiating trajec- tories at y = ±2 and terminating the integration upon reaching the event where y = ±1. The implementation employs MATLAB’s ode45 differential equation solver equipped with an event detection mechanism that halts the numerical integration precisely when the so- lution trajectory intersects the boundaries y = −1 or y = 1. This occurs when trajectories evolve toward periodic attractors located on the stable manifold, preventing them from ever returning to the fold boundaries. The MATLAB implementation incorporates ex- ception handling to identify and manage these special cases through appropriate break statements. We determine frequency locking behavior through computational analysis of the Poincaré return map F iterations, where the map itself is constructed by numerically integrating the desingularized system. Our computational results, presented in Figure 7, demonstrate the iterative dynamics of the Poincaré return map for representative trajectories. To ensure convergence to attracting behavior, we discard an initial transient period of 100 iterations before analyzing the system’s period N and identifying attractor points.