EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 1, Article Number 5787 ISSN 1307-5543 – ejpam.com Published by New York Business Global The Asymmetric Periodically Forced Van Der Pol Oscillator Ibrahim Alraddadi Department of Mathematics, Faculty of Science, Islamic University of Madinah, Medina 42210, KSA Abstract. We review geometric singular perturbation theory (GSPT) which has been used to explain the behaviour of the singular slow-fast system near the singular limit. In particular, we follow the analysis of Guckenheimer et al. [10] for the periodically forced symmetric van der Pol oscillator (β = 0), then we constructed the Poincaré return map for studying the bifurcation phenomena of this model. We generalise to a asymmetric forced van der Pol oscillator for β ̸= 0. We show that the forced asymmetric van der Pol oscillator can become frequency locked due to the forcing. Then, we extend this analysis to show how the symmetry breaking parameter β in a periodically forced van der Pol oscillator influences the width of Arnold tongues (also known as frequency locking regions), and we find these frequency locking regions in the parameter space (a, ω). 2020 Mathematics Subject Classifications: 34C23, 37C25, 37-04, 37C35, 34C25, 34C26 Key Words and Phrases: Dynamical system, Periodic solutions, bifurcations Slow-fast nonlinear oscillator models can have periodic orbits with alternating slow and fast motion that are called relaxation oscillations [11]. These oscillations occur at different times in slow-fast dynamical systems. Relaxation oscillators have been used to understand a wide range of biological problems such as heartbeat (van der Mark and van der Pol model [29]), neuronal activity (the Fitz-Hugh-Nagumo model and the Morris-Lecar model [13]), and population cycles of predator-prey type [23]. Since 1920, the van der Pol oscillator was introduced to illustrate the behaviour ob- served in electrical circuits by Balthazar van der Pol and van der Mark [27]. The van der Pol oscillator is a type of a relaxation oscillator. Van der Pol and Mark investigated that the relaxation oscillation of the van der Pol oscillator is influenced by external peri- odic forcing. They found that the period of the relaxation oscillation was a proportion of the forcing period over wide parameter regimes, called “Frequency demultiplication” [28]. This phenomenon is now known as frequency locking, phase locking, entrainment or mode locking [16]. There has been a lot of significant research done on the forced van der Pol oscillator (see e.g.[5, 10, 18]). The van der Pol oscillator is applied in several fields, such DOI: https://doi.org/10.29020/nybg.ejpam.v18i1.5787 Email address: ialraddadi@iu.edu.sa (I. Alraddadi) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 2 of 23 as oscillatory processes in physics, electronics, neurobiology, and the dynamics of glacial cycles [3, 6–8, 12, 14, 20]. Forced nonlinear oscillator models can reproduce the oscillations in the climate record. Examples of two-dimensional Pleistocene ice age models are the forced van der Pol and the forced van der Pol Duffing oscillators [2]. Crucifix used in [6] a forced van der Pol oscillator (VDP) to consider astronomical forcing and asymmetry between the phases of ice forming and melting during the late Pleistocene. Additionally, the van der Pol oscillator has been modified to explain the different levels of ice volume that indicate a glacial or interglacial case by the signal of oscillates[7]. De Saedeleer et al. [7] used the van der Pol oscillator as a possible low-order model for the first time identification generalised synchronisation (see more [1, 20, 24]) between ice age cycles and astronomical forcing. The forced VDP oscillator is used by Ditlevsen and Ashwin as a conceptual model to describe the dynamics of glacial cycles and possible dynamical causes of the middle Pleistocene transition (MPT)[8]. In this paper, we use geometric singular perturbation theory (GSPT) to explain the behaviour of the singular slow-fast oscillator near the singular limit (see more [19]). In particular, we followed the analysis of Guckenheimer et al. [10] for the periodically forced symmetric van der Pol oscillator (β = 0). Section 3 constructs the Poincaré return map for studying the bifurcation phenomena of this model. In Section 4, we show how the symmetry breaking parameter β in a periodically forced van der Pol oscillator influences the width of Arnold tongues (also known as frequency locking regions). 1. The asymmetric van der Pol Oscillator We consider the modified van der Pol oscillator [2] which was proposed as a low-order model of the ice-age cycles in [6, 7]. The generalised van der Pol model has the following form [2]: τ2κ2 d2y dt2 − ατκ(1− y2) dy dt + y − γF (t) + β = 0 (1) oscillations occur even in the unforced case γ = 0. The term −ατκ(1− y2) increases the oscillations when y2 < 1 and damps the oscillations when y2 > 1. For large enough α, the van der Pol (1) has relaxation oscillations as can be seen in the system of equations. By using the Liénard transformation x = y − y3 3 − ẏ/τκα, we transform (1) into the system ẋ = 1 τκ (γk sin(2πωt)− β − y) ẏ = α τκ (y − y3 3 + x) (2) which is a slow-fast system if α ≫ 1 . The dynamics of (2) involves that the slow variable x represents the deviation of ice volume and the fast variable represents some feedback mechanism with hysteresis. In Table 1, the default parameters used for the I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 3 of 23 (a) (b) Figure 1: Phase portrait (a) and time series (b) showing periodically forced van der Pol oscillator (3) which has relaxation oscillations with slow variable x and fast variable y. The parameters values are given in Table 1. The green line is the y−nullcline, and the blue line is the x−nullcline. In (a), there is an unstable equilibrium point which is surrounded by stable limit cycle. The arrows show the direction of the vector field of (3). forced van der Pol oscillator (2) are taken from [2], ω is the frequency (obliquity forcing is ω ≈ 1 41) and k represents the amplitude of the forcing. 2. Applications to the symmetric forced van der Pol equation 2.1. The periodically forced van der Pol oscillator Note that the van der Pol oscillator for small ε has a periodic orbit with attracting and slow and fast motions. Geometric singular perturbation theory was used in [10] to understand the full system. The following discussion is based on [10, 11, 30, 31]. Here, we consider the symmetric van der Pol oscillator (2) where (β = 0) with periodic forcing [2]: ẋ = 1 τκ (γk sin(2πωt)− β − y) ẏ = α τκ (y − y3 3 + x) (3) is a non-autonomous system as shown in Figure 1. We set τκ = 1 through a suitable rescaling of time. We define a new parameter ε = 1 α , a = γ, k1 = 1 and θ = ωt of the form ẋ = a sin(2πθ)− y − β εẏ = y − y3 3 + x θ̇ = ω (4) where ẋ ≡ dx dt , ẏ ≡ dy dt and θ̇ ≡ dθ dt . When 0 < ε ≪ 1, the system has a relaxation oscillation as shown in Figure 1. The slow variable is x and the fast variable is y. We transform the slow time scale t to fast time scale T by rescaling the time t = εT I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 4 of 23 dx dT = ε(a sin(2πθ)− y − β) dy dT = y − y3 3 + x dθ dT = εω. (5) We now study the two systems in the singular limit ε = 0 as follows - namely the slow subsystem: ẋ = a sin(2πθ)− y − β 0 = y − y3 3 + x θ̇ = ω (6) and the fast subsystem: dx dT = 0 dy dT = (y − y3 3 + x− β) dθ dT = 0 (7) 2.2. The slow system The flow of (6) is called the slow flow on the critical/slow manifold. The critical manifold is in this case given by S := {(x, y, θ) ∈ R3| y − y3 3 + x = 0} (8) It is apparent that x = y3 3 − y defines the critical manifold which is the set of equilibrium points for the fast subsystem (7). Note that (4) has a repelling sheet Sr = S∩{−1 < y < 1} and two attracting sheets Sa = S ∩ {y < −1}, Sa = S ∩ {y > 1}, and trajectories of the DAEs on the critical manifold are the slow trajectories that are wholly defined up to the time when a trajectory hits a fold on the critical manifold. There are fold points on the two lines: L− = { ( − 1, 2/3, θ) : θ ∈ (0, 2πn)} and L+ = { ( 1,−2/3, θ) : θ ∈ (0, 2πn)}, where n is positive integer Z+. 2.3. The layer problem The family of differential equations (7) is called the layer problem. The solutions of differential equations dx dT = 0 and dθ dT = 0 can solved analytically to give: x(T ) = C1, θ(T ) = C2 ∀C1,2 ∈ R (9) I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 5 of 23 Figure 2: Three-dimensional view of the periodically forced van der Pol equation for β = 0. The surface represents a two-dimensional critical manifold S. A typical trajectory (blue) of (4) remains close to the stable branch (y ≤ −1) of S before jumping to other stable branch (y ≥ 1) at the lines of fold points L±. The trajectories of the fast subsystem (layer) problem (7) are computed to dy dT = g(y), where g(y) = y − y3 3 + x (10) with x acting as a parameter (see Figure 3). It is simple to find the equilibrium points of g(y) since it is dependent on the variable y. The equilibria of the equation is given by the solutions of the equation g(y) = 0 ⇐⇒ x = y3 3 − y. At the bifurcation points, the derivative of g(y) with respect to y is ∂ ∂y g(y) = 0 ⇐⇒ (1− y2) = 0. Hence, the bifurcation points (called singular points) are (x, y) = (±2 3 ,∓1). 2.4. The desingularised reduced system For ε = 0, the system (6) can be reduced to the system of ODEs on normally hyperbolic critical manifolds. The trajectories of the reduced system are good approximations to the solutions of the full system (4) near these manifolds. The projection of the system can be defined as x = φ(y, θ) on the critical manifold and the slow flow represented in the terms (y, θ). Since the critical manifold S is a curve with the function x = φ(y, θ) = y3 3 − y one can see that this curve has fold points when y = ±1. By differentiating the critical manifold to obtain ẋ = (y2 − 1)ẏ, the Implicit Function Theorem means that g(φ(y, θ), y, θ) = 0. This implies that a chain rule gives the relationship I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 6 of 23 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 1.25 x -3 -2 -1 0 1 2 3 y Figure 3: Phase Portrait showing the solutions of the fast system (7) (sold line and dashed line) and the fast flow subsystem indicated by the arrows. The blue line and dashed red line are the critical manifold S of VDP. For −1 < y < 1, there is an unstable submanifold (f ′(y) > 0) of S while the stable submanifolds (f ′(y) < 0) of S are two disjoint branches where y < 1 and y > 1. If |x| < 2 3 , then a system (7) has a single stable equilibrium point at the origin point (0, 0). In this figure, there are two stable equilibrium points (black circle) (±2 3 ,∓1) and one unstable equilibrium point (blue circle) (0, 0) for |x| > 2 3 . Note this phase portrait is independent of the slow evolution of θ. ∂g ∂y ẏ + ∂g ∂x ẋ+ ∂g ∂θ θ̇ = 0 Since gy is non zero, this equation can be solved for ẏ. gyẏ = −(gxf(x, y, θ) + gθω), where gy = 1− y2, gx = 1, gθ = 0 and f(x, y, θ) = a sin(2πθ)− y − β. Then, (1− y2)ẏ = −(a sin(2πθ)− y − β) (11) We rescale the time by t = (y2 − 1)s and substitute into the RHS of (4) to obtain the reduced system as follows: dθ ds = −ωgy, dy ds = gxf(x, y, θ) + gθω. we rewrite the system as dθ ds = ω(y2 − 1), dy ds = a sin(2πθ)− y − β. (12) which is called the desingularised reduced system of forced van der Pol equation (4). The desingularised reduced system is a time-reparameterised slow flow. The time t and the time s both move in the same direction on the stable branch of the critical manifold, but they move in opposite directions on the unstable branch of the critical manifold. In this I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 7 of 23 section we focus on the symmetric VDP (4) when β = 0, following the analysis in [10]. The equilibria of (12) are folded singularities which lie on the fold lines L±. Trajectories of (12) are located on the slow flow along stable sheets of the critical manifold until it reaches the boundary of this stable sheet on (θ,±1). These points (y = ±1) are the boundaries of two stable sheets. Then, the trajectories jump from the fold line to another stable sheet and they are crossing with y > 1 or y < −1. The equilibrium points (called folded singularities) of (12) are θ∗ = sin−1(±1 a ) 2π and y∗ = ±1. There is no solution for a < 1 as shown in Figure 5(a). We get two folded singularities on the line L±1 at (y, θ) = (±1, 1/4) for a = 1, and we have four folded singularities for a > 1 on the line L±1 as shown below p1,2(θ ∗ 1,2, y ∗) = ( sin−1( 1a) 2π , 1 ) and ( sin−1(−1 a ) 2π ,−1 ) ,where θ∗1,2 < 1/4 p3,4(θ ∗ 3,4, y ∗) = ( sin−1( 1a) 2π , 1 ) and ( sin−1(−1 a ) 2π ,−1 ) ,where θ∗3,4 > 1/4 To analyse the stability of the equilibrium points, we calculate the Jacobian matrix of (12) at (θ∗, y∗): J(θ∗, y∗) = [ 0 2ωy 2πa cos(2πθ) −1 ] and the eigenvalues of (J − λI) at (θ, y) where I is the identity matrix: λ1,2 = −1± √ 1 + 16ωπa cos(2πθ) 2 The classification of equilibria according to λ1,2 : a = 1 and θ = 1/4, J has λ1 = 0 and λ2 = −1 of which the equilibrium points are folded saddle-nodes on the line L1 or L−1. For a > 1, two equilibrium points p1,4 are folded saddles (see Figure 4). When 1 < a < √ 1 + ( 1 16πω ) 2, two other equilibrium points p2,3 are stable node (see Figure 5(b)). For a = √ 1 + ( 1 16πω ) 2, these equilibrium points p2,3 are folded nodes, and for a > √ 1 + ( 1 16πω ) 2 p2,3 are folded foci. Also, we compute stable and unstable manifolds Ws and Wu of the folded saddle (θs,±1) in the desingularised slow flow system (12) for β = 1.2 as shown in Figure 6. A trajectory Wu denotes the first intersection with y = ±1. In the backward time on the stable manifold branch, a trajectory Ws denotes the first intersection with y = ±2. In this following analysis, we compute a Poincaré map for (4) that has a return mechanism via a folded critical manifold. I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 8 of 23 (a) (b) Figure 4: Trajectories of the desingularised slow flow (12) for β = 0. (a) The parameter values are a = 2.3 and ω = 1, (b) a = 20 and ω = 5. There are four folded singularities obtained in both cases. The two folded saddles equilibria (θs,±1) are indicated by the cross and the two folded foci equilibria (θn,±1) are indicated by the circle. The stable Ws (green) and unstable Wu (blue) manifolds of the folded saddles are on the fold lines y = ±1. (a) (b) Figure 5: Trajectories of the desingularised slow flow (12) for β = 0 and ω = 0.54559. (a) shows no folded equilibria for a = 0.9. (b) Four folded singularities obtained for a = 1.00041. The two folded saddles equilibria (θs,±1) indicated by the cross and the two folded nodes equilibria (θn,±1) indicated by the circle. The stable Ws (green) and unstable Wu (blue) manifolds of the folded saddles are on the fold lines y = ±1. This case shows the connection between the folded saddle and the folded node. I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 9 of 23 (a) (b) (c) Figure 6: Trajectories of the desingularised slow flow (12) for β = 1.2. (a-b) show the four folded saddle equilibria on the fold lines y = ±1. The two folded saddles equilibria (θs,±1) are indicated by the cross and the two folded foci equilibria (θn,±1) are indicated by the circle. In (c), there are only two folded equilibria on the fold line y = −1. These folded equilibria are folded saddles and folded foci. The stable Ws (green) and unstable Wu (blue) are manifolds of the folded saddles on the fold lines y = ±1. 3. Return Map for the periodically forced singular VDP oscillator The dynamics of return map in slow-fast systems were studied by Szmolyan and Wech- selberger [26], and Guckenheimer[9–11]. In this section we study the bifurcation analysis of periodic orbits of (12) via a Poincaré return map. These bifurcations correspond to the bifurcation of periodic orbits in the periodically forced singular (VDP) oscillator (4). We follow the construction of the return map of the desingularised slow flow system as found in [10], together with our analysis of (12) in subsection (2.4). We defined the first return map from y = 2 (S2(θ, 2)) to itself for the desingularised slow flow system (12): F : S2 → S2 I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 10 of 23 by the flow of (12), S2(θ, 2) as the Poincaré section where the first intersection of the trajectory starting at S2(θ, 2) with θ being mod 1. The map F can be decomposed into four maps (discontinuous): F = J− ◦ P− ◦ J+ ◦ P+ (13) where P+ : S2 → S1, J+ : S1 → S−2, P− : S−2 → S−1, J− : S−1 → S2, Note that S±1 = L±. Two maps P±(θ, y) are well defined on the stable slow manifold S±2 (from y = ±2 to the fold lines y = ±1) and J±(θ, y) represent the transition map from fold lines L± to S∓2. In Figure 2, we note that trajectories of the periodically forced singular (VDP) oscillator (4) continue on the stable branch of the critical manifold S until it reaches a fold line L±1, then jump to another stable branch of S . This jumping behaviour is well defined by the maps: J+(θ, 1) = (θ,−2) and J−(θ,−1) = (θ, 2). In the singular case where ε → 0, there are regions where no trajectories arrive at the fold line L±1, while for sufficiently small ε, some trajectories of the nonsingular system (4) cross the fold lines near the folded saddle equilibria and continue on the unstable branch of the critical manifold S. These trajectories are called canards [15], and can appear in attractors near the folded singularity. The properties of canard trajectories in Forced van der pol systems were studied by Guckenheimer, Hoffman and Weckesser and later by Szmolyan and Wechselberger [4, 25, 26, 31]. Moreover, Guckenheimer et al. observed that there can be can be chaotic attractors for return map F (13) by analysing the behaviour associated with canard trajectories [10]. They examined various aspects of the dynamics of the return map. 3.1. Numerical approximation of the first return map The first return map F (13) is computed by using a Matlab numerical approximation scheme [17]. Following the computation technique of the first return map, we compute P+ and P− separately, while J+ and J− are explicit as known in (13). These maps P+ and P− are approximated by solving system (12): starting at y = ±2 and finishing at an event where y = ±1. we use the Matlab numerical ode45 solver with an associated termination event that stops the integration of (12) where the flow intersects y = −1 or y = 1. The map F (13) depends on three parameter values a, ω and β because they are found by solving (12). For some parameter value a (given that |a| is big enough), the map P−(θ, y) or P+(θ, y) may not be defined because it evolves towards a periodic orbit on the stable branch of the critical manifold and it can never reach the fold again. The Matlab code catches this case as special by using a break. I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 11 of 23 3.2. Dynamics and bifurcations of the first return map The return map discussed above can be iterated to understand the dynamics of tra- jectories in this singular system (4) in the limit ε → 0. These trajectories are governed by θn+1 = F (θn) for a given initial condition on S2. The intersection point of the flow (12) with Poincaré section S2 is a fixed point for the return map F (13). In particular, there will be fixed points of the Poincaré return map where F (θ) = θ. As in Figure 8, fixed points of the Poincaré return map occur where the graph of F (θ) and the diagonal line intersect point F (θ) = θ. This figure shows an example of the stability of the fixed point which is indicated by the cobwebs converging to the stable fixed point. More generally, a periodic orbit of the original ODE gives a periodic orbit of the Poincaré return map. In other words, periodic orbits of the map F correspond to periodic orbits of (4). Recall that θ will be a periodic point for F of period n if Fn(θ) = θ. Moreover, stable fixed point or stable periodic orbits for the return map that correspond to a stable periodic orbit of the flow. Figure 10 shows two examples of the first return map (13) for β = 0 and ω = 0.54559. For a < 1, the return map (13) becomes an invertible circle map (Figure 10(a)). For a > 1, the return map (13) become noninvertible (Figure 10(b)). The flow of (13) starts at S2 for various initial conditions θ. For a < 1, the graph of the first return map (13) is the smooth curve that appears in (Figure 10(a)). Whereas for a > 1, there is a gap region in the graph of F because of the presence of folded equilibria in S±1 as shwon in Figure 9. In this case the trajectories cannot leave the folded equilibria region between the folded saddle and unstable manifold of the folded saddle (see Figure 7) . This creates a discontinuity in the graph of the first return map (13) (Figure 10(b)). I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 12 of 23 (a) (b) Figure 7: Schematic diagram shows vector field of the desingularised system (12) on the upper stable part of manifold for any trajectory starts at y = 2. (a) for a < 1, there is no folded equilibria, trajectory can reach everywhere on the line y = 1. In this case the map from y = 2 to y = 1 is a homeomorphism. (b) for a > 1, A is folded saddle equilibria and B is folded sink, and C where is the unstable manifold of the folded saddle hits the fold. Note that trajectories reach everywhere on the line y = 1 except the region between A and C, this creating a gap. The reason for this gap is due to the folded equilibria region where no trajectory exits. This region is located between folded saddle equilibria (A) and the unstable manifold of the folded saddle equilibria (C). Figure 11 shows a saddle-node bifurcation of the first return map (13) for a = 1.02 and three different value of ω = 0.18, 0.19, 0.20. we see a saddle-node bifurcation of the first return map that is passing through a saddle point in (b). The saddle-node bifurcation of the periodic orbit occurs on the boundary of the frequency locking region as discussed in more details (see section 3.3). Also, the saddle-node bifurcation of the first return map (13) occurs as β varies (see Figure 13). I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 13 of 23 (a) (b) Figure 8: (a) The first return map F (13) for a = 1.1, ω = 0.35 and β = 0. The red points are the graph of the first return map and the black line is the diagonal line. The two fixed points of (13) are indicated by the intersection points of F (θ) (red) with the diagonal θ (black). The stability of these fixed points depends on the slope of the function F (θ). (b) A cobweb diagram shows a blue trajectory attracted to a stable fixed point of F (13). (a) (b) Figure 9: The first return map F (13) for β = 0. (a) corresponds to the parameter choice a = 2.3, ω = 1 in the desingularised system (12) in Figure 4(a), and (b) to a = 20, ω = 5 in Figure 4(b). Note the presence of fixed points and discontinuities. I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 14 of 23 (a) (b) Figure 10: The first return map F (13) for β = 0 and ω = 0.54559. (a) Corresponds to a = 0.9, the parameter choice in the desingularised system (12) in Figure 5(a), and (b) to a = 1.00041, the parameter choice in Figure 5(b). The first return map (13) is an invertible circle map in (a). (b) The first return map (13) becomes noninvertible. Note the gap appearing in (b). The reason for this gap is due to the folded equilibria region where no trajectory exits. This region is located between the folded saddle and the unstable manifold of the folded saddle. (a) (b) (c) Figure 11: The first return map F (13) for ω = 0.35 with different values of a (a-c) where a = 1.05, a = 1.08 and a = 1.1 respectively. There is a saddle-node bifurcation of the first return map at a = 1.08 where the red curve of F becomes tangent to the diagonal line (black). 3.3. Frequency locking in the iteration of the return map Frequency locking can occur for iterations of the Poincaré map and there are transitions between periodic (also called phase locking) and quasiperiodic motion [22]. Both types of behaviour can be distinguished by computing the rotation number N [20]. If we find that the rotation number N of the Poincaré return map is rational, i.e. N = p/q with p and q integers, then we say the motion is phase-locked. Otherwise, we say the motion is I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 15 of 23 quasiperiodic. Frequency locking regions (also called Arnold tongues) identify a periodic solution which generically persists as the amplitude a and the frequency ω of the forcing are varied [22]. Boundaries of Arnold tongues correspond to saddle-node bifurcations of these N -periodic orbits[21]. In this subsection, we show the frequency locking regions (or Arnold tongues) in the parameter space (a, ω). Subsequently, we show that a periodic solution persists for varying a and ω and that this solution is frequency locked (forming Arnold tongue). As shown in Figure 11, the saddle-node bifurcation of the first return map occurs on the boundary of frequency locking in a parameter space (a, ω) of the system for relatively weak forcing. In parameter space (a, ω), frequency locking regions exist where the period N of the oscillator remains constant under perturbation. We identify the frequency locking by estimating the iterations of the Poincaré return map F and we numerically solved the desingularised system (4) to construct the return map F . From our computation, shown in Figure 12, we show the iterations of the Poincaré return map for a typical point after ignoring a transient of 100 iterations so as to identify the period N of the oscillator to get an attractor point of the system. We describe our method as the following: (i) For different values of a and ω, we numerically solved the desingularised system (4) to construct the Poincaré return map F . (ii) Pick a random initial condition θ∗. (iii) We set θ0 = Fnt(θ∗) mod 1, where nt iterates 500 times for the Poincaré return map F (13). (iv) We use the iterations of the Poincaré return map for a typical point after ignoring a transient of 100 iterations so as to identify a trajectory attractor of the system. (v) We try to find the average period of θi where i ≤ nmax by searching for smallest i ∈ {1, . . . , nmax} such that |(F i(θ0) mod 1)− θ0| < ϵ (a) If this is satisfied, then we set N = F i(θ0)− θ0. (b) If it is not satisfied, then we set N = 0. 4. Symmetry breaking In Sections 2 and 3 we reviewed the results of [10], the symmetric forced van der Pol oscillator (4) for β = 0. Here we generalise to a asymmetric forced van der Pol oscillator for β ̸= 0. In this section, we show that the forced asymmetric van der Pol oscillator can become frequency locked due to the forcing. We also show how the symmetry breaking parameter β in a periodically forced van der Pol oscillator influences the width of Arnold I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 16 of 23 Figure 12: (a, ω) parameter plane showing regions of the existence of a stable period of N orbits for the return map F in (13) for the symmetric van der Pol oscillator. The colorbar indicates the rotation numbers N of the return map. Observe in this case β = 0, (from down to up) there are phase-locking regions in N -periodic orbits of (13) (in differ- ent colours). These regions correspond well with standard Arnold tongues type pictures (frequency-locking): the frequency ratio is ω and the forcing amplitude is a. The large regions correspond to Arnold tongues for frequency-locking 1:1 (black), 3:1(light blue), 5:1(pink), 7:1(dark blue). Note that there are bistability regions between these tongues depending on the initial conditions. The boundary of the Arnold tongues is typically the saddle-node bifurcation of a periodic point of map F (13) (compare to [10, Figure 7.1]). tongues (also known as frequency locking regions), and we find these frequency locking regions in the parameter space (a, ω). (a) (b) (c) Figure 13: The first return map F (13) when a = 1.02 and ω = 1.19 with different values of β. (a-c) show β = 0.8, 0.81 and 0.82 respectively, a saddle-node bifurcation of the first return map F (13). Figure 14 illustrates the frequency-locking regions of the iteration of F (13) in the I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 17 of 23 two-parameter plane (a, ω) for different values of β = 0.25, 0.5, 0.8, 1. The Arnold tongue regions appear in this case β ̸= 0. Note that the width of Arnold tongues depends on the parameter β. There are large Arnold tongue regions of period 1, 3, 5, 8 for small β = 0.5 (Figure 14(a)). Increasing β to 1, leads to slightly smaller Arnold tongue regions of periods between 1 and 25 as shown in Figure 14(b-d). Moreover, (Figure 14(d)) shows the bistability regions, but these regions are very narrow as compared with β = 0 in Figure 12. There is a reason for the white wedge, where no rotation numbers are shown. In other words, trajectories start to get non returning points in the white wedge. (a) (b) (c) (d) Figure 14: The (a, ω) parameter plane showing regions of the stable periodic orbit of the return map F (13) as β ̸= 0. The colorbar indicates the rotation period N of (13) for β = 0.5, 0.8, 1 and 1.2. The Arnold tongues regions appear for β ̸= 0. (a) The large regions correspond to Arnold tongues for frequency-locking 1:1 (brown), 3:1 (purple), 5:1 (light blue) and 8:1 (brown). As β increases, (b-d) there are smaller Arnold tongues regions for frequency-locking between periods 1 and 15 as a region where the return is no longer defined appears for small a. (d) there are bistability regions but these regions are very narrow as compared with β = 0 in Figure 12. Figure 15 shows the graph of the first return map (13) for β = 1 and a with different I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 18 of 23 values of ω. The saddle-node bifurcation of the first return map (13) occurs as ω varies. In Figure 14(c), this case corresponds to the attractor point of periodic orbits moving forward from inside a frequency locking region to the outside of the frequency locking region. As ω increases in Figure 16(a), we note this discontinuity in the first return map F (13) hits the periodic orbit. In Figure 17, the black circle corresponds to the value of (a, ω, β) in Figure 16. (a) (b) (c) Figure 15: The graph of the first return map F (13) for β = 1 and a = 1.65 with different values of ω. (a-c) shows a saddle-node bifurcation of the first return map for β = 1 as ω presents throughout ω = 1.81. I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 19 of 23 (a) (b) Figure 16: (a) The graph of the first return map F (13) for β = 1 and a = 1.65 and ω = 1.92. (b) The flow of (13) on the desingularised slow flow system (12) for two different initial values of θ. The value of θ is 0.271 (blue) and 0.273 (red). For (a) θ = 0.271 is on the left curve of the diagonal line, and θ = 0.273 is on the right curve. For (b) The blue and red trajectories start close together, hit the upper fold line, and then jump to another branch. The red trajectory goes to the left of the folded saddle in the lower folded line. While the blue trajectory comes close to the folded saddle. Vertical lines indicate the endpoint of the red and blue trajectories. (b) corresponds to the black circle in Figure 17. (b) For periodic orbits at discontinuity where ε → 0. While periodic orbits at the discontinuity for sufficiently small ε > 0, some layer trajectories of the nonsingular system (4) cross the fold lines near the folded saddle equilibria and continue the unstable branch of the critical manifold S then move to another stable branch of the critical manifold S. Periodic orbits without discontinuity where ε ̸= 0, there is a region of canard trajectory where no trajectories jump at the fold line. I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 20 of 23 parameter value description α 11.11 Fast/slow timescale separation β 0.25 Symmetry breaking γ 0.75 Effective forcing amplitude κ 35.09 Unforced oscillation timescale has period 100 kyr for time scaling τ = 1 τ 1 Time scaling Table 1: shows the default parameters for the forced van der Pol oscillators (2) (taken from [2]). (a) (b) Figure 17: The (a, ω) parameter plane showing regions of the stable periodic orbit of the return map F (13) for β = 1 as shown in Figure 14(c). (b) shows a small range region for a and ω. Note the white region is not identified as a low period periodic orbit of the system as shown between frequency-locking regions 4:1 (yellow) and 5:1 (pink). The black circle represents a = 1.65 and ω = 1.92 as shown in Figure 16. 5. Conclusion We use geometric singular perturbation theory (GSPT) to explain the behaviour of the singular slow-fast system near the singular limit. In particular, we followed the analysis of Guckenheimer et al. [10] for the periodically forced symmetric van der Pol oscillator (β = 0), then we construct a Poincaré return map to study the bifurcation phenomena of this model. Subsequently, we show that a periodic solution persists for varying a and ω and that this solution is frequency locked (forming Arnold tongues). As shown in Figure 11, we show that the saddle-node bifurcation of the first return map occurs on the boundary of frequency locking in a parameter space (a, ω) of the system for relatively weak forcing. In parameter space (a, ω), frequency locking regions exist where the period N of the oscillator remains constant under perturbation. In this thesis, we identify the frequency locking by estimating the iterations of the return map F and we numerically solve the desingularised system (12) to construct the return map F . From our computations such as shown in Figure 12, we can identify the rotation numbers of attractors of the system. For β ̸= 0, we I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 21 of 23 show that the forced asymmetric van der Pol oscillator can become frequency locked due to the forcing. We also show how the symmetry breaking parameter β in a periodically forced van der Pol oscillator influences the width of Arnold tongues (also known as frequency locking regions), and we showed the frequency locking regions (or Arnold tongues) in the parameter space (a, ω). Data availability statement: All data generated or analysed during the research investigation have been included in the paper. References [1] Henry DI Abarbanel, Nikolai F Rulkov, and Mikhail M Sushchik. Generalized syn- chronization of chaos: The auxiliary system approach. Physical review E, 53(5):4528, 1996. [2] 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. [3] Alexander Balanov, Natalia Janson, Dmitry Postnov, and Olga Sosnovtseva. From simple to complex. Springer. [4] Katherine Bold, Chantal Edwards, John Guckenheimer, Sabyasachi Guharay, Kath- leen Hoffman, Judith Hubbard, Ricardo Oliva, and Warren Weckesser. The forced van der Pol equation II: Canards in the reduced system. SIAM Journal on Applied Dynamical Systems, 2(4):570–608, 2003. [5] Mary L. Cartwright and JE Littlewood. On non-linear differential equations of the second order. J. London Math. Soc, 20, 1945. [6] Michel Crucifix. Oscillators and relaxation phenomena in Pleistocene climate theory. Trans. R. Soc. A, 370:1140–1165, 2012. [7] 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. [8] 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. [9] John Guckenheimer. Return maps of folded nodes and folded saddle-nodes. Chaos: An Interdisciplinary Journal of Nonlinear Science, 18(1):015108, 2008. [10] 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. [11] John Guckenheimer, Kathleen Hoffman, and Warren Weckesser. Bifurcations of Re- laxation Oscillations Near Folded Saddles. I. J. Bifurcation and Chaos, 15(11):3411– 3421, 2005. I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 22 of 23 [12] Chihiro Hayashi. Nonlinear oscillations in physical systems. Princeton University Press, 2014. [13] Sze-Bi Hsu and Junping Shi. Relaxation oscillation profile of limit cycle in predator- prey system. Discrete & Continuous Dynamical Systems-B, 11(4):893, 2009. [14] Eugene M. Izhikevich. Dynamical Systems in Neuroscience: The Geometry of Ex- citability and Bursting. The MIT Press, 07 2006. [15] Christian Kuehn. Multiple time scale dynamics, volume 191. Springer, 2015. [16] Yuri A. Kuznetsov. Elements of Applied Bifurcation Theory, Second Edition. Library, page 591, 1998. [17] Stephen Lynch. Dynamical Systems with Applications using MATLAB®. Springer International Publishing, 2014. [18] R Mettin, U Parlitz, and W Lauterborn. Bifurcation structure of the driven van der Pol oscillator. International Journal of Bifurcation and Chaos, 3(06):1529–1555, 1993. [19] Karl H M Nyman, Peter Ashwin, and Peter D Ditlevsen. Bifurcation of critical sets and relaxation oscillations in singular fast-slow systems. Nonlinearity, 33(6):2853– 2904, apr 2020. [20] Arkady Pikovsky, Jurgen Kurths, Michael Rosenblum, and Jürgen Kurths. Synchro- nization: a universal concept in nonlinear sciences, volume 12. Cambridge university press, 2003. [21] Arkady Pikovsky, Michael Rosenblum, and Jürgen Kurths. Synchronization: A Uni- versal Concept in Nonlinear Sciences. Cambridge Nonlinear Science Series 12, 2003. [22] Arkady S Pikovsky, Ulrike Feudel, and Sergey P Kuznetsov. Strange nonchaotic attractors: Dynamics between order and chaos in quasiperiodically forced systems, volume 56. World Scientific, 2006. [23] C. Rocsoreanu, A. Georgescu, and N. Giurgiteanu. The FitzHugh-Nagumo Model: Bi- furcation and Dynamics. Mathematical Modelling: Theory and Applications. Springer Netherlands, 2012. [24] Nikolai F Rulkov, Mikhail M Sushchik, Lev S Tsimring, and Henry DI Abarbanel. Generalized synchronization of chaos in directionally coupled chaotic systems. Phys- ical Review E, 51(2):980, 1995. [25] Peter Szmolyan and Martin Wechselberger. Canards in R3. Journal of Differential Equations, 177(2):419–453, 2001. [26] Peter Szmolyan and Martin Wechselberger. Relaxation oscillations in R3. Journal of Differential Equations, 200(1):69–104, 2004. [27] Balth Van der Pol. On “relaxation-oscillations”. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11):978–992, 1926. [28] Balth Van der Pol and Jan Van Der Mark. Frequency demultiplication. Nature, 120(3019):363–364, 1927. [29] Balth Van Der Pol and Jan Van Der Mark. LXXII. The heartbeat considered as a relaxation oscillation, and an electrical model of the heart. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 6(38):763–775, 1928. [30] Ferdinand Verhulst and Taoufik Bakri. The dynamics of slow manifolds. Journal of I. Alraddadi / Eur. J. Pure Appl. Math, 18 (1) (2025), 5787 23 of 23 the Indonesian Mathematical Society, 13:1–16, 2007. [31] Martin Wechselberger. À propos de canards (Apropos canards). Transactions of the American Mathematical Society, 364(6):3289–3309, 2012.