EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 4, 2024, 3708-3729 ISSN 1307-5543 – ejpam.com Published by New York Business Global Periodic Solutions and Bifurcations in a Long-Timescale Oceanic Carbon Cycle Model Ibrahim Alraddadi 1 Department of Mathematics, Faculty of Science, Islamic University of Madinah, Medina 42210, KSA Abstract. We study the dynamics of a recent model of Rothman for long timescale carbon cycle. We reproduce and extend various results of the Rothman model. We present numerical results showing that the model exhibits both stable and unstable limit cycles via Hopf bifurcations as the parameters are varied. We numerically find normal forms of Bautin bifurcations to confirm their criticality. We also extend the analysis of the normal form coefficients to identify where the fold limit cycle bifurcation occurs 2020 Mathematics Subject Classifications: 34C23, 37C25, 37-04, 37C35,3 4C25, 34C26 Key Words and Phrases: Dynamical system, Periodic solutions, bifurcations 1. Introduction A modern view of climate variability represents a qualitative change in the dynamics of the Earth system as shown in [15]. On short time scales, the Earth’s carbon cycles are between photosynthesis, which converts carbon dioxide (CO2) to organic carbon, and respiration, which converts organic carbon back to CO2. Human activities have caused changes in atmospheric CO2 levels that have recently increased by 50% compared to the previous years in ages past, most of which has been absorbed by the oceans [8]. There is a lot of interest in forecasting how the climate will respond to these changes, as shown in [15]. Over long timescales, other reasons of carbon can vary, in particular in the ocean. Significant disruptions from the global carbon cycle are caused by the changes in the concentration of carbon in the ocean. Rothman [12] studied the evolution of the carbonate system in the upper ocean as part of the marine carbon cycle. These changes affect the marine carbonate system, which is a major element of fundamental importance to the Earth’s carbon cycle. Environmental perturbations can significantly result in considerable disruptions in the Earth’s carbon cycle. To understand these observations, Rothman [12] suggested a mathematical model to describe the characteristic rate of change in the marine carbon cycle. This model formed DOI: https://doi.org/10.29020/nybg.ejpam.v17i4.5458 Email address: ialraddadi@iu.edu.sa (I. Alraddadi) https://www.ejpam.com 3708 Copyright: © 2024 The Author(s). (CC BY-NC 4.0) I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3709 as a two-equation dynamical system to describe the dynamics of carbon and the response feedback mechanisms of shallow-oceans in the marine carbon cycle. Rothman’s carbon cycle model builds upon a foundation of dynamical systems approaches in climate science. Recent work by Crucifix [1] on oscillators and relaxation phenomena in Pleistocene cli- mate theory provides context for understanding the significance of Rothman’s model in capturing complex carbon cycle dynamics. Additionally, Lenton et al. [5] introduced the concept of tipping elements in the Earth system, many of which involve the carbon cycle, underscoring the importance of identifying critical thresholds and bifurcations in climate models. The theory of dynamical systems provides many useful techniques to analyse the observed behaviour of the model, such as bifurcation phenomena. This approach can help us to understand the important features of the marine carbonate system. In this article presents a study of oscillations in a long-timescale carbon cycle model (7) introduced by Rothman [12]. We briefly introduce the carbon cycle model, analyse the stability of the equilibrium point and show the region of stability in parameter space. We present numerical results showing that the model exhibits both stable and unstable limit cycles via Hopf bifurcations as the parameters are varied. We study the normal form of Hopf and Bautin bifurcations for Rothman’s carbon cycle model (7). We review the technique for computing Lyapunov coefficients. We use MAPLE to compute the first and second Lyapunov coefficients l1 and l2 in the normal form at a bifurcation point where l0 = 0. We also extend the analysis of the normal form coefficients to identify where the fold limit cycle bifurcation occurs. 2. Rothman’s Carbon Cycle Model A model of the marine carbon cycle illustrates the interactions of shallow-ocean res- piration with fluxes of carbon into and out of the shallow ocean. We consider the upper ocean to be a well-mixed open system. The main carbonate input is dissolved calcium carbonate (CaCO3) which is transported to the shallow ocean by rivers and the output from the shallow ocean to the sediments. Rothman’s carbon cycle model builds upon his previous studies on carbon cycle dynamics [10, 11, 13]. The system can be written in the following mathematical form [12]: ȧ = 2[jin −B(a,w)], ẇ = (1 + ν)jin −B(a,w) +R(a,w)− (w − w0)/τw (1) where a is total alkalinity, w is total dissolved inorganic carbon, and other terms represent various fluxes and parameters. Key advancements in this model (1) include the explicit representation of ocean carbonate chemistry, allowing for a more detailed analysis of ma- rine carbon dynamics. The incorporation of nonlinear feedbacks captures complex interac- tions within the carbonate system, enabling the model to produce rich dynamical behavior. While Rothman’s earlier works [10, 13] primarily addressed either long-term (million-year) carbon cycle stability or conceptual frameworks, and other models in the field often fo- cused on short-term (annual to decadal) fluctuations, this model uniquely bridges the gap I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3710 by exploring carbon cycle variations on intermediate timescales, ranging from thousands to hundreds of thousands of years. Notably, the model demonstrates the potential for self- sustained oscillations, providing a mechanism for intrinsic carbon cycle variability that does not depend on external forcing. While Saltzman and Maasch’s low-order models [7] focus on glacial-interglacial cycles, Rothman’s approach specifically addresses the dynam- ics of the upper ocean carbonate system on intermediate timescales, providing a more detailed treatment of carbonate chemistry feedbacks. Paillard and Parrenin’s conceptual model [9] relies heavily on threshold mechanisms; in contrast, Rothman’s model demon- strates how continuous nonlinear feedbacks can lead to complex dynamics, including limit cycles. This model’s significance lies in its ability to bridge short-term variability and long- term trends in the carbon cycle, offering new insights into intrinsic carbon cycle variability and potential explanations for observed historical variations. The model (1) expresses the concentrations of total alkalinity (a) and total dissolved inorganic carbon (DIC) w in units of µmot.kg−1. Figure 1 illustrates a schematic diagram of (1), jin is the rate equivalent to the change of dissolved CaCO3 concentration added by the rivers in oceans. The flux B(a,w) indicates the rate at which CaCO3 exits from the system to be buried in sediments in the deep ocean. The flux R(a,w) shows the rate of the upper ocean production of CO2 by the respiration process. The dimensionless term ν corresponds to the strength of an external source of CO2 (e.g. volcanic emissions). The functional specification of feedback (B and R) is expressed below B(a,w) = B[c(a,w)] = b ∗ jins(c, cp), s(c, cp) = cγ cγ + cγp (2) where s is a sigmodial function, and bjin denotes the maximum rate of carbonate burial B. For s = 0.5, cp represents the transitional carbonate ion concentration. The exponent γ is the sigmoid sharpness of the transition R(a,w) = R[c(a,w)] = θjins̄(c, cx), (3) where s̄ = 1− s, θjin is the maximum value of carbonate respiration R, and cx represents the crossover concentration at the midway point that exists when R transits from a high state to a low state. The differential equation in (4) below describes the rate of change of c which depends on a and w dc dt = da dt ∂c ∂a + dw dt ∂c ∂w (4) To approximate the magnitude of the partial derivatives ∂c ∂a and ∂c ∂w , we use the ”buffer function” f(c) [12](see Figure 2): ∂c ∂a ≈ − ∂c ∂w ≡ f(c), f(c) = f0 cβ cβ + cβf (5) where f0, β and cβf are specific values for the equilibrium of the carbonate system which reduces equation (5) to the expression in (6): ċ = f(c)(ȧ− ẇ) (6) I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3711 Figure 1: Schematic diagram of the evolution of the model (1) as shown in [12]. The left and right panels show how total alkalinity a and total dissolved inorganic carbon (DIC) w changed over the chemical state of the upper ocean. This diagram has three blue lines: the wavy line shows the air and water interface, and a thick line separates the shallow sea from the deep sea. The dashed line represents the sediment–seawater interface. Concentration fluxes are shown by arrows that only move in one direction. The dashed line between the sediment and seawater indicates no dynamic difference between the deep sea and sediment. This feedback B(a,w) and R(a,w) determines the extent of a and w accumulation before leaving the system. where c is the concentrations of CO2− 3 . By substituting (2) and (3) into model (1), we obtain ċ as a function of c and w: ċ/f(c) = µ[1− bs(c, cp)− θs̄(c, cx)− ν] + w − w0 ẇ = µ[1− bs(c, cp) + θs̄(c, cx) + ν]− w + w0 (7) is called Rothman’s carbon cycle model [12] with respect to the new timescale dT = dt/τw. The value of parameters µ, b, θ, cp, cx, cf , w0, γ and β are listed in Table 1. These parameters correspond to the properties of the modern ocean. We change the value of cx, b, µ, cp, θ and ν to show the stability region of (7) in the parameter space. 3. Steady state and stability analysis In this section, we analyse the behaviour of Rothman’s carbon cycle model (7) and in particular, the stability of the equilibrium point. The nullclines (ċ = ẇ = 0) of (7) are given by w0 = µ[1− bs(c, cp)− θs̄(c, cx)− ν] + w or f(c) =0 w = µ[1− bs(c, cp) + θs̄(c, cx) + ν] + w0 (8) I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3712 Figure 2: The buffer function f(c) characterises the equilibrium chemistry of the carbonate system. The smooth curve is approximated by (5). The parameters f0 = 0.694, β = 1.7 and cf = 43.9 are specific values for the equilibrium of the carbonate system. where s̄(c, cx) = 1 − s(c, cx). To find the equilibrium point (c∗, w∗) of (7), we substitute w into (ċ = 0): µ[2− 2bs(c, cp)] = 0 bs(c, cp) = 1 b cγ cγ + cγp = 1 c∗ = (b− 1) −1 γ cγp (9) from (3), w∗ = µ[θ(1− s(c, cx)) + ν] + w0 w∗ = w0 + µ[θ(1− cγ cγ + cγx ) + ν] w∗ = w0 + µ[θ − θcγp cγp + (b− 1)cγx + ν] (10) Now, we introduce new variables v and u for shifting an equilibrium point (c∗, w∗) to the origin (0, 0): v = c− c∗ u = w − w∗ is an equilibrium point for a two-dimensional dynamical system: v̇ = f(v, u) u̇ = g(v, u) (11) where g and f are nonlinear, corresponding to the linear form[ v̇ u̇ ] = [ fv fu gv gu ] [ v u ] (12) where J(v, u) = [ fv fu gv gu ] I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3713 Figure 3: The phase portrait (a) and time series (b-c) show a stable limit cycle of Roth- man’s carbon cycle model (7) for cx = 70. The initial condition and other parameters are given values in Table 1. The black circle represents an unstable equilibrium point (c∗, w∗). The nullclines represent the red line ẇ = 0 and black line ċ = 0. is the Jacobian matrix of (11) at the equilibrium point (v, u). The trace τ(v, u) and the determinant ∆(v, u) of J are τ = trJ = fv + gu ∆ = detJ = fvgu − gvgu To analyse the stability of the equilibrium point, we find the eigenvalues of (J − λI) at (u, v) where I is the identity matrix: λ1,2 = τ ± √ τ2 − 4∆ 2 according to τ and ∆ of J , the stable equilibrium point is (τ < 0) and unstable equilibrium point is (τ > 0). Now, we return to the Rothman carbon cycle model (7), τ = −1± ∆ 2 (1− θbcγpc γ x cγp + [(b− 1)cγx]2 ) (13) I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3714 Figure 4: This figure shows the stability region of the equilibrium point (c∗, w∗) for Roth- man’s carbon cycle model (7) in parameter space of cp and b. Other parameters are given in Table 1. The stability region is divided into four regions. The yellow and green regions represent the stable and unstable node respectively, where (c∗, w∗) has real eigenvalues. The blue and red regions represent a stable spiral and an unstable spiral where (c∗, w∗) has complex eigenvalues. Note Hopf bifurcations occur on the boundary between the blue and red regions where (c∗, w∗) has purely imaginary eigenvalues. where, ∆ = 2µγf0(b− 1) 1 γ +1 cβ−1 p b(b− 1)β/γcβf + bcβp (14) since b > 1 and all parameters are positive, the sign of τ depends on θ and µ (i.e. τ > 0 when µ ≫ 0 and τ < 0 when θ > 0). λ1,2 are pure imaginary when τ = 0 (let µ = µc where τ = 0). The stability of this equilibrium point (c∗, w∗) is determined by the sign of the eigenvalues λ1,2 of the Jacobian Matrix. We use (13) and (14) to find numerically (using MATLAB) λ1,2. In Figure 5, we examine the stability of (c∗, w∗) over a range of values of parameters b (a-b) and cp (c-d), respectively. Other parameter values of (7) are given in Table 1. The colored branch shows the location and the stability of (c∗, w∗) when one parameter is varied. Green/yellow branches indicate that a stable/unstable node equilibrium point has positive /negative real eigenvalues, respectively. Red/blue branches indicate that a stable/unstable spiral equilibrium point has negative/positive real parts of complex eigenvalues, respectively. As discussed in [12], Rothman identifies supercritical and subcritical Hopf bifurcation branches and the stable equilibrium point are in the parameter space. Here, we review some results from [12] to extend in the analysis of the equilibrium point in the parameter spaces as shown in Figure 6. The stability of (c∗, w∗) is indicated by the yellow region (stable node), the green region (unstable node), the blue region (stable spiral) and the red region (unstable node). A Hopf bifurcation occurs in the boundary between the red and blue regions where the eigenvalues have pure imaginary values as shown in [12, Fig.3]. I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3715 Figure 5: The location and stability of the equilibrium point are shown for Rothman’s carbon cycle model (7) with parameters as given in Table 1. Panels (a) and (b) show the location and the stability of the equilibrium point (c∗, w∗) for varying the parameter b. Similarly, panels (c) and (d) show the location and stability, but for the parameter cp instead. In all cases, other parameter values are given in Table 1. The stability of the equilibrium point (c∗, w∗) is determined by the sign of the eigenvalues λ1,2. The location and the stability of (c∗, w∗) is indicated by blue (stable spiral) and red (unstable spiral). The analysis and explanation of Rothman’s carbon cycle model is discussed in more detail in [12]. The model (7) exhibits both stable and unstable limit cycles via Hopf bifurcations as the parameter cx is varied. There is only one stable equilibrium point (c∗, w∗) of (7) for cx < 55.89 (See Figure 7). An unstable limit cycle appears between the stable equilibrium point and the stable limit cycle when 55.89 < cx < 62.61 (See Figure 8). There is a stable limit cycle and an unstable equilibrium point when cx > 62.61 (See Figure 3). 4. Bifurcations of Rothman’s carbon cycle model This section illustrates the normal form of Hopf and Bautin bifurcations, which include the saddle-node bifurcation of limit cycles. Section 4.1 illustrates the conditions and highlights the Lyapunov coefficients that appear in this normal form. In addition, we review the technique for computing Lyapunov coefficients. We use Theorem 1 in Section 4.1 and Theorem 2 in Section 4.2 to identify the topological normal form for Hopf and I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3716 Figure 6: The stability region of the equilibrium point (c∗, w∗) is shown for Rothman’s carbon cycle model (7). In the parameter space (a-d), two parameters are varying and other parameters are given in Table 1. (a-c) show the stability diagram for the parameter spaces of µ against b, cx and θ, (d) shows the parameter space of θ and cx. The stability region is divided into four regions: yellow (stable node) and green (unstable node) regions, which indicate where (c∗, w∗) has real eigenvalues. Stable and unstable spiral equilibrium points are indicated by blue and red regions where (c∗, w∗) has complex eigenvalues. Note Hopf bifurcations occur at the boundary between blue and red regions where (c∗, w∗) has purely imaginary eigenvalues. Bautin bifurcations. Section 5 derives the normal form of Hopf and Bautin bifurcations for Rothman’s carbon cycle model (7). In this section, we identify the Bautin bifurcation and classify its criticality by finding appropriate coefficients in the normal form. In addition, we use the normal form coefficients to identify the saddle-node bifurcation of limit cycles when it occurs. 4.1. Computation of the critical normal form coefficients The coefficients of the normal form at the bifurcation point can be calculated in several ways. We follow a method that is used in [4] and then we apply it to an example of ODEs dependent on one parameter. We consider the autonomous system ẋ = f(x, α) (15) I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3717 Figure 7: The phase portraits show solutions of Rothman’s carbon cycle model (7) for ν = 0.4. The parameter cx = 55.88 (a), cx = 60 (b) and cx = 62.62 (c), and other parameters are given in Table 1. (b-c) show solutions of (7) with different initial conditions. The model (7) exhibits both stable and unstable limit cycles via Hopf bifurcations as the parameter cx is varied. In (a-b), there is a stable equilibrium point (black point) (c∗, w∗). (b) An unstable limit cycle appears between the stable equilibrium point and the stable limit cycle. In (c), there is a stable limit cycle and an unstable equilibrium point (black circle) when cx = 62.62. with smooth right hand side f and the vector of state variables x ∈ Rn and parameters α ∈ Rm, having at (x, α) = (0, 0) an equilibrium point with eigenvalues λ1,2 = ±iω0, ω0 > 0. Consider the system (15) can be written as ẋ = A(α)x+ F (x, α) (16) where A(α) is the Jacobian matrix of (15) and F is a smooth vector field with zero Jacobian, so the Taylor expansion in x starts with the quadratic term as follows: F (x, 0) = 1 2! B(X,X) + 1 3! C(X,X,X) + 1 4! D(X,X,X,X) +O(||X||)5 (17) I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3718 Figure 8: The phase portrait (a) and time series (b-c) are shown for Rothman’s carbon cycle model (7) when cx = 55.88. The initial conditions and parameters are given in Table 1. The black point represents the stable equilibrium point (c∗, w∗) and the nullclines are shown by the red line ẇ = 0 and black line ċ = 0. where B,C and D are multilinear symmetric functions [4]: B(X,Y ) = n∑ j,k=1 ∂2Fi(ξ) ∂ξj∂ξk ∣∣∣∣ ξ=0 XjYk, i = 1, . . . , n (18) C(X,Y, Z) = n∑ j,k,l=1 ∂3Fi(ξ) ∂ξj∂ξk∂l ∣∣∣∣ ξ=0 XjYkZl, i = 1, . . . , n (19) C(X,Y, Z, U) = n∑ j,k,l,r=1 ∂4Fi(ξ) ∂ξj∂ξk∂ξl∂ξr ∣∣∣∣ ξ=0 XjYkZlUr, i = 1, . . . , n (20) The Jacobian matrix A of (15) at the equilibrium point (x, α) = (0, 0) has pure imaginary eigenvalues λ1,2 = ±iω0 (the Hopf bifurcation condition is satisfied). I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3719 parameter value description µ 250 characteristic concentration jinτw b 4 maximum CaCO3 burial rate θ 5 maximum respiration feedback rate cx 70 respiration cp 110 burial ν 0 injection rate w0 2000 reference DIC concentration γ 4 sigmoid sharpness index f0 0.694 maximum buffer factor cf 43.9 buffering β 1.7 sigmoid sharpness index Table 1: This table shows the default parameters for Rothman’s carbon cycle model (7) [12, SI Appendix, Table S1]. Let q(α) ∈ C2 be an eigenvector of A corresponding to λ(α), and p(α) ∈ C2 be an eigenvector of the transpose AT corresponding to λ(α). It follows that: A(α)q(α) = λ(α)q(α), AT (α)p(α) = λ(α)q(α) ⟨p, q⟩ = 2∑ i=1 p̄iqi = 1 We introduce the complex variable z = (p, x), where x = zq + z̄q̄ System (15) can be written as the formula [4, Lemma 3.3, p.66]: ż = λ(α)z + g(z, z, α), (21) where g(z, z̄, α) = ⟨p, F (zq + z̄q̄, α)⟩ is a smooth function and can be written as a formal Taylor series in two complex variables (z and z̄): g(z, z̄, α) = ∑ k+l≥2 1 k!l! gkl(α)z kz̄l, k, l = 0, 1, ... (22) where gkl(α) = ∂k+l ∂zk∂z̄l ⟨p, F (zq+ z̄q̄, α)⟩ ∣∣ z=0 . To simplify (15), we remove all the quadratic terms by coordinate changes as shown by the following [[4],Lemma(3.6) p.96]: ż = λz + g30 6 z3 + g21 2 z2z̄ + g12 2 zz̄2 + g03 6 z̄3 +O(|z|)4 (23) where λ = µ(α) ± iω(α), µ(0) = 0 ,ω(0) = ω0 > 0 and gij = gij(α), can be transformed by an invertible parameter-dependent change of complex coordinates, I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3720 z = w + h20 2 w2 + h11ww̄ + h02 2 w̄2 + h03 6 w3 + h12 2 ww̄2 + h03 6 w̄3 h20 = g20 λ , h11 = g11 λ , h02 = g02 2λ̄− λ , h03 = g03 2λ̄− λ h12 = g12 2λ̄ (24) for all sufficiently small |α|, for an equation with only one cubic term: ẇ = λw + c1w 2w̄ +O(|w|)4 (25) where c1 = g20g11(2λ+λ̄) 2|λ|2 + |g11|2 λ + |g02|2 2(2λ−λ̄) + g21 2 . From [4, Lemma 3.7, p 98], suppose µ(0) ̸= 0 and Re(c1(0)) ̸= 0. Then, the equation can be transformed by a parameter-dependent linear coordinate transformation, a time rescaling, and a nonlinear time reparametrization into an equation of the form du dθ = (β + i)u+ l1(β) |l1(β)| u|u|2 +O(|u|4) = (β + i)u+ su|u|2 +O(|u|4) (26) where u is a new complex coordinate, and θ is the time and β is a new parameter, and s = sign l1(0) = sign Re(c1(0)) = ±1. The real function l1(β) is called the first Lyapunov coefficient [4, Definition 3.3,p 99]. To describe bifurcation analytically, we present an example of the first Lyapunov coefficient l1 in the complex form of the dynamical system near the original equilibrium point. Example 1. Consider the two differential equations which depend on one parameter: ẋ =αx− y − x(x2 + y2) ẏ =x+ αy − y(x2 + y2) (27) This equation has the equilibrium (x, y) = (0, 0) with eigenvalues λ1,2 = α ± i. For bifurcations of equilibria, topological normal forms present general bifurcation diagrams. Together with system (27), let us introduce a new variable z = x + iy, z̄ = x − iy, |z|2 = zz̄ = x2 + y2 and the differential equation as ż = ẋ+ iẏ = α(x+ iy) + i(x+ iy)− (x+ iy)(x2 + y2) (28) Now we can rewrite system (27) in the following complex form: ż = (α+ i)z − z|z|2. (29) I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3721 This is the normal form of (27). To present the polar form , let z = reiθ(t) ż = ṙeiθ(t) + riθ̇eiθ(t) (30) Substituting into (6.12), we obtain the polar form: ṙeiθ(t) + riθ̇eiθ(t) = reiθ(α+ i− r2) (31) Now we can rewrite (6.14) as: ṙ =r(α− r2) θ̇ =1 (32) where r and θ are independent in uncouple system (32). The first equation has the equi- librium of (32) at r = 0 and α = 0. This equilibrium is a stable spiral for α < 0 and an unstable spiral for α > 0. It is possible to notice from (32) that the system has a unique and stable limit cycle for any α > 0 of radius r = √ α. More generally, when α passes through zero, under certain non-degeneracy conditions, we have a bifurcation in system (27) called the Andronov-Hopf bifurcation [14]. A cubic coefficient (called the first Lyapunov coefficient) in the complex form (29), determines the behaviour of (27) in the neighbourhood of a Hopf bifurcation point, if it is non-zero. We recall this in the following theorem: Theorem 1. [4, Topological normal form for the Hopf bifurcation, p100] Suppose two- dimensional dynamical system, ẋ = f(x, α), x ∈ R2, α ∈ R (33) with smooth f , has for all sufficiently small |α| the equilibrium x = 0 with eigenvalues λ1,2 = µ(α)± iω(α) where µ(0) = 0, ω(0) = ω0 > 0. Let the following conditions be satisfied: (H.1) l1(0) ̸= 0 , where l1 is the first Lyapunov coefficient (H.2) µ′(0) ̸= 0. Then there are invertible coordinate and parameter changes and a time reparametrization transforming (33) into the complex form ż = (l0 + i)z + l1z|z|2 +O(|z|4) (34) where l0 = 0 and l1 ̸= 0. I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3722 4.2. The normal form of the Bautin bifurcation The Bautin bifurcation (or generalised Hopf bifurcation) is codimension-2 of an equilib- rium. The Bautin point separates branches of subcritical and supercritical Andronov-Hopf bifurcations in the parameter plane. In particular, the equilibrium has purely imaginary eigenvalues and the first Lyapunov coefficient for Hopf bifurcation vanishes as Bautin bi- furcation conditions. Note when two limit cycles which collide and disappear via a saddle- node bifurcation of limit cycles. We can review the analysis of the Bautin bifurcation by stating the following theorem. Theorem 2. [4, Topological normal form for the Bautin bifurcation, p308] Suppose two- dimensional dynamical system, ẋ = f(x, α), x ∈ R2, α ∈ R2 (35) with smooth f , has for all sufficiently small |α| the equilibrium x = 0 with eigenvalues λ1,2 = µ(α)± iω(α) where µ(0) = 0, ω(0) = ω0 > 0 for α = 0, let Bautin bifurcation conditions hold: µ(0) = 0, l1(0) = 0 where l1(α)is the first Lyapunov coefficient. Let the following conditions be satisfied: (B.1) l2(0) ̸= 0 , where l2 is the second Lyapunov coefficient. (B.2) The map → (µ(α), l1(α)) T is regular at α = 0. Then there are invertible coordinate and parameter changes and a time reparametrization transforming (35) into the complex form ż = (l0 + i)z + l1z|z|2 + l2z|z|4 +O(|z|6) (36) where l0 = l1 = 0, l2 ̸= 0. Computation of Lyapunov numbers When a Hopf bifurcation point is detected as l0 = 0, the type of Hopf bifurcation of (33) corresponds to the sign of the first Lyapunov coefficient l1. We review the technique for computing the first and second Lyapunov coefficients at the bifurcation point. Also, the Taylor coefficients formula gkl of the quadratic terms in g(z, z̄, 0) can be expressed by the following [4, p311]: l1(α) = 1 2ω2 0 Re(ig20g11 + ω0g21), g20 = ⟨p,B(q, q)⟩, g11 = ⟨p,B(q, q̄)⟩, g02 = ⟨p,B(q̄, q̄)⟩ (37) I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3723 and g21 = ⟨p, C(q, q, q̄)⟩ − 2⟨p,B(q, A−1B(q, q̄)))⟩+ ⟨p,B(q̄, (2iω0E −A)−1B(q, q))⟩ + 1 iω0 ⟨p,B(q, q)⟩⟨p,B(q, q̄)⟩ − 2 iω0 |⟨p,B(q, q̄)⟩|2 − 1 3iω0 |⟨p,B(q̄, q̄)⟩|2 At the Bautin bifurcation point l1(0, 0) = 0, the formula gives l2(0, 0): 12l2(0, 0) = 1 ω0 Reg32 + 1 ω2 0 Im [ g20ḡ31 − g11(4g31 + 3ḡ22)− 1 3 g02(g40 + ḡ13)− g30g12 ] + 1 ω3 0 { Re [ g20 ( ḡ11(3g12 − ḡ30) + g02 ( ḡ12 − 1 3 g30 ) + 1 3 ḡ02g03 ) +g11 ( ḡ02 ( 3 5 ḡ30 + 3g12 ) + 1 3 g02ḡ03 − 4g11g30 )] + 3Im(g20g11)Img21 } + 1 ω4 0 { Im [ g11ḡ02(ḡ 2 20 − 3ḡ20g11 − 4g211) ] + Im(g20g11) [ 3Re(g20g11)− 2|g02|2 ]} (38) where g12 = ⟨p, C(q, q, q̄)⟩, g30 = ⟨p, C(q, q, q)⟩, g03 = ⟨p, C(q̄, q̄, q̄)⟩ This explains the general method for determining the critical normal form coefficients, which include third-order coefficients for the Hopf bifurcation and fifth-order coefficients for the Bautin bifurcation. The other formula for computation of the coefficients of the normal form are derived in [3]. 5. The normal form of bifurcations for Rothman’s carbon cycle model In this section, we present the normal form of Hopf and Bautin bifurcations for Roth- man’s carbon cycle model (7). Our analysis is based on the normal form theory of Hopf and Bautin bifurcations (see Subsection 5.1 and section 4.2). We use the conditions and highlight the Lyapunov coefficients that appear in this normal form. We use theorems 1 and 2 to identify the topological normal form for Hopf and Bautin bifurcations, and we aim to identify the Bautin bifurcation and classify its criticality by finding the appropriate coefficients in the normal form. We rewrite (7) here in the following complex form (36): ż = (l0 ± iw)z + l1z|z|2 + l2z|z|4 +O(|z|6) (39) where z ∈ C is a complex variable, and l0, l1,and l2 are real parameters. We now use the method discussed in sections 4.1-4.2 to determine the Bautin bifurcation in Rothman’s model (7). We implement this formula by using MAPLE to compute the normal form coefficients of Hopf and Bautin bifurcations for (7) [6]. I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3724 Figure 9: This shows the first l1 and second l2 Lyapunov coefficients along the two Hopf branches of (7) for varying cx ((a)-(b)) and µ ((c)-(d)). We use MAPLE to compute the l1 and l2 in the normal form at Hopf bifurcation point where l0 = 0. The type of Hopf bifurcation corresponds to the sign of the l1. The black stars indicate Bautin bifurcation points. We follow the technique in the previous section 4.1 to compute numerical values of of l1 and l2. We write the right-hand side of (7) as functions F1,2. We fix all parameters except µ or cx as bifurcation parameters, and F1,2 is then solved in order to get an equilibrium point (c, w). Then we compute the Jacobian matrix J and its transpose K at this point to get a value µ when the equilibrium point (c, w) has purely imaginary eigenvalues ±iω. We find the critical eigenvectors p of J and the critical eigenvectors q of K such that the normalisation ⟨p, q⟩ = 1. We introduce new complex variables X1,2 to evaluate the complex function H(z, z̄) as presented in (39). We need to expand H(z, z̄) to the fifth- order at (z, z̄) = (0, 0). Using (37) and (38) we are able to calculate numerical values of l1 and l2 are evaluated at the Hopf bifurcation point. Figure 9 shows the first l1 and second l2 Lyapunov coefficients along the two Hopf branches of (7) for varying cx (a)-(b)) and µ ((c)-(d)). Figure 10 shows the first l1 against second l2 Lyapunov coefficients along the two Hopf branches of (7) for varying cx and µ. In Figure 11, the first l1 and second l2 Lyapunov coefficients of (7) at a Hopf bifurcation point when cx and µ as the bifurcation parameters. we highlight the red curve (l1 < 0) and the blue curve (l1 > 0) correspond to the supercritical Hopf bifurcation and the subcritical Hopf bifurcation respectively. The two Bautin bifurcation points are represented by black I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3725 Figure 10: This shows the first l1 against second l2 Lyapunov coefficients along the two Hopf branches of (7) for varying cx and µ. We use MAPLE to compute the l1 and l2 in the normal form at Hopf bifurcation point where l0 = 0. The type of Hopf bifurcation corresponds to the sign of the l1. The two black stars indicate Bautin bifurcation points. points at l1 = 0, and we display these by leading at the sign of l2. 5.1. Saddle-node bifurcation of limit cycles This complex form (39) undergoes saddle-node bifurcation of limit cycles if l1l2 < 0 and l21 − 4l0l2 = 0 (40) Then, we attempt to approximate the location of the saddle-node limit cycles. In Figure 12, we used a numerical continuation software XPPAUT [2] to plot the bifurcation diagram of Rothman’s carbon cycle model (7) in the parameter plane (cx, µ), and we illustrated the fold limit cycle and the Hopf bifurcation of (7). We illustrate the approximation using the two methods in Figure 13. Figure 13 shows in the parameter plane (cx, µ) the curves of bifurcation of periodic orbits near the Bautin bifurcation (indicated by the red star). The saddle-node bifurcation of limit cycles are both of the Hopf curve. Figures 13(a) and 13(b) consider the left and right side of the Bautin bifurcation from Figure 12. These are computed in two ways: for the Hopf bifurcation curve, the XPPAUT calculations and the MAPLE calculations agree very well. On the other hand, the curve of the saddle-node bifurcation of limit cycles is calculated by two methods. For the Hopf curve, two method calculations have the same criticality and are lying on the same side of the Hopf curve in Figures 13(a) and 13(b), showing that both calculations indicate this is the subcritical Hopf bifurcation curve which lies to the left and to the right of the Bautin bifurcation from Figure 12. I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3726 60 70 80 90 100 110 cx 50 100 150 200 250 300 350 400 450 500 550 600 supercritical hopf bifurcation l 1 <0 subcritical hopf bifurcation l 1 >0 l 2 >0 bautin bifurcation l 1 =0 Figure 11: This illustrates the bifurcation diagram of (7) for the parameter plane (cx, µ). Two types of Hopf bifurcation are presented by blue (subcritical) and red (supercritical) curves. Bautin bifurcation points are indicated by black stars. l2 is negative along red parts of curve and l2 is positive along yellow parts of curve. All calculations were computed with MAPLE then MATLAB was used to plot the bifurcation diagram in the parameter space (cx, µ). In Figure 13(a), there is reasonable agreement between the XPPAUT and MAPLE calculations. For the Hopf bifurcation curve, the MAPLE calculations are closer to the Hopf line than the XPPAUT calculations which are presumably more accurate. The curve of the saddle-node bifurcation of limit cycles is plausible by using MAPLE calculations. This curve ends up being tangent to the XPPAUT curve that is presumably the correct curve, and then they meet at the Bautin bifurcation point. In Figure 13(b), the Maple calculations appear to be reasonable for the Hopf curve that agrees very well with XPPAUT. However, the MAPLE calculations of saddle-node bifurcation of limit cycles are less capable of resolving their behaviour near the second Bautin bifurcation. This could possibly be due to the normal form coefficients being much smaller in this case. 6. Conclusion These nonlinear oscillator models is the long timescale carbon cycle model of Rothman for long term oscillations in the global marine carbon cycle. We have confirmed the bifurcation behaviour found by Rothman [12] using analytical approximation. In Section 4, we derive the normal form of Hopf and Bautin bifurcations for Rothman’s carbon cycle model (7). Our analysis is based on the normal form theory of Hopf and Bautin bifurcations (see Subsection 5.1 and 4.2). We use the conditions and highlight the Lyapunov coefficients that appear in this normal form. We use Theorems 1 and 2 to I. Alraddadi / Eur. J. Pure Appl. Math, 17 (4) (2024), 3708-3729 3727 Figure 12: Bifurcation diagram for (7) computed with XPPAUT in the parameter plane (cx, µ). Hopf bifurcation for (7) is represented by blue curve. Two branches of saddle-node bifurcation of limit cycles for (7) are represented by black lines. identify the topological normal form for Hopf and Bautin bifurcations. From our results in this subsection, we are able to identify the Bautin bifurcation and classify its criticality by finding the appropriate coefficients in the normal form. We review the technique for computing Lyapunov coefficients. We use MAPLE to compute the first and second Lyapunov coefficients l1 and l2 in the normal form at a bifurcation point where l0 = 0. The type of Hopf bifurcation corresponds to the sign of the first Lyapunov coefficient l1. In Figure 11, we highlight the red curve (l1 < 0) and the blue curve (l1 > 0) correspond to the supercritical Hopf bifurcation and the subcritical Hopf bifurcation respectively. The Bautin bifurcation points are represented by black points at l1 = 0, and calculate the sign of l2 at these points to determine the Bautin bifurcation scenario. The main results are compared in Figure 13, showing the two-parameter bifurcation diagrams of Rothman’s carbon cycle model (7), with one of the diagrams computed with Maple and the other with XPPAUT. We use MAPLE to compute l2 at both points and find its negative for the parameters used, which confirms that both of the Bautin bifurca- tions are of supercritical type, then we confirm this using XPPAUT. Also, we attempt to approximate the location of the saddle-node bifurcation of limit cycles. We illustrate the approximation using the two methods in Figure 13. We compute the Hopf bifurcation and saddle-node bifurcation of limit cycles using XPPAUT and also using the value of l2 com- puted using MAPLE, and then we find that for one of these, the saddle-node bifurcation of limit cycles was approximated reasonably well in MAPLE using the fifth order normal form for one of them (the saddle-node bifurcation of limit cycles) but not the other. REFERENCES 3728 Figure 13: This shows (cx, µ) two-parameter bifurcation diagrams of (7). Left (a) and right (b) show curves of bifurcations of periodic orbits near the Bautin bifurcation (indicated by the red star). These bifurcation curves are computed using two methods: MAPLE and XPPAUT. Hopf bifurcation and saddle-node bifurcation of limit cycles are computed using MAPLE and represented by dashed (blue) and red lines. From the XPPAUT calculations, the Hopf bifurcation and saddle-node bifurcation of limit cycles are respectively shown using blue and red circles. Conflict of interest: The authors of the paper state that there is no conflict of interest. Data availability statement: All data generated or analyzed during the research investigation have been included in the paper. References [1] Michel Crucifix. Oscillators and relaxation phenomena in Pleistocene climate theory. Trans. R. Soc. A, 370:1140–1165, 2012. [2] Bard Ermentrout. Simulating, analyzing, and animating dynamical systems: a guide to XPPAUT for researchers and students. SIAM, 2002. [3] Yu. A. Kuznetsov. Numerical Normalization Techniques for All Codim 2 Bifurcations of Equilibria in ODE’s. SIAM Journal on Numerical Analysis, 36(4):1104–1124, 1999. REFERENCES 3729 [4] Yuri A. Kuznetsov. Elements of Applied Bifurcation Theory, Second Edition. Library, page 591, 1998. [5] Timothy M. Lenton, Hermann Held, Elmar Kriegler, Jim W. Hall, Wolfgang Lucht, Stefan Rahmstorf, and Hans Joachim Schellnhuber. Tipping elements in the earth’s climate system. Proceedings of the National Academy of Sciences, 105(6):1786–1793, 2008. [6] Stephen Lynch. Dynamical Systems with Applications using MATLAB®. Springer International Publishing, 2014. [7] Kirk A. Maasch and Barry Saltzman. A low-order dynamical model of global climatic variability over the full Pleistocene. Journal of Geophysical Research: Atmospheres, 95(D2):1955–1963, 1990. [8] J.Murray Mitchell. An overview of climatic variability and its causal mechanisms. Quaternary Research, 6(4):481–493, 1976. [9] Didier Paillard and Frédéric Parrenin. The Antarctic ice sheet and the triggering of deglaciations. Earth and Planetary Science Letters, 227(3):263–271, 2004. [10] Daniel Rothman. Earth’s carbon cycle: A mathematical perspective. Bulletin of the American Mathematical Society, 52(1):47–64, 2015. [11] Daniel H. Rothman. Thresholds of catastrophe in the earth system. Science Advances, 3(9):e1700906, 2017. [12] Daniel H. Rothman. Characteristic disruptions of an excitable carbon cycle. Proceed- ings of the National Academy of Sciences, 116(30):14813–14822, 2019. [13] Daniel H. Rothman, John M. Hayes, and Roger E. Summons. Dynamics of the neopro- terozoic carbon cycle. Proceedings of the National Academy of Sciences, 100(14):8124– 8129, 2003. [14] H Strogatz Steven. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering, 1994. [15] Anna S von der Heydt, Peter Ashwin, Charles D Camp, Michel Crucifix, Henk A Dijkstra, Peter Ditlevsen, and Timothy M Lenton. Quantification and interpretation of the climate variability record. Global and Planetary Change, 197:103399, 2021.