Leadership Parameters and Bifurcation of Political Unrest: a Mathematical Formalism with Cases Study

Leadership Parameters and Bifurcation of Political Unrest: A Mathematical Formalism with Case Studies of Soeharto, SBY, Jokowi, and Prabowo

Abstract

This paper develops a mathematical framework to analyze political unrest through the lens of bifurcation theory, integrating leadership parameters with nonlinear dynamic modeling. We propose seven leadership parameters---consensus, legitimacy, crisis management, narrative control, economic stability, elite management, and repression versus consensus---as structural determinants that modulate the coefficients of a coupled differential equation system describing trust (T), economic stress (E), protest intensity (P), and the emergence of alternative leaders (H). By embedding leadership parameters into the formalism, we derive analytical conditions for bifurcation, identify critical thresholds of external economic shocks (_E), and demonstrate how leadership quality shifts these thresholds. Numerical simulations are conducted for four Indonesian leaders---Soeharto, Susilo Bambang Yudhoyono (SBY), Joko Widodo (Jokowi), and Prabowo Subianto---illustrating distinct bifurcation patterns and the potential rise of "black horse" contenders under conditions of low legitimacy and high unrest. The results show that leadership qualities are not only qualitative traits but also mathematical modulators of systemic stability, providing new insights for the study of political dynamics in complex societies.

Main Motivation

Existing mathematical models of unrest focus primarily on external variables (economic shocks, trust erosion, protest contagion).

Leadership quality is usually treated qualitatively, not as a formal parameter in bifurcation analysis.

Our contribution is to synthesize leadership parameters with bifurcation theory, showing mathematically how governance style shifts critical thresholds for unrest and alters the likelihood of black horse emergence.

Indonesia provides an ideal case: four leaders with contrasting leadership profiles (Soeharto, SBY, Jokowi, Prabowo), spanning authoritarian, transitional democratic, and contemporary populist regimes.

Outline 

1. Introduction

Background on mathematical modeling of political unrest.

Gap: absence of leadership parameters in bifurcation models.

Contribution: formalism linking seven leadership parameters to systemic stability.

2. Mathematical Formalism

Definition of state variables: Trust (T), Economic Stress (E), Protest Intensity (P), Black Horse Potential (H).

Definition of seven leadership parameters and normalization.

Formal integration of leadership parameters into ODE system.

Derivation of normal form for bifurcation (saddle-node, Hopf).

3. Analytical Derivation

Fixed point analysis.

Jacobian derivation and eigenvalue conditions.

Critical thresholds (_E, _crit) as function of leadership parameters.

4. Numerical Simulations

Parameter calibration for Soeharto, SBY, Jokowi, Prabowo.

Bifurcation diagrams (P* vs _E).

Sensitivity analysis: how legitimacy, narrative, and repression shift _crit.

Simulation of black horse emergence (H).

5. Case Study: Indonesia

Soeharto: high consolidation, low legitimacy abrupt tipping.

SBY: high legitimacy, moderate coalition resilient.

Jokowi: strong narrative, elite cooptation stable until sudden shift.

Prabowo: low legitimacy, high repression low _crit, unstable under shocks.

6. Discussion

Theoretical implications: leadership as mathematical modulator.

Policy implications: stability trade-offs between repression and legitimacy.

Limitations and directions for future work.

7. Conclusion

Leadership quality is quantifiable as parameters in nonlinear dynamics.

Integration of leadership parameters with bifurcation theory provides novel insights into political unrest and regime transitions.

I. Introduction

A. Background on Mathematical Modeling of Political Unrest

The study of political unrest has long attracted scholars from political science, sociology, and economics. Yet in recent decades, the analytical tools of mathematics---particularly nonlinear dynamics and bifurcation theory---have offered a fresh perspective for understanding how societies transition from stability to instability. Unlike purely descriptive approaches, mathematical modeling seeks to formalize the underlying mechanisms that drive collective action, state fragility, and the sudden emergence of mass protests.

Traditional models of unrest often treat social dynamics as systems sensitive to external shocks, such as abrupt economic crises, food price inflation, or policy missteps that trigger widespread grievances. In such frameworks, unrest is viewed as the amplification of discontent through social contagion, where individual grievances accumulate until a threshold is exceeded. This threshold-based perspective has been mathematically expressed through differential equations that capture feedback loops between trust in government, economic stress, and protest mobilization. Once a tipping point is reached, the system can undergo a bifurcation---a qualitative change in its equilibrium behavior---resulting in the rapid escalation of demonstrations or even regime change.

Pioneering contributions in this domain have drawn upon complex systems theory, treating political unrest as analogous to phase transitions in physics. In these models, small perturbations can have negligible effects in stable regimes but can produce disproportionate consequences when the system is near criticality. This insight explains why protests sometimes dissipate harmlessly, while at other times they ignite waves of political upheaval. Furthermore, the inclusion of stochastic terms in these models highlights the role of noise-induced transitions, where seemingly minor incidents (e.g., the death of a protester, a controversial policy announcement) can catalyze broader instability if the system is already fragile.

The mathematical formalism of unrest has therefore evolved to emphasize two intertwined dimensions:

1. Structural variables -- such as economic conditions, institutional resilience, and trust in governance.

2. Dynamic processes -- including protest contagion, state response (repression or accommodation), and elite defection.

By combining these dimensions, bifurcation theory provides a rigorous lens through which political scientists and applied mathematicians can jointly analyze regime stability. The central idea is that societal stability is not linear; small incremental changes in underlying variables can push the system toward critical thresholds where qualitative transformations occur.

In recent years, this approach has been further refined by introducing concepts such as the "black horse" phenomenon, where emergent political actors gain sudden traction amidst systemic instability. These models capture how the interplay of declining legitimacy, rising grievances, and weakened elite cohesion creates opportunities for outsider figures to rapidly rise to prominence.

Taken together, the mathematical modeling of political unrest moves beyond narrative explanation and provides formal predictive tools. It allows not only for the identification of tipping points but also for the simulation of policy interventions that may shift critical thresholds and alter the trajectory of instability.

B. Gap: Absence of Leadership Parameters in Bifurcation Models

While bifurcation models of political unrest have provided powerful insights into how societies transition from order to instability, most of these frameworks share a common limitation: they treat leadership as an exogenous or background factor, rather than as a measurable and integral component of the system. In classical formulations, regime stability is modeled almost exclusively in terms of economic stress, institutional resilience, and social contagion, while the qualities of political leaders---such as their ability to consolidate coalitions, communicate narratives, or manage crises---remain outside the formal mathematical structure.

This omission creates a conceptual gap. In practice, leadership is often the decisive variable that shapes whether external shocks escalate into systemic breakdown or are absorbed without destabilization. For instance, two societies facing comparable economic stress may diverge dramatically in outcomes depending on the legitimacy of leadership, the degree of repression versus consensus, or the effectiveness of narrative control. Yet, in most bifurcation models, these differences are implicitly assumed rather than explicitly parameterized.

The absence of leadership parameters also constrains the explanatory power of existing models. They can describe when a tipping point might occur under certain economic or institutional conditions, but they cannot fully explain why similar shocks produce unrest in one regime and not in another. Nor can they simulate the role of leadership style in shifting critical thresholds. This limitation becomes evident when comparing political leaders across different eras or contexts: authoritarian rulers, democratic consensus-builders, and populist strategists all face economic shocks, but their governance qualities profoundly alter the dynamics of unrest.

In addition, ignoring leadership overlooks the feedback loop between perception and stability. Leadership is not simply a passive background condition; it actively shapes trust (T), modifies elite cohesion (K), influences protest intensity (P) through narrative framing, and even conditions the space for emergent "black horse" actors. Without incorporating these dimensions, mathematical models risk reducing complex socio-political processes to overly mechanical systems, detached from the lived reality of governance.

Therefore, there is a pressing need to extend bifurcation models by embedding leadership qualities as formal parameters. By doing so, we can move beyond models that focus only on structural stress and develop a more holistic framework in which both external shocks and internal leadership dynamics jointly determine systemic stability. This integration is particularly vital for understanding contemporary politics, where leadership style can mean the difference between resilience and collapse under similar external conditions.

C. Contribution: Formalism Linking Seven Leadership Parameters to Systemic Stability

To address the gap outlined above, this study introduces a formal mathematical synthesis that integrates leadership qualities directly into bifurcation models of political unrest. We operationalize leadership not as a vague qualitative trait, but as a structured set of seven parameters that can be normalized, weighted, and systematically embedded within a nonlinear dynamical system. These parameters are:

1. Consensus-building capacity,

2. Legitimacy and public trust,

3. Crisis management effectiveness,

4. Narrative and communication control,

5. Economic governance and stability,

6. Elite management and coalition maintenance, and

7. Balance of repression versus consensus.

Each of these parameters is mathematically mapped onto the coefficients of coupled differential equations that describe the evolution of trust (T), economic stress (E), protest intensity (P), and the potential rise of black horse leaders (H). By doing so, leadership no longer remains exogenous but becomes a quantifiable driver of systemic dynamics.

This integration yields several contributions. First, it enables the derivation of analytical conditions for bifurcation in which leadership parameters shift the critical threshold (_crit) of external shocks required to trigger unrest. For example, higher legitimacy and effective narrative control increase the system's resilience by raising _crit, while overreliance on repression may lower it, making unrest more likely even under moderate shocks.

Second, the formalism allows for comparative simulations across leadership profiles. By calibrating the seven parameters to empirical cases, such as Soeharto, Susilo Bambang Yudhoyono (SBY), Joko Widodo (Jokowi), and Prabowo Subianto, we can illustrate how different governance styles produce distinct bifurcation patterns. This demonstrates not only the explanatory power of the model but also its empirical relevance to real-world political trajectories.

Third, the inclusion of leadership parameters provides a bridge between political science and applied mathematics. It offers a framework in which qualitative assessments of leadership can be systematically translated into quantitative models, facilitating interdisciplinary dialogue and predictive capacity.

Finally, this contribution advances the study of black horse emergence. By explicitly linking leadership failures---such as weak legitimacy or poor elite management---to the growth potential of outsider figures, the model captures the nonlinear interaction between internal governance dynamics and systemic fragility.

In summary, the formalism developed here reframes leadership as a mathematical modulator of political stability. It demonstrates that the fate of regimes under economic or social stress cannot be understood without accounting for how leaders manage legitimacy, elites, narratives, and repression. This synthesis moves the field toward a more comprehensive and predictive science of political unrest.

II. Mathematical Formalism

A. Definition of State Variables: Trust (T), Economic Stress (E), Protest Intensity (P), and Black Horse Potential (H)

To construct a rigorous mathematical model of political unrest, we begin by formalizing the state variables that capture the dynamic evolution of societal stability. These variables are defined on normalized scales, bounded between 0 and 1, to ensure comparability across contexts and facilitate numerical simulation.

1. Trust (T)
Trust represents the degree of public confidence in governing institutions and political leadership. A value of T = 1 indicates complete trust and widespread legitimacy, while T = 0 reflects total erosion of confidence, often preceding systemic collapse. Trust functions as a stabilizing variable, exerting downward pressure on protest intensity and buffering the effects of economic stress. It evolves endogenously through feedback loops with both leadership parameters (e.g., legitimacy, narrative control) and external shocks.
2. Economic Stress (E)
Economic stress captures the level of socioeconomic pressure faced by society, such as inflation, unemployment, or income inequality. Unlike trust, which reflects perception and legitimacy, economic stress reflects material conditions. E = 0 denotes a state of economic ease or prosperity, while E = 1 corresponds to severe economic crisis. This variable is also directly influenced by exogenous shocks, represented by a control parameter (_E), making it a natural driver of bifurcation.
3. Protest Intensity (P)
Protest intensity measures the extent and persistence of collective action against the regime. It is conceptualized as a dynamic population-level response, encompassing demonstrations, strikes, or other forms of political mobilization. P = 0 signifies total quiescence, while P = 1 represents maximal unrest approaching revolutionary levels. Protest intensity grows with economic stress and declining trust, but is counteracted by leadership mechanisms such as consensus-building, repression, and elite cooptation.
4. Black Horse Potential (H)
Black horse potential refers to the latent probability of emergent outsider leaders or alternative elites gaining traction during instability. Unlike the other variables, H is not directly observable at all times but manifests strongly during moments of systemic weakness. H = 0 means no viable outsider figures exist or gain traction, while H = 1 corresponds to a scenario in which alternative leaders become dominant challengers to the regime. The growth of H is stimulated by high protest intensity and low trust, while being suppressed when leadership successfully maintains legitimacy and elite cohesion.
Together, these four variables form the core dynamic system. They are interconnected through nonlinear feedbacks: economic stress undermines trust, low trust amplifies protest intensity, rising protests open opportunities for black horse actors, and the success or failure of these processes depends on the embedded leadership parameters. By treating these variables as dynamical states, we lay the foundation for deriving the formal equations governing systemic stability.

B. Definition of Seven Leadership Parameters and Normalization

To integrate leadership quality into a mathematical model of political unrest, we define seven core parameters. Each parameter represents a distinct dimension of governance and is normalized to a continuous scale between 0 (absence of the quality) and 1 (maximum strength). This normalization facilitates both cross-leader comparison and mathematical embedding within the dynamical system.

1. Consensus-Building Capacity (C)
Definition: The ability of leaders to mediate conflicts, forge compromises, and include diverse societal actors in policymaking.

Normalization: C = 0 represents total exclusionary politics, while C = 1 reflects highly inclusive and participatory governance.

2. Legitimacy and Public Trust (L)
Definition: The degree of acceptance of leadership by the governed, grounded in perceived fairness, transparency, and historical credibility.

Normalization: L = 0 implies no legitimacy (pervasive distrust), while L = 1 implies universal trust and acceptance.

3. Crisis Management Effectiveness (M)
Definition: The competence of leaders in addressing sudden shocks such as natural disasters, pandemics, or political scandals.

Normalization: M = 0 denotes ineffectiveness and paralysis, while M = 1 denotes rapid, coordinated, and effective responses.

4. Narrative and Communication Control (N)
Definition: The capacity to frame political discourse, shape public perception, and dominate narratives through media and rhetoric.

Normalization: N = 0 indicates fragmented communication and loss of narrative control, while N = 1 indicates strong, coherent, and persuasive messaging.

5. Economic Governance and Stability (G)
Definition: The leader's ability to maintain macroeconomic stability, ensure growth, and mitigate inequalities that drive discontent.

Normalization: G = 0 signals economic mismanagement, while G = 1 reflects robust and equitable economic stewardship.

6. Elite Management and Coalition Maintenance (E_c)
Definition: The ability to manage political elites, prevent defection, and maintain coalition discipline in support of the regime.

Normalization: E_c = 0 indicates fractured elite networks and frequent defections, while E_c = 1 indicates cohesive and loyal elite coalitions.

7. Repression versus Consensus Balance (R)
Definition: The strategic mix between coercive repression and consensual inclusion in handling dissent.

Normalization: R = 0 represents reliance on pure repression (high risk of backfire), while R = 1 represents a balanced and calibrated approach where coercion is minimal and consensus prevails.

Normalization Rationale

All seven parameters are bounded between 0 and 1 for comparability.

Calibration can be based on qualitative coding (expert judgment, historical records) or quantitative proxies (indices of governance, corruption perception, freedom scores).

These parameters are not independent; they interact nonlinearly within the model, shaping coefficients that determine how trust (T), economic stress (E), protest intensity (P), and black horse potential (H) evolve.

By defining leadership in this normalized, multidimensional way, the model transforms leadership qualities from descriptive attributes into mathematical levers. This allows us to formally analyze how leadership alters bifurcation thresholds and the stability of political systems under external shocks.

C. Formal integration of the seven leadership parameters into an ODE system

In this subsection we present a concrete, transparent, and analytically tractable way to embed the seven leadership parameters defined in Section 2.B into a coupled system of ordinary differential equations (ODEs) that governs the time evolution of the state variables (trust), (economic stress), (protest intensity), and (black-horse potential). The presentation below is written to satisfy the precision required for international peer review: every variable and coefficient is defined, assumptions are stated, and the mapping from leadership parameters to model coefficients is explicit.

1. Notation and preliminaries

State variables (all normalized to the unit interval):

T(t)[0,1]T(t)\in[0,1]T(t)[0,1] --- public trust / legitimacy,
E(t)[0,1]E(t)\in[0,1]E(t)[0,1] --- economic stress (0 = none, 1 = extreme stress),
P(t)[0,1]P(t)\in[0,1]P(t)[0,1] --- protest intensity (0 = quiescence, 1 = maximal unrest),
H(t)[0,1]H(t)\in[0,1]H(t)[0,1] --- black-horse emergence potential.
Leadership parameters (normalized to [0,1][0,1][0,1]; Section 2.B):

CCC --- consensus-building capacity,
LLL --- legitimacy (public trust baseline),
MMM --- crisis management effectiveness,
NNN --- narrative/communication control,
GGG --- economic governance effectiveness,
EcE_cEc --- elite management / coalition maintenance,
RRR --- repression vs. consensus balance (higher = more consensus, lower = more repression).
External control parameter:

ER0\mu_E\in\mathbb{R}_{\ge0}ER0 --- exogenous economic shock magnitude (driver of EEE); treated as a bifurcation parameter in bifurcation diagrams.

We introduce short-hand mappings from leadership parameters to model coefficients. The mapping is chosen to be monotonic and dimensionless; it is sufficiently general to be adapted or re-calibrated for empirical work.

2. Model structure --- core assumptions

1. Trust growth/decay. Trust increases via economic recovery and successful narrative/capacity to manage crises, and decreases with protest intensity and perceived repression/backfire. Leadership parameters L,N,M,GL,N,M,GL,N,M,G act to bolster trust; repression RRR modulates backfire effects.
2. Economic stress dynamics. Economic stress is driven by exogenous shocks E\mu_EE and endogenous feedback from protests (disruption of economic activity) and poor crisis management. Good economic governance GGG and effective crisis management MMM speed damping of EEE.
3. Protest dynamics. Protest follows logistic-style growth (capacity-limited mobilization) and is positively driven by economic stress and lowered trust; it is suppressed by elite cooptation KKK and by legitimacy LLL and narrative NNN. Repression (low RRR) can suppress protests short-term but produces legitimacy loss and thus can increase protest indirectly (backfire).
4. Black-horse dynamics. HHH grows when protests are sustained and trust is low; elite fragmentation (low EcE_cEc) and poor consensus capacity CCC raise the growth potential for outsider leaders. Strong elite management EcE_cEc and high LLL suppress HHH.
5. Saturation and nonlinearity. Many processes saturate (limited population to mobilize, limited trust ceiling). Sigmoidal or saturating functions (e.g., logistic, Hill functions, tanh\tanhtanh) are used to capture these effects.
6. Optional stochastic forcing. To capture real-world idiosyncratic shocks and noise-induced tipping, we allow for additive white noise terms in each equation; the deterministic skeleton below is our primary object for bifurcation analysis.

3. Explicit ODE system (deterministic skeleton)

We propose the following system (all time derivatives are ):

\boxed{%

\begin{aligned}

\dot T &= \; \underbrace{\alpha_T \, S_E(E;\theta_E,\kappa_E)\,(L + \sigma_N N + \sigma_M M)\,(1-T)}_{\text{trust recovery via economy, narrative, crisis mgmt}} 

\;-\; \underbrace{\beta_T\, P}_{\text{trust erosion by protest}} 

\;-\; \underbrace{\gamma_T(1-R)\,P}_{\text{repression backfire term}} \\

\dot E &= \; \mu_E \;-\; \underbrace{\lambda_E(G,M)\,E}_{\text{damping by economic policy \& crisis mgmt}}

\;+\; \underbrace{\phi_P\, P}_{\text{protest worsens economy}} \\

\dot P &= \; \underbrace{\rho_P\,P\left(1-\dfrac{P}{P_{\max}}\right)}_{\text{intrinsic mobilization}} 

\;+\; \underbrace{\eta_E\, S_E(E;\theta_E,\kappa_E)}_{\text{economic grievance channel}} 

\;-\; \underbrace{\delta_T\, T}_{\text{trust dampening}} 

\;-\; \underbrace{\kappa_K\,K(C,E_c,R)\,P}_{\text{elite cooptation / suppression}} \\

\dot H &= \; \underbrace{\alpha_H\,P}_{\text{opportunity from protest}} \;+\; \underbrace{\beta_H\,E}_{\text{economic opening}} \;-\; \underbrace{\gamma_H\,T}_{\text{trust suppression}} \;-\; \underbrace{\lambda_H\,E_c}_{\text{elite containment}}.

\end{aligned}

}

Each term and coefficient is described in detail below.

4. Functional forms and coefficient mappings

(a) Economic trust / grievance coupling

We use a smooth saturating function for the economic channel:

S_E(E;\theta_E,\kappa_E) \;=\; \dfrac{E^{\kappa_E}}{E^{\kappa_E} + \theta_E^{\kappa_E}},

\qquad \theta_E\in(0,1),\ \kappa_E>0.

is close to 0 for small and approaches 1 for large ; it is conveniently differentiable for Jacobian calculation.

(b) Elite cooptation function .
Elite cooptation reduces protest pressure by absorbing potential elite challengers and tying resources to the incumbent. We propose a bounded saturating functional form:

K(C,E_c,R) \;=\; K_{\max}\,\tanh\big(\rho_K\,[w_C\,C + w_{E}\,E_c + w_R\,R]\big),

\qquad K_{\max}>0,\ \rho_K>0,

(c) Damping of economic stress .
Economic governance and crisis management accelerate recovery:

\lambda_E(G,M) \;=\; \lambda_{0} \;+\; \lambda_{G}\,G \;+\; \lambda_{M}\,M,

\qquad \lambda_{0}>0,\ \lambda_G,\lambda_M\ge0.

(d) Trust recovery amplitude and narrative / crisis weightings.
We set

\alpha_T>0,\qquad \text{and}\qquad \alpha_T(L,N,M) \equiv \alpha_T\, (L + \sigma_N N + \sigma_M M),

(e) Repression backfire coefficient.
The term encodes that the more repressive the regime (smaller ), the larger the trust erosion per unit protest (backfire). .

(f) Protest suppression by trust and cooptation.
sets how strongly trust reduces protest. sets the strength of suppression from elite cooptation ; if repression is used in combination with cooptation, K may be large short term but lead to long-term increases in grievance; this trade-off is captured when calibrating and the dependence of other coefficients on .

(g) Black-horse growth terms.
respectively measure sensitivity of to protest, economic stress, trust, and elite cohesion. For example, is large when elite management strongly prevents outsider leader formation.

5. Interpretation of coefficient dependence on leadership parameters

For transparency we give explicit example mappings (linear monotonic maps, to be adapted empirically):

\begin{aligned}

\delta_T &= \delta_0 + \delta_L\, L + \delta_N\, N,\\

\kappa_K  &= \kappa_0 + \kappa_{E}\, E_c + \kappa_C\, C - \kappa_R\,(1-R),\\

\lambda_E &= \lambda_0 + \lambda_G\, G + \lambda_M\, M,\\

\alpha_T &= \alpha_0\,(L + \sigma_N N + \sigma_M M),\\

\gamma_T &= \gamma_0\,(1-R),\\

\lambda_H &= \lambda_{H0}\, E_c.

\end{aligned}

These mappings ensure:

Higher generally increase stabilizing coefficients (trust recovery, economic damping, cooptation efficiency), thereby shifting the system toward resilience.

Lower (higher repression) increases short-term suppression but increases (backfire) which erodes faster given protests.

6. Boundedness and invariance

Because all state variables are normalized to [0,1][0,1][0,1], the model must ensure forward invariance of the unit cube. For this we choose coefficient values and saturating functions so that the vector field points inward on the boundaries T=0,1T=0,1T=0,1, E=0,1E=0,1E=0,1, P=0,1P=0,1P=0,1, H=0,1H=0,1H=0,1. In practice this is enforced by (i) logistic form for PPP with carrying capacity Pmax1P_{\max}\le1Pmax1, (ii) linear damping for EEE with sufficient E\lambda_EE, and (iii) multiplicative factors (1T)(1-T)(1T) in trust recovery.

7. Stochastic extension (optional for realism)

To model noise-induced transitions and to conduct ensemble simulations, we consider the stochastic system

dX=F(X;E,)dt+(X)dWt,X=(T,E,P,H),dX = F(X;\mu_E,\Theta)\,dt \;+\; \Sigma(X)\, dW_t, \qquad X=(T,E,P,H)^\top,dX=F(X;E,)dt+(X)dWt,X=(T,E,P,H),

where FFF is the deterministic vector field given above, \Theta denotes the vector of leadership-dependent parameters, WtW_tWt is a standard multi-dimensional Wiener process, and (X)\Sigma(X)(X) is a noise amplitude matrix (often taken diagonal with small entries). Stochastic simulations are used to study noise-induced tipping near critical parameter values Ecrit\mu_E \approx \mu_{\mathrm{crit}}Ecrit.

8. Fixed points and Jacobian (prelude to bifurcation analysis)

Fixed points X=(T,E,P,H)X^*=(T^*,E^*,P^*,H^*)X=(T,E,P,H) satisfy F(X;E,)=0F(X^*;\mu_E,\Theta)=0F(X;E,)=0. The Jacobian matrix J(X;E,)=F/XXJ(X^*;\mu_E,\Theta)=\partial F/\partial X\big|_{X^*}J(X;E,)=F/XX is a 444\times444 matrix whose spectrum determines linear stability. Explicit symbolic expressions for the Jacobian entries are straightforward to compute given the chosen functional forms (logistic, tanh\tanhtanh, Hill functions). For bifurcation analysis we examine how eigenvalues i(E,)\lambda_i(\mu_E,\Theta)i(E,) cross the imaginary axis as E\mu_EE (or other control parameters) vary:

Saddle-node (fold) bifurcation: a real eigenvalue crosses zero and a pair of fixed points collides/annihilates --- typically associated with abrupt tipping in PPP.
Hopf bifurcation: a complex conjugate pair crosses the imaginary axis small amplitude oscillations / protest waves may appear.
Using center manifold reduction and normal form calculations (standard textbook procedures), one can reduce the dynamics near a nonhyperbolic fixed point to canonical scalar forms such as x=r+x2\dot x = r + x^2x=r+x2 (fold) or z=(+i)z+z2z\dot z = (\alpha + i\omega) z + \beta |z|^2 zz=(+i)z+z2z (Hopf). In these reduced coordinates the bifurcation parameter rrr is an explicit function of E\mu_EE and the leadership-dependent coefficients \Theta; hence leadership parameters shift the location and nature of bifurcations.

9. Practical calibration and identifiability remarks

Calibration sources. Leadership parameter values C,L,M,N,G,Ec,RC,L,M,N,G,E_c,RC,L,M,N,G,Ec,R can be informed by expert coding, composite governance indices (e.g., World Bank governance indicators scaled), media sentiment measures, or historical qualitative coding. Coefficients such as G,E,0\lambda_G,\kappa_E,\delta_0G,E,0 are calibrated by fitting observed time series of protest incidence and economic indicators where available, using likelihood methods or Bayesian inference.
Identifiability. The mapping from leadership parameters to coefficients must avoid collinearity that compromises parameter identifiability; regularization and sensitivity analysis are recommended.
10. Summary of the mapping logic

Leadership parameters enter the ODE system by modifying coefficient functions that govern rates of trust recovery/deterioration, strength of cooptation, damping of economic stress, and suppression or facilitation of protest and black-horse emergence.
Analytically, this mapping makes the leadership vector =(C,L,M,N,G,Ec,R)\Theta=(C,L,M,N,G,E_c,R)=(C,L,M,N,G,Ec,R) a secondary control set: while E\mu_EE is the primary exogenous bifurcation parameter, \Theta determines crit()\mu_{\mathrm{crit}}(\Theta)crit(), the critical shock amplitude required to cause bifurcation.
The formalism therefore yields explicit, testable statements of the form:
 crit=(),\mu_{\mathrm{crit}} = \Phi(\Theta),crit=(),
 where \Phi is obtained by solving the fixed-point conditions and locating parameter values where the Jacobian acquires zero or purely imaginary eigenvalues.
This completes the formal presentation of how the seven leadership parameters are integrated into a mathematically precise ODE model for political unrest. The next logical steps (Section 3) are (i) fixed-point computation and Jacobian derivation in closed form for a reduced parametrization, (ii) normal-form reduction to obtain analytic expressions for in leading order, and (iii) numerical continuation and bifurcation diagrams to visualize how different leadership profiles (Soeharto, SBY, Jokowi, Prabowo) shift critical thresholds. If you wish, I will proceed immediately to (i) and (ii) with explicit symbolic derivations for a reduced three-variable core (e.g., ) to make the bifurcation analysis fully explicit.

Below we derive the normal forms for the two generic bifurcations we care about (saddle-node / fold and Hopf) for our unrest system. I keep the derivation general (valid for the 3--4D system introduced earlier) but give an explicit workflow and the exact algebraic expressions you must compute to obtain the reduced normal forms. That makes the derivation verifiable and ready to be implemented symbolically (e.g., in Mathematica/Maple/SymPy) or numerically.

Notation & setup

State vector XRnX\in\mathbb{R}^nXRn with n=3n=3n=3 or 444 (we typically use the reduced core X=(T,E,P)X=(T,E,P)^\topX=(T,E,P) or full X=(T,E,P,H)X=(T,E,P,H)^\topX=(T,E,P,H)).
Deterministic vector field X=F(X;,) \dot X = F(X;\mu,\Theta)X=F(X;,), where \mu is the primary scalar bifurcation parameter (here =E\mu=\mu_E=E) and \Theta denotes the leadership parameter vector.
Fixed point X(,)X^*(\mu,\Theta)X(,) satisfies F(X;,)=0F(X^*;\mu,\Theta)=0F(X;,)=0.
Jacobian J(X;,)=DXFXJ(X^*;\mu,\Theta) = D_XF|_{X^*}J(X;,)=DXFX. Eigenvalues i()\lambda_i(\mu)i() of JJJ determine linear stability.
We will derive the normal form in two cases:

a. Saddle-node (fold): a simple real eigenvalue passes through zero. The reduced normal form near the bifurcation is x=r+x2+...\dot x = r + \alpha x^2 + \ldotsx=r+x2+... (one-dimensional).
b. Hopf: a complex conjugate pair crosses the imaginary axis at i0\pm i\omega_0i0. The reduced normal form on the center manifold is z=(+i0)z+1zz2+...\dot z = (\eta + i\omega_0)z + \ell_1 z|z|^2 + \ldotsz=(+i0)z+1zz2+... where the real part of 1\ell_11 (the first Lyapunov coefficient) determines super/subcritical nature.
a. Saddle-node (fold) --- reduction & normal form

1. Existence / linear condition

Assume at =c\mu=\mu_c=c there is an equilibrium X=XcX^*=X_cX=Xc such that:

F(Xc;c)=0F(X_c;\mu_c)=0F(Xc;c)=0.
Jc:=J(Xc;c)J_c := J(X_c;\mu_c)Jc:=J(Xc;c) has a single simple zero eigenvalue, and all other eigenvalues have nonzero real parts. Let vvv be a right nullvector and www the corresponding left nullvector (row vector) normalized so that wv=1w^\top v = 1wv=1.
2. Nondegeneracy & transversality conditions

Two standard nondegeneracy conditions must hold for a classical saddle-node:

Nondegeneracy (quadratic nonlinearity):
a:=12wD2F(Xc)[v,v]0.a := \tfrac{1}{2}\, w^\top D^2F(X_c)[v,v] \;\neq\; 0.a:=21wD2F(Xc)[v,v]=0.
Here D2F(Xc)[v,v]D^2F(X_c)[v,v]D2F(Xc)[v,v] is the vector whose iii-th component is j,k2Fixjxk(Xc)vjvk\sum_{j,k}\frac{\partial^2 F_i}{\partial x_j\partial x_k}(X_c)\,v_j v_kj,kxjxk2Fi(Xc)vjvk.
Transversality (parameter dependence pushes eigenvalue through zero):
b:=wF(Xc;c)0.b := w^\top \partial_\mu F(X_c;\mu_c) \;\neq\; 0.b:=wF(Xc;c)=0.

If a0a\neq0a=0 and b0b\neq0b=0, the local dynamics on the one-dimensional center manifold reduce (after smooth coordinate change and parameter re-scaling) to:

y=b(c)+ay2+higher order terms.\dot y \;=\; b(\mu-\mu_c) + a\, y^2 + \text{higher order terms}.y=b(c)+ay2+higher order terms.

Usual normal form scaling yields the canonical fold:

x=r+x2,rb(c),xy.\dot x = r + x^2,\qquad r \propto b(\mu-\mu_c),\;x\propto y.x=r+x2,rb(c),xy.

Interpretation: the sign of aba bab tells whether equilibria are created/destroyed for \mu increasing past c\mu_cc. For our problem, computing aaa and bbb (explicitly in terms of partial derivatives) gives the analytic expression for crit()\mu_{\mathrm{crit}}(\Theta)crit() locally.

3. How to compute aaa and bbb concretely

Solve Jcv=0J_c v = 0Jcv=0 for right nullvector vvv.
Solve wJc=0w^\top J_c = 0wJc=0 for left nullvector www, normalize wv=1w^\top v = 1wv=1.
Compute the vector D2F(Xc)[v,v]D^2F(X_c)[v,v]D2F(Xc)[v,v] with components
(D2F(Xc)[v,v])i=j,k=1n2Fixjxk(Xc)vjvk.\big(D^2F(X_c)[v,v]\big)_i \;=\; \sum_{j,k=1}^n \frac{\partial^2 F_i}{\partial x_j\partial x_k}(X_c)\, v_j v_k.(D2F(Xc)[v,v])i=j,k=1nxjxk2Fi(Xc)vjvk.

Then a=12wD2F(Xc)[v,v]a = \tfrac{1}{2}\, w^\top D^2F(X_c)[v,v]a=21wD2F(Xc)[v,v].
Compute F\partial_\mu FF (partial derivative of vector field w.r.t. \mu); evaluate at (Xc,c)(X_c,\mu_c)(Xc,c) and form b=wF(Xc;c)b = w^\top \partial_\mu F(X_c;\mu_c)b=wF(Xc;c).
If aaa and bbb satisfy the inequalities above, you have the fold normal form.

Practical note: For our model the only explicit dependence on \mu is in E\dot EE via the additive E\mu_EE term, so F=eE\partial_\mu F = e_EF=eE unit vector in EEE-direction; thus b=wEb = w_Eb=wE (the EEE-component of the left nullvector), making the transversality test computationally straightforward.

b. Hopf bifurcation --- reduction & normal form

1. Linear condition

Assume at =h\mu=\mu_h=h there is equilibrium X=XhX^*=X_hX=Xh such that:

F(Xh;h)=0F(X_h;\mu_h)=0F(Xh;h)=0.
Jacobian JhJ_hJh has a simple pair of purely imaginary eigenvalues 1,2(h)=i0\lambda_{1,2}(\mu_h)=\pm i\omega_01,2(h)=i0 with 0>0\omega_0>00>0; all other eigenvalues have nonzero real parts. The pair is simple (algebraic multiplicity 1).
2. Transversality & nondegeneracy

Transversality (nonzero drift of eigenvalue): ddRe()=h0.\left. \dfrac{d}{d\mu} \mathrm{Re}\,\lambda(\mu)\right|_{\mu=\mu_h} \neq 0.ddRe()=h=0. This can be checked using the formula:
ddh=wF(Xh;h)wv,\left.\dfrac{d\lambda}{d\mu}\right|_{\mu_h} = \dfrac{w^\top \partial_\mu F(X_h;\mu_h)}{w^\top v},ddh=wvwF(Xh;h),
where vvv is the complex right eigenvector Jhv=i0vJ_h v = i\omega_0 vJhv=i0v and www is the corresponding left eigenvector normalized wv=1w^\top v = 1wv=1. The real part of this derivative must be nonzero.
First Lyapunov coefficient 1\ell_11 (nondegeneracy): compute 1\ell_11 at (Xh,h)(X_h,\mu_h)(Xh,h). If Re10\mathrm{Re}\,\ell_1 \neq 0Re1=0, the Hopf is generically nondegenerate. The sign of Re1\mathrm{Re}\,\ell_1Re1 determines supercritical (Re1<0\mathrm{Re}\,\ell_1<0Re1<0) or subcritical (Re1>0\mathrm{Re}\,\ell_1>0Re1>0) Hopf.
3. Formula for the first Lyapunov coefficient 1\ell_11

Follow the standard multilinear form notation. Let A=JhA = J_hA=Jh. Define bilinear and trilinear maps:

B(u,v)=D2F(Xh)[u,v]B(u,v) = D^2F(X_h)[u,v]B(u,v)=D2F(Xh)[u,v] vector in Rn\mathbb{R}^nRn.
C(u,v,w)=D3F(Xh)[u,v,w]C(u,v,w) = D^3F(X_h)[u,v,w]C(u,v,w)=D3F(Xh)[u,v,w].
Let vCnv\in\mathbb{C}^nvCn be the right eigenvector for eigenvalue i0i\omega_0i0, and wCnw\in\mathbb{C}^nwCn be the left eigenvector normalized so that wv=1w^\ast v = 1wv=1 (here ww^\astw is conjugate transpose). Then the first Lyapunov coefficient 1\ell_11 is (one common expression --- see Kuznetsov / Guckenheimer & Holmes):

1=120Re(wC(v,v,v)2wB(v,(A2i0I)1B(v,v))+wB(v,(A0I)1B(v,v))),\ell_1 \;=\; \dfrac{1}{2\omega_0} \operatorname{Re}\Big( w^\ast C(v,\bar v, v) - 2 w^\ast B\big(v, (A - 2 i\omega_0 I)^{-1} B(v,v)\big) + w^\ast B\big(\bar v, (A - 0\cdot I)^{-1} B(v,\bar v)\big) \Big),1=201Re(wC(v,v,v)2wB(v,(A2i0I)1B(v,v))+wB(v,(A0I)1B(v,v))),

where v\bar vv is the complex conjugate of vvv, and (AI)1(A - \lambda I)^{-1}(AI)1 denotes the resolvent on the subspace complementary to the center (in practice one solves appropriate linear systems for the vectors appearing).

This expression may look forbidding; operationally we do:

Compute vvv and www for JhJ_hJh at h\mu_hh, normalize wv=1w^\ast v = 1wv=1.
Compute B(v,v)B(v,v)B(v,v), then solve (A2i0I)q=B(v,v)(A - 2 i\omega_0 I)q = B(v,v)(A2i0I)q=B(v,v) for qqq (unique provided 2i02i\omega_02i0 is not an eigenvalue).
Compute B(v,v)B(v,\bar v)B(v,v), then solve Ar=B(v,v)A r = B(v,\bar v)Ar=B(v,v) for rrr (unique as 0 is not an eigenvalue on the complement).
Compute the inner products wC(v,v,v)w^\ast C(v,\bar v,v)wC(v,v,v), wB(v,q)w^\ast B(v,q)wB(v,q), wB(v,r)w^\ast B(\bar v,r)wB(v,r).
Assemble 1\ell_11 with the formula above.
If Re(1)<0\mathrm{Re}(\ell_1) <0Re(1)<0 supercritical Hopf (stable small-amplitude limit cycle emerges). If Re(1)>0\mathrm{Re}(\ell_1) >0Re(1)>0 subcritical (unstable cycle, dangerous: possible sudden large amplitude oscillations).

4. Practical simplification for our model

The vector field terms are low-order polynomials/sigmoids; you can compute D2FD^2FD2F and D3FD^3FD3F symbolically. For example, the logistic term in P\dot PP supplies quadratic nonlinearity, SE(E)S_E(E)SE(E) supplies nonlinear terms in EEE. The explicit derivatives (partial derivatives up to third order) are straightforward to compute symbolically once you fix the functional forms (we used logistic / Hill / tanh\tanhtanh earlier).
Because our primary control parameter \mu enters E\dot EE linearly, the transversality derivative F\partial_\mu FF is simple (unit vector in EEE-direction). This simplifies the evaluation of dd\dfrac{d\lambda}{d\mu}dd.
If we work with the reduced 3D core (T,E,P)(T,E,P)(T,E,P), the Jacobian is 333\times333; numerically computing eigenvectors and solving the small linear systems for q,rq,rq,r is routine and stable.
c. Suggested workflow to produce explicit analytic expressions & numbers

1. Choose reduced order --- recommended: start with n=3n=3n=3 (T,E,P) because it captures the essential bifurcation (trust, economy, protest) and is algebraically tractable. Keep HHH as a slave variable if desired (or include it later for two-stage analysis). Document the timescale assumptions if you eliminate variables.
2. Compute fixed points analytically (if possible) or solve numerically for X(,)X^*(\mu,\Theta)X(,). For analytic expressions, you may need to adopt simplifying approximations (e.g., small PPP, linearize SES_ESE near threshold).
3. Form Jacobian J(X;,)J(X^*;\mu,\Theta)J(X;,) and identify c\mu_cc where det(J)\det(J)det(J) or appropriate characteristic polynomial meets the bifurcation condition (zero eigenvalue for fold; pair at i0\pm i\omega_0i0 for Hopf).
4. Compute nullvectors / eigenvectors v,wv,wv,w at the bifurcation point.
5. Evaluate multilinear forms D2F,D3FD^2F, D^3FD2F,D3F at XcX_cXc. These are small arrays of partial derivatives; for a 333-dimensional system you need up to 27 third-order partials (but many are zero for standard terms).
6. Compute normal form coefficients:
For fold: compute a=12wD2F[v,v]a = \tfrac12 w^\top D^2F[v,v]a=21wD2F[v,v] and b=wFb = w^\top \partial_\mu Fb=wF.
For Hopf: compute 1\ell_11 using the algorithm above.
7. Classify the bifurcation using signs of aba bab (fold) and Re1\mathrm{Re}\,\ell_1Re1 (Hopf). Produce normal form approximations and asymptotic estimates for amplitude and period (Hopf) and for the location of saddle-node folds.
8. Produce bifurcation diagrams (P steady states vs \mu) using numerical continuation (AUTO, MATCONT or custom continuation) --- this validates normal-form approximations and shows global structure.
9. Interpretation: express crit\mu_{\mathrm{crit}}crit implicitly as a function of \Theta by locating \mu such that the fold condition holds. Where possible, perform local asymptotic expansion to show crit()=0+iii+O(2)\mu_{\mathrm{crit}}(\Theta) = \mu_{0} + \sum_i \beta_i \Theta_i + O(\|\Theta\|^2)crit()=0+iii+O(2) --- the i\beta_ii are interpretable sensitivities.
d. Worked symbolic template

We are now:

1. Take the explicit reduced T,E,PT,E,PT,E,P deterministic skeleton used earlier:
{T=TSE(E)(L+NN+MM)(1T)(T+T(1R))P,E=EE(G,M)E+PP,P=PP(1P/Pmax)+ESE(E)TTKK(C,Ec,R)P,\begin{cases} \dot T = \alpha_T S_E(E)(L+\sigma_N N+\sigma_M M)(1-T) - (\beta_T+\gamma_T(1-R))P,\\[4pt] \dot E = \mu_E - \lambda_E(G,M) E + \phi_P P,\\[4pt] \dot P = \rho_P P(1-P/P_{\max}) + \eta_E S_E(E) - \delta_T T - \kappa_K K(C,E_c,R)\,P, \end{cases}T=TSE(E)(L+NN+MM)(1T)(T+T(1R))P,E=EE(G,M)E+PP,P=PP(1P/Pmax)+ESE(E)TTKK(C,Ec,R)P,
with chosen concrete smooth forms E(E)=E/(E+)S_E(E)=E^{\kappa}/(E^\kappa+\theta^\kappa)SE(E)=E/(E+) and K=tanh(K())K=\tanh(\rho_K(\cdot))K=tanh(K()).
2. Symbolically compute:
Fixed point equations F(X;)=0F(X;\mu)=0F(X;)=0 and solve numerically for X()X^*(\mu)X().
Jacobian JJJ and compute condition for fold (detJ=0\det J=0detJ=0 with one zero eigenvalue) to locate c\mu_cc numerically for each leadership profile \Theta.
Compute v,wv,wv,w, then aaa and bbb per the formulas given.
Compute Hopf candidate points by solving Re()=0\mathrm{Re}(\lambda)=0Re()=0 and Im()=\mathrm{Im}(\lambda)=\pm\omegaIm()= and then compute 1\ell_11.
3. Produce:
Normal form coefficients a,b,1a,b,\ell_1a,b,1 as explicit numerical values for Soeharto, SBY, Jokowi and Prabowo (using the parameter mappings we already discussed).
Short analytic expansion for crit()\mu_{\mathrm{crit}}(\Theta)crit() (via sensitivity calculation) and plots validating the normal form.

we ran a fold (saddle-node) analysis on the reduced model for the four leadership profiles and produced numerical estimates of the candidate fold points plus the fold normal-form coefficients.

III. Analytical Derivation

A. Fold (saddle-node) analysis --- Numerical results & normal-form extraction

3.A.1 Overview and numerical method

We numerically investigated saddle-node (fold) bifurcations of the reduced three-dimensional model (T,E,P)(T,E,P)(T,E,P) introduced in Section 2. The primary bifurcation parameter is the exogenous economic shock magnitude E\mu_EE. For each leadership profile \Theta (Soeharto, SBY, Jokowi, Prabowo), we performed the following computational steps:

1. Steady-state computation. For a fixed E\mu_EE and leadership vector \Theta, the deterministic vector field F(T,E,P;E,)F(T,E,P;\mu_E,\Theta)F(T,E,P;E,) was integrated forward until approximate convergence to a steady state X(E,)X^*(\mu_E,\Theta)X(E,). Time-integration used an explicit Euler step with sufficiently small step size and long horizon to ensure convergence inside the invariant cube [0,1]3[0,1]^3[0,1]3.
2. Jacobian evaluation. At the numerically found steady state XX^*X we evaluated the Jacobian J(X;E,)J(X^*;\mu_E,\Theta)J(X;E,) by finite differences. The linear stability spectrum {i(E,)}\{\lambda_i(\mu_E,\Theta)\}{i(E,)} is given by the eigenvalues of JJJ.
3. Fold detection (coarse sweep + bisection). We scanned E[0,1]\mu_E\in[0,1]E[0,1] on a coarse grid to detect sign changes of the smallest real part eigenvalue miniRei\min_i \operatorname{Re}\lambda_iminiRei. Upon detecting an interval containing a sign change we refined the root by bisection to estimate the candidate fold value crit()\mu_{\mathrm{crit}}(\Theta)crit() where miniRei0\min_i \operatorname{Re}\lambda_i\approx 0miniRei0.
4. Normal-form coefficient computation. At the candidate point (Xc,c)(X_c,\mu_c)(Xc,c) we computed the right nullvector vvv and left nullvector www associated with the near-zero eigenvalue (via eigenvectors of JJJ and JJ^\topJ). Approximating the directional second derivative D2F(Xc)[v,v]D^2F(X_c)[v,v]D2F(Xc)[v,v] by symmetric finite differences, we computed the fold normal-form coefficients
a=12wD2F(Xc)[v,v],b=wEF(Xc;c),a \;=\; \tfrac{1}{2}\,w^\top D^2F(X_c)[v,v],\qquad b \;=\; w^\top \partial_{\mu_E} F(X_c;\mu_c),a=21wD2F(Xc)[v,v],b=wEF(Xc;c),
where EF\partial_{\mu_E}FEF reduces to the unit vector in the EEE-direction because E\mu_EE enters the model additively in E\dot EE. The numerical values of aaa and bbb verify the fold nondegeneracy conditions a0a\neq 0a=0 and b0b\neq 0b=0.
Implementation notes. The finite-difference computations and eigenvector solves were performed with double precision; step sizes and perturbations for numerical derivatives were chosen conservatively to balance truncation and round-off error. All computations used the coefficient mapping from the seven leadership parameters described in Section 2.C. The full code and interactive outputs are archived in the analysis notebook.

3.A.2 Numerical results (summary)

Table 1. Fold detection and normal-form coefficients (coarse numerical results).

The interactive table produced in the analysis session ("Fold analysis results (coarse)") lists, for each leader, the estimated crit\mu_{\mathrm{crit}}crit, the steady-state coordinates (Tc,Ec,Pc)(T_c,E_c,P_c)(Tc,Ec,Pc) at the bifurcation point, the approximate zero eigenvalue 0\lambda_00, and the normal-form coefficients aaa and bbb. (The full numeric table from the session is included as Table 1 in the manuscript submission files.)

Qualitative summary of the main findings (from the computed table and diagnostic plot):

SBY --- highest estimated crit\mu_{\mathrm{crit}}crit. Under the chosen coefficient mapping, SBY's leadership vector (high legitimacy LLL, high crisis management MMM) produces the most resilient system: larger exogenous economic stress is required to drive the system to fold tipping. In normal-form terms, the associated b>0b>0b>0 and a0a\neq0a=0 satisfy the nondegeneracy conditions; the fold is classical.

Jokowi --- crit\mu_{\mathrm{crit}}crit slightly lower than SBY but still large. Strong narrative control NNN and substantial elite management EcE_cEc increase stabilizing coefficients, shifting the critical shock to the right; however, the slope of the branch past the fold is relatively steep (sharp tipping once the threshold is passed).
Soeharto --- intermediate crit\mu_{\mathrm{crit}}crit. Although Soeharto's high consensus/cooptation produces a nontrivial cooptation term KKK that suppresses protests at low shocks, low legitimacy LLL decreases trust recovery capacity: the fold occurs at moderate E\mu_EE, and the post-fold jump in PPP is abrupt (classic fold / tipping behavior).
Prabowo --- lowest crit\mu_{\mathrm{crit}}crit among the four. The combination of relatively low legitimacy LLL, weaker crisis management MMM, and repressive tilt (low RRR) yields a fragile configuration: even moderate E\mu_EE drives the system through the fold. Normal-form coefficients show aaa and bbb of magnitudes consistent with a robust (nondegenerate) saddle-node.
Diagnostic figure. Figure 1 (displayed in the analysis session) plots miniRei(E)\min_i \operatorname{Re}\lambda_i(\mu_E)miniRei(E) for each leadership profile. The zero-crossings correspond to the crit\mu_{\mathrm{crit}}crit entries in Table 1 and visually illustrate the relative positions of the fold points.

3.A.3 Local normal-form approximation and validation

At each computed fold point we constructed the one-dimensional normal form

x=b(Ec)+ax2,\dot x \;=\; b(\mu_E-\mu_c) + a\, x^2,x=b(Ec)+ax2,

with aaa and bbb taken from the finite-difference computations. To validate the normal form approximation:

We integrated the reduced normal form for parameter perturbations E=c\mu_E=\mu_c\pm\varepsilonE=c with 1\varepsilon\ll11 and compared the equilibrium and transient behavior of the scalar model to trajectories of the full 333-dimensional model initialized near XcX_cXc.
1. The normal form reproduces the leading-order scaling of the post-fold jump and the qualitative behavior (creation/annihilation of equilibria). Figure 2 (constructed in the session) overlays the normal-form bifurcation skeleton with the numerically computed P(E)P^*(\mu_E)P(E) branch for each leader; agreement is excellent close to c\mu_cc, with departures at larger Ec|\mu_E-\mu_c|Ec due to higher-order terms and global nonlinearities.
2. Interpretation. The sign of aaa and the sign of bbb together determine whether equilibria are created or annihilated as E\mu_EE is varied upward. For all four leadership profiles the product aba bab had a sign consistent with a standard saddle-node unlocking unrest at E>c\mu_E>\mu_cE>c in our coefficient mapping. This confirms the qualitative mechanism that economic shocks above a leader-dependent critical value induce abrupt transitions from low to high protest intensity.

3.A.4 Robustness checks and sensitivity

We performed sensitivity checks to verify that fold locations shift coherently with changes in specific leadership parameters:

Legitimacy LLL: increasing LLL by +0.1 shifts crit\mu_{\mathrm{crit}}crit to larger values for all leaders (system becomes more resilient). The magnitude of the shift is largest for profiles with initially low LLL (e.g., Prabowo).
Elite management EcE_cEc: increasing EcE_cEc raises KKK and increases the effective stability; in many cases this shifts crit\mu_{\mathrm{crit}}crit to the right, but at the cost of making post-fold jumps sharper if RRR is low (cooptation + repression tradeoff).
Repression balance RRR: lowering RRR (more repressive) can transiently increase short-term suppression (raising apparent local stability) but increases backfire effects (via T(1R)\gamma_T(1-R)T(1R)), reducing trust and thus lowering crit\mu_{\mathrm{crit}}crit in the medium term --- a quantitatively important trade-off revealed by the fold computations.
These sensitivity tests confirm that the leadership parameter vector \Theta acts as a secondary control set that shifts the primary bifurcation parameter threshold crit=()\mu_{\mathrm{crit}}=\Phi(\Theta)crit=().

3.A.5 Limitations and computational caveats

The analysis used finite differences for second derivatives; publication-grade symbolic derivatives and numerical continuation (AUTO / MATCONT) are recommended to refine bifurcation curves and produce rigorous continuation branches.
Results are conditional on the chosen coefficient mappings from leadership parameters to model coefficients. Empirical calibration (likelihood / Bayesian fitting) would reduce subjectivity.
The reduced T,E,PT,E,PT,E,P model captures the essential fold mechanism but omits explicit dynamical coupling with HHH (black-horse potential). Section 4 provides full simulations including HHH.
3.A.6 Concluding remarks for the fold analysis

The fold analysis demonstrates that leadership quality, encoded in the seven-dimensional \Theta, systematically modulates the critical economic shock amplitude crit\mu_{\mathrm{crit}}crit required for abrupt social unrest. The four case studies illustrate a clear ordering of resilience (SBY Jokowi Soeharto Prabowo) under the adopted mappings. The extracted normal-form coefficients validate that the transitions are standard saddle-node bifurcations and provide interpretable sensitivity measures for policy: increasing legitimacy and crisis management robustness yields the largest shifts in crit\mu_{\mathrm{crit}}crit, raising the regime's resilience to economic perturbations.

B Hopf normal-form calculations and limit-cycle analysis (augmented model)

3.B.1 Why augment the model?

In the reduced (T,E,P)(T,E,P)(T,E,P) system the numerics showed fold (saddle-node) tipping as the dominant route to unrest; a coarse Hopf search returned no complex-pair crossings. To allow oscillatory protest waves (small-amplitude cycles born at Hopf), we introduce an explicit policy/response lag that closes a negative-feedback loop with delay-like dynamics---well known to generate Hopf when gain and lag exceed damping.

We add a slow state GGG (aggregate policy response), yielding a 4D ODE that preserves our leadership linkages:

T=TS(E)(L+NN+MM)(1T)(T+T(1R))P,E=EdEE+PPTTGG,P=PP(1PPmax)+ES(E)TTc4PGG,G=G+kpPktT.\begin{aligned} \dot T &= \alpha_T\,S(E)\,\big(L+\sigma_N N+\sigma_M M\big)\,(1-T)\;-\;(\beta_T+\gamma_T (1\!-\!R))\,P,\\ \dot E &= \mu_E - d_E E + \phi_P P - \phi_T T \;\;-\; \phi_G G,\\ \dot P &= \rho_P P\Big(1-\tfrac{P}{P_{\max}}\Big) + \eta_E S(E) - \delta_T T - c_4 P \;\;-\; \chi_G G,\\ \dot G &= \frac{-G + k_p P - k_t T}{\tau}. \end{aligned}TEPG=TS(E)(L+NN+MM)(1T)(T+T(1R))P,=EdEE+PPTTGG,=PP(1PmaxP)+ES(E)TTc4PGG,=G+kpPktT.

Here S(E)=EE+S(E)=\dfrac{E^\kappa}{E^\kappa+\theta^\kappa}S(E)=E+E is the same stress activation as before. The leadership parameters =(C,L,M,N,ES,EM,R)\Theta=(C,L,M,N,ES,EM,R)=(C,L,M,N,ES,EM,R) modulate coefficients as previously defined; new couplings obey:

Gain: kpk_p\uparrowkp with M,EMM, EMM,EM; ktk_t\uparrowkt with M,LM, LM,L.
Lag: \tau\uparrow when CC\downarrowC and MM\downarrowM (slower, more inertial state reaction).
Policy impacts: G,G\phi_G,\chi_G\uparrowG,G with LLL and ESESES (legitimate policies help both economy and de-escalation).
Linear Hopf conditions (cubic sub-block)

Linearizing at an equilibrium X=(T,E,P,G)X^*=(T^*,E^*,P^*,G^*)X=(T,E,P,G) gives Jacobian JJJ. Because GGG only couples with (T,E,P)(T,E,P)(T,E,P), the oscillatory onset is governed by a cubic characteristic polynomial on the (E,P,G)(E,P,G)(E,P,G) feedback loop (trust acts as a stabilizing leak and coupling):

3+a12+a2+a3=0,\lambda^3 + a_1 \lambda^2 + a_2 \lambda + a_3 =0,3+a12+a2+a3=0,

with a1=trJEPGa_1 = -\operatorname{tr}J_{EPG}a1=trJEPG, a2a_2a2 the sum of principal minors of JEPGJ_{EPG}JEPG, and a3=detJEPGa_3=\det J_{EPG}a3=detJEPG.
A (generic) Hopf occurs when the Routh--Hurwitz equalities are met:

a1>0,a2>0,a3>0,anda1a2=a3,a_1>0,\quad a_2>0,\quad a_3>0,\quad\text{and}\quad a_1 a_2 = a_3,a1>0,a2>0,a3>0,anda1a2=a3,

with ddERe=h0 \tfrac{d}{d\mu_E}\operatorname{Re}\lambda\big|_{\mu=\mu_h} \neq 0dEdRe=h=0. The first Lyapunov coefficient 1\ell_11 (computed from multilinear terms) sets the criticality: 1<0\ell_1<0\Rightarrow1<0 supercritical (small stable cycles); 1>0\ell_1>0\Rightarrow1>0 subcritical (unstable cycles and hysteresis).

3.B.2 Numerical search and outcome (what we found)

I implemented the 4D model and performed a coarse Hopf scan over E[0,1]\mu_E\in[0,1]E[0,1] for the four leadership profiles (Soeharto, SBY, Jokowi, Prabowo), then attempted 1\ell_11 estimation via finite-difference multilinear forms. With conservative (empirically plausible) gains and lags tied to leadership:

No Hopf crossings were detected in the scanned range; the real parts of the most oscillatory eigenpairs remained negative where imag parts were nonzero.
Interpretation: under the baseline mapping, the feedback loop (PGE,P)(P \to G \to E,P)(PGE,P) is not quite "hot" enough (insufficient gain and/or lag relative to damping c4,dEc_4, d_Ec4,dE) to overturn fold-dominated dynamics.
This result does not rule out Hopf; it diagnoses that the current gain--lag--damping combination sits on the stable side of the Routh--Hurwitz boundary.

3.B.3 How to make Hopf appear (testable parameter windows)

From the linear conditions and inspection of JJJ, three levers move a1a2a_1 a_2a1a2 toward a3a_3a3:

1. Increase effective lag or inertia
Raise \tau (slower policy) or insert an additional first-order lag in the GGG channel. Politically: sluggish, procedural or fragmented response.
2. Raise loop gain
Increase G\phi_GG and G\chi_GG (policy acts strongly on EEE and PPP), and/or kpk_pkp (policy reacts aggressively to PPP). High-gain + delay = overshoot oscillations.
3. Reduce intrinsic damping
Lower c4c_4c4 (less protest damping via cooptation/repression blend) or dEd_EdE (slower economic dissipation).
A practical continuation plan for the manuscript:

Treat :=(GGkp)\zeta:=\tau\cdot (\phi_G\,\chi_G\,k_p):=(GGkp) as a Hopf driver. Continue equilibria in (E,)(\mu_E,\zeta)(E,) and track the Hopf curve H()\mathcal{H}(\Theta)H() where a1a2=a3a_1 a_2=a_3a1a2=a3.
Expect Hopf to emerge for high \zeta (slow + strong policy loop) and diminish for large ktk_tkt (trust-sensitive policy reducing overshoot).
Policy reading: fast, trust-weighted responses (high ktk_tkt, moderate kpk_pkp, moderate G,G\phi_G,\chi_GG,G, small \tau) suppress protest waves; slow, over-reactive responses amplify them.

3.B.4 Normal form and interpretation

Once a Hopf point (X,h)(X^*,\mu_h)(X,h) is located, the center-manifold reduction yields

z=(+i0)z+1zz2+O(z4),\dot z = (\sigma + i \omega_0) z + \ell_1 z |z|^2 + \mathcal{O}(|z|^4),z=(+i0)z+1zz2+O(z4),

with =(Eh)\sigma=\alpha(\mu_E-\mu_h)=(Eh). The sign of 1\ell_11 determines whether cycles are born stable (supercritical, 1<0\ell_1<01<0, small periodic protest waves) or born unstable (subcritical, 1>0\ell_1>01>0, dangerous hysteresis lobes where finite shocks kick the system onto large excursions). In our baseline scans, the system remained on the pre-Hopf side; if you want, we can increase \zeta systematically and produce (h,0,1)(\mu_h,\omega_0,\ell_1)(h,0,1) tables plus limit-cycle amplitude curves.

3.B.5 Methods --- targeted Hopf search in the augmented model

To test whether oscillatory protest dynamics (small-amplitude limit cycles) could appear in our framework we augmented the reduced model (T,E,P)(T,E,P)(T,E,P) with an explicit policy/response variable GGG. The augmented deterministic system is

T=TS(E)(L+NN+MM)(1T)(T+T(1R))P,E=EdEE+PPTTGG,P=PP(1PPmax)+ES(E)TTc4PGG,G=G+kpPktT,\begin{aligned} \dot T &= \alpha_T S(E)\,(L+\sigma_N N + \sigma_M M)\,(1-T)\;-\;(\beta_T+\gamma_T(1-R))\,P,\\[4pt] \dot E &= \mu_E - d_E E + \phi_P P - \phi_T T - \phi_G G,\\[4pt] \dot P &= \rho_P P\Big(1-\frac{P}{P_{\max}}\Big) + \eta_E S(E) - \delta_T T - c_4 P - \chi_G G,\\[4pt] \dot G &= \dfrac{-G + k_p P - k_t T}{\tau}, \end{aligned}TEPG=TS(E)(L+NN+MM)(1T)(T+T(1R))P,=EdEE+PPTTGG,=PP(1PmaxP)+ES(E)TTc4PGG,=G+kpPktT,

where S(E)S(E)S(E) is the same saturating stress activation used previously. Leadership parameters \Theta modulate the coefficients as in Section 2C; additional GGG-loop coefficients (kp,kt,,G,Gk_p,k_t,\tau,\phi_G,\chi_Gkp,kt,,G,G) are set by plausible mappings from the leadership vector (e.g., stronger crisis management and elite coordination increase kpk_pkp; lower consensus increases \tau).

Because Hopf bifurcations require both loop gain and effective delay/inertia, we carried out a targeted numerical experiment that systematically raised a scalar multiplier zzz (the "gain lag" multiplier) that scales the GGG-loop gain and lifts the effective lag:

For each leader (Soeharto, SBY, Jokowi, Prabowo) we scanned z{1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0}z \in \{1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0\}z{1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0} and E[0,1]\mu_E\in[0,1]E[0,1] on a coarse grid.
For each (z,E)(z,\mu_E)(z,E) we computed the steady state numerically (forward integration to convergence), evaluated the Jacobian, and tested for an eigenvalue pair with near-zero real part and nonzero imaginary part (the Hopf linear condition).
When a candidate crossing was found we refined E\mu_EE locally and attempted to compute the first Lyapunov coefficient 1\ell_11 via finite-difference multilinear forms to assess criticality (super/subcritical).
All numerical derivatives used conservative finite-difference increments; runs were intentionally coarse to explore the parameter space quickly. The code and parameter mapping used are archived in the analysis notebook.

3.B.6 Results --- targeted search outcomes

Summary (numeric run). The targeted sweep (z up to 5) did not detect Hopf bifurcations for any of the four leadership profiles under the coefficient mappings employed. The run summary CSV is available as:

/mnt/data/Targeted_Hopf_detection__increasing_z_.csv

This negative finding indicates that, with the chosen baseline maps from leadership parameters to model coefficients and with moderate increases of the policy loop gain/lag (up to z=5), the linearization of the augmented system remains in the non-oscillatory regime. In dynamical language: the Routh--Hurwitz conditions for cubic E--P--G subloop stability were not violated for the scanned parameter window; the dominant route to instability remained a saddle-node (fold) in E\mu_EE, as shown in Section 3.1.

Interpretation. The absence of Hopf in the scanned region implies one of three possibilities (or their combination):

1. The effective loop gain lag (the product of kp,G,G,k_p,\phi_G,\chi_G,\taukp,G,G,) required to produce Hopf is larger than the range we scanned (i.e., \zeta must exceed the tested maximum).
2. The baseline damping terms (e.g., c4c_4c4 --- protest damping; dEd_EdE --- economic dissipation) are sufficiently large that even large zzz cannot overcome them.
3. Alternative structural elements are needed to produce Hopf, for example: an extra delay (explicit delay differential equation), stronger nonlinear gain, or an additional state that closes a different feedback loop (media contagion, law enforcement dynamics, etc.).
Because the Hopf detection requires precise eigenvalue crossing, coarse scans and finite-difference multilinear calculations can miss narrow windows. Thus the negative result should be interpreted as no Hopf in the explored, plausible parameter wedge, not a universal impossibility.

3.B.7 Example of a Hopf-generating route (recommendations & tested levers)

Analytical inspection of the linearized E--P--G subblock and the Routh--Hurwitz determinant suggests three straightforward ways to create a Hopf regime in model experiments:

Amplify loop gain: increase kpk_pkp, G\phi_GG, and G\chi_GG simultaneously. Politically, this corresponds to a response policy that is both strong and directly couples to protests and the economy.
Lengthen effective lag: increase \tau (slower policy), or introduce an explicit second lagging state. Sluggish feedback produces phase lag, which together with gain yields oscillation.
Reduce damping: decrease c4c_4c4 and/or dEd_EdE so the loop is less attenuated.
Practical numerical bounds (for reproduction): in exploratory trials we found no Hopf up to z=5z=5z=5 where zzz multiplies kpk_pkp and scales \tau (our mapping used moderate tau scaling). Empirically, to produce Hopf one would likely need to push zzz significantly larger than 5 or set 1\tau\gtrsim 11 with large gains --- parameter windows that must be justified by domain evidence if used in an applied study.

3.B.8 Figures and table (to include)

Table 3. Hopf search summary (targeted z sweep). Include columns: Leader | z tested | Hopf detected (Y/N) | h\mu_hh (if found) | 1\ell_11 (if computed) | 0\omega_00 (if computed). (The CSV created by the run is saved at /mnt/data/Targeted_Hopf_detection__increasing_z_.csv.)
Figure 3. Stability trace around candidate regions. If/when Hopf is found, plot Re() of the most oscillatory eigenpair vs E\mu_EE for several z values to show the crossing; overlay imaginary parts (0).
Figure 4 (representative) --- Limit cycle beyond Hopf. If Hopf emerges for a chosen (z,h)(z,\mu_h)(z,h), integrate the full deterministic system slightly beyond h\mu_hh to produce (i) time series of P,E,TP,E,TP,E,T, and (ii) phase portrait EEE vs PPP. The code attempted to save such a figure to /mnt/data/limit_cycle_example.png if a Hopf was found and a durable cycle was observed.
Note: in the present run the CSV indicates no Hopf was detected (all found = False for the four leaders), so Figure 3/4 are placeholders until a Hopf is located with a refined search.

3.B.9 Robustness and recommended next steps

For a publication-quality treatment of Hopf regimes we recommend the following immediate steps:

1. Continuation-based search: Use AUTO or MATCONT to continue equilibria and detect Hopf loci in the (E,z)(\mu_E,z)(E,z) plane. Continuation is far more reliable than coarse grid scanning and will produce accurate h(z)\mu_h(z)h(z) curves and 1\ell_11 computations. (I can prepare MATCONT/AUTO input files if you want.)
2. Parameter identification: Bound the plausible ranges of kp,,G,Gk_p,\tau,\phi_G,\chi_Gkp,,G,G using historical policy response times and intensities. Avoid arbitrary large zzz unless politically justified.
3. Symbolic derivatives: Replace finite-difference multilinear forms with symbolic Jacobian/Hessian/Tensor derivatives (SymPy/Autograd) to compute 1\ell_11 robustly.
4. If Hopf is found: produce amplitude vs parameter plots (limit cycle amplitude as function of Eh\mu_E-\mu_hEh) and phase portraits for each leader to interpret oscillation amplitude and frequency differences arising from leadership vectors.

IV. Numerical Simulation

A Fixed-point analysis (analytical derivation)

1. Model (deterministic skeleton, restated)

We work with the reduced deterministic skeleton (same notation as Section 2.C). For clarity I rewrite the ODEs in compact notation:

\begin{aligned} \dot T &= \; \alpha_T\,S_E(E)\,B_L\,(1-T)\;-\;B_P\,P, \tag{1}\\[4pt] \dot E &= \; \mu_E \;-\; \lambda_E\,E \;+\; \phi_P\,P \;-\; \phi_T\,T, \tag{2}\\[4pt] \dot P &= \; \rho_P\,P\Big(1-\dfrac{P}{P_{\max}}\Big) \;+\; \eta_E\,S_E(E) \;-\; \delta_T\,T \;-\; \kappa_K\,K\,P, \tag{3} \end{aligned}

where (to condense notation)

SE(E)=EE+S_E(E)=\dfrac{E^\kappa}{E^\kappa+\theta^\kappa}SE(E)=E+E (Hill function),
BL(L+NN+MM)B_L \equiv (L + \sigma_N N + \sigma_M M)BL(L+NN+MM) (trust recovery baseline),
BP(T+T(1R))B_P \equiv (\beta_T + \gamma_T(1-R))BP(T+T(1R)) (trust erosion from protest),
E=E(G,M)\lambda_E=\lambda_E(G,M)E=E(G,M), K=K(C,Ec,R)\kappa_K=\kappa_K(C,E_c,R)K=K(C,Ec,R), and other coefficients are functions of \Theta as specified in Section 2.B--2.C. All parameters are nonnegative; Pmax1P_{\max}\le1Pmax1.
All variables T,E,P[0,1]T,E,P\in[0,1]T,E,P[0,1].

2. Steady-state equations

Set T=E=P=0\dot T=\dot E=\dot P=0T=E=P=0. Denote steady states by T,E,PT^*,E^*,P^*T,E,P. The steady-state system is:

\begin{aligned} 0 &= \alpha_T\,S_E(E^*)\,B_L\,(1-T^*) \;-\; B_P\,P^*, \tag{4a}\\[4pt] 0 &= \mu_E \;-\; \lambda_E\,E^* \;+\; \phi_P\,P^* \;-\; \phi_T\,T^*, \tag{4b}\\[4pt] 0 &= \rho_P\,P^*\Big(1-\dfrac{P^*}{P_{\max}}\Big) \;+\; \eta_E\,S_E(E^*) \;-\; \delta_T\,T^* \;-\; \kappa_K\,K(E^*,\Theta)\,P^*. \tag{4c} \end{aligned}

K(E,)K(E^*,\Theta)K(E,) denotes the elite-cooptation term evaluated at leadership parameters (the particular functional form for KKK was given in Section 2.C).

3. Algebraic reduction --- express TT^*T and eliminate

Equation (4a) is linear in TT^*T. Solve for TT^*T:

TSE(E)BL(1T)=BPP1T=BPTSE(E)BLP.\alpha_T\,S_E(E^*)\,B_L\,(1-T^*) = B_P P^* \quad\Longrightarrow\quad 1 - T^* = \dfrac{B_P}{\alpha_T\,S_E(E^*)\,B_L}\,P^*.TSE(E)BL(1T)=BPP1T=TSE(E)BLBPP.

Hence

T=1BPTBLSE(E)P(5)\boxed{ \; T^* \;=\; 1 \;-\; \dfrac{B_P}{\alpha_T\,B_L\,S_E(E^*)}\, P^* \; } \tag{5}T=1TBLSE(E)BPP(5)

(Important physical constraint: 0T10\le T^*\le10T1 restricts admissible pairs (E,P)(E^*,P^*)(E,P); in particular the denominator must be positive and the right-hand side must be in [0,1][0,1][0,1].)

Substitute (5) into (4b) and (4c) to eliminate TT^*T.

(a) Substitute into economic steady state (4b)

From (4b):

EEE+PPTT=0.\mu_E - \lambda_E E^* + \phi_P P^* - \phi_T T^* = 0.EEE+PPTT=0.

Using (5):

EEE+PPT(1BPTBLSE(E)P)=0.\mu_E - \lambda_E E^* + \phi_P P^* - \phi_T\left(1 - \dfrac{B_P}{\alpha_T B_L S_E(E^*)} P^*\right) = 0.EEE+PPT(1TBLSE(E)BPP)=0.

Rearrange to isolate E\mu_EE:

E=EEPP+TTBPTBLSE(E)P(6)\boxed{ \; \mu_E \;=\; \lambda_E E^* \;-\; \phi_P P^* \;+\; \phi_T \;-\; \phi_T \dfrac{B_P}{\alpha_T B_L S_E(E^*)} P^* \; } \tag{6}E=EEPP+TTTBLSE(E)BPP(6)

Equation (6) will be our expression for E\mu_EE at equilibrium as a function of (E,P)(E^*,P^*)(E,P) and parameters \Theta. In bifurcation analysis E\mu_EE is the primary control parameter; this equation describes the locus of equilibria in (E,P,E)(E,P,\mu_E)(E,P,E) space.

(b) Substitute into protest steady state (4c)

Plug (5) into (4c):

0=PP(1PPmax)+ESE(E)T(1BPTBLSE(E)P)KK(E,)P.0 = \rho_P P^*\Big(1-\frac{P^*}{P_{\max}}\Big) + \eta_E S_E(E^*) - \delta_T\left(1 - \frac{B_P}{\alpha_T B_L S_E(E^*)} P^*\right) - \kappa_K K(E^*,\Theta) P^*.0=PP(1PmaxP)+ESE(E)T(1TBLSE(E)BPP)KK(E,)P.

Collect terms to make a relationship between PP^*P and EE^*E. Write it as:

PP(1PPmax)intrinsic logisticKK(E,)Pcooptation damping+ESE(E)T+TBPTBLSE(E)P=0.(7)\underbrace{\rho_P P^*\Big(1-\frac{P^*}{P_{\max}}\Big)}_{\text{intrinsic logistic}} \; - \; \underbrace{\kappa_K K(E^*,\Theta) P^*}_{\text{cooptation damping}} \; + \; \eta_E S_E(E^*) \; - \; \delta_T \; + \; \delta_T \dfrac{B_P}{\alpha_T B_L S_E(E^*)} P^* \;=\; 0. \tag{7}intrinsic logisticPP(1PmaxP)cooptation dampingKK(E,)P+ESE(E)T+TTBLSE(E)BPP=0.(7)

Rearrange grouping terms in powers of PP^*P:

0=PPmax(P)2+(PKK+TBPTBLSE(E))P+(ESE(E)T).\begin{aligned} 0 &= -\frac{\rho_P}{P_{\max}}(P^*)^2 \;+\; \left(\rho_P - \kappa_K K + \delta_T \dfrac{B_P}{\alpha_T B_L S_E(E^*)}\right) P^* \;+\; \left(\eta_E S_E(E^*) - \delta_T\right). \end{aligned}0=PmaxP(P)2+(PKK+TTBLSE(E)BP)P+(ESE(E)T).

Define for compactness (all quantities evaluated at EE^*E when relevant):

A2(E)PPmax,A1(E)PKK(E,)+TBPTBLSE(E),A0(E)ESE(E)T.\begin{aligned} A_2(E^*) &\equiv -\frac{\rho_P}{P_{\max}},\\[4pt] A_1(E^*) &\equiv \rho_P - \kappa_K K(E^*,\Theta) + \delta_T \dfrac{B_P}{\alpha_T B_L S_E(E^*)},\\[4pt] A_0(E^*) &\equiv \eta_E S_E(E^*) - \delta_T. \end{aligned}A2(E)A1(E)A0(E)PmaxP,PKK(E,)+TTBLSE(E)BP,ESE(E)T.

Then (7) becomes the scalar quadratic in PP^*P:

A2(E)(P)2+A1(E)P+A0(E)=0.(8)\boxed{ \; A_2(E^*) (P^*)^2 + A_1(E^*) P^* + A_0(E^*) \;=\; 0. \; } \tag{8}A2(E)(P)2+A1(E)P+A0(E)=0.(8)

Interpretation: For each admissible EE^*E, the possible PP^*P satisfy (8). The number of real, nonnegative roots PP^*P of this quadratic (subject to 0PPmax0\le P^*\le P_{\max}0PPmax and that TT^*T from (5) lies in [0,1][0,1][0,1]) determines how many equilibria the full system has for that EE^*E. Multiplicity (double root) occurs when discriminant vanishes --- this is the fold condition.

4. Fold condition (saddle-node) in algebraic form

A saddle-node (fold) occurs when two equilibria coalesce; algebraically this is a double root of (8) with respect to PP^*P for some EE^*E. The double-root condition is:

\begin{aligned} F(P^*,E^*) &\equiv A_2(E^*) (P^*)^2 + A_1(E^*) P^* + A_0(E^*) = 0, \tag{9a}\\[4pt] \frac{\partial F}{\partial P}(P^*,E^*) &\equiv 2 A_2(E^*) P^* + A_1(E^*) = 0. \tag{9b} \end{aligned}

Equations (9a)--(9b) eliminate the PP^*P dependence: solving (9b) gives the candidate double root:

P=A1(E)2A2(E).P^* = -\dfrac{A_1(E^*)}{2 A_2(E^*)}.P=2A2(E)A1(E).

Substitute into (9a) to get a single scalar equation in EE^*E (and parameters \Theta):

A2(E)(A1(E)2A2(E))2+A1(E)(A1(E)2A2(E))+A0(E)=0.A_2(E^*)\left(-\frac{A_1(E^*)}{2A_2(E^*)}\right)^2 + A_1(E^*)\left(-\frac{A_1(E^*)}{2A_2(E^*)}\right) + A_0(E^*) = 0.A2(E)(2A2(E)A1(E))2+A1(E)(2A2(E)A1(E))+A0(E)=0.

This simplifies to the discriminant condition:

(E)=A1(E)24A2(E)A0(E)=0.(10)\boxed{ \; \Delta(E^*) \;=\; A_1(E^*)^2 \;-\; 4 A_2(E^*) A_0(E^*) \;=\; 0. \; } \tag{10}(E)=A1(E)24A2(E)A0(E)=0.(10)

Thus the fold locus is given by all EE^*E solving (10). For each such EE^*E we can compute PP^*P via (9b) and TT^*T via (5), and then find E\mu_EE from (6). The resulting E\mu_EE is the candidate crit()\mu_{\mathrm{crit}}(\Theta)crit() where the fold occurs.

Remark: Because AiA_iAi depend on SE(E)S_E(E^*)SE(E), K(E)K(E^*)K(E), etc., (10) is generally a nonlinear equation in EE^*E that must be solved numerically or symbolically (if closed forms are available). But it is explicit and amenable to algebraic manipulation (e.g., resultant elimination) if required.

5. Jacobian at equilibrium --- explicit partial derivatives

To analyze local stability and bifurcation nondegeneracy, compute the Jacobian matrix J=F/XJ=\partial F/\partial XJ=F/X evaluated at the equilibrium X=(T,E,P)X^*=(T^*,E^*,P^*)X=(T,E,P).

From (1)--(3) the Jacobian entries Jij=Xi/XjJ_{ij}=\partial \dot X_i/\partial X_jJij=Xi/Xj are:

J=(TTETPTTEEEPETPEPPP)(T,E,P).J = \begin{pmatrix} \partial_T \dot T & \partial_E \dot T & \partial_P \dot T\\[4pt] \partial_T \dot E & \partial_E \dot E & \partial_P \dot E\\[4pt] \partial_T \dot P & \partial_E \dot P & \partial_P \dot P \end{pmatrix}_{(T^*,E^*,P^*)}.J=TTTETPETEEEPPTPEPP(T,E,P).

Compute each entry:

From T=TSE(E)BL(1T)BPP\dot T = \alpha_T S_E(E) B_L (1-T) - B_P PT=TSE(E)BL(1T)BPP:
TT=TSE(E)BL,ET=TBL(1T)SE(E),PT=BP.\begin{aligned} \partial_T \dot T &= -\alpha_T S_E(E^*) B_L,\\[4pt] \partial_E \dot T &= \alpha_T B_L (1-T^*) S_E'(E^*),\\[4pt] \partial_P \dot T &= -B_P. \end{aligned}TTETPT=TSE(E)BL,=TBL(1T)SE(E),=BP.

From E=EEE+PPTT\dot E = \mu_E - \lambda_E E + \phi_P P - \phi_T TE=EEE+PPTT:
TE=T,EE=E,PE=P.\begin{aligned} \partial_T \dot E &= -\phi_T,\\[4pt] \partial_E \dot E &= -\lambda_E,\\[4pt] \partial_P \dot E &= \phi_P. \end{aligned}TEEEPE=T,=E,=P.

From P=PP(1P/Pmax)+ESE(E)TTKK(E)P\dot P = \rho_P P(1-P/P_{\max}) + \eta_E S_E(E) - \delta_T T - \kappa_K K(E) PP=PP(1P/Pmax)+ESE(E)TTKK(E)P:
Compute derivative carefully: logistic term derivative: P[PP(1P/Pmax)]=P(12P/Pmax)\partial_P [\rho_P P(1-P/P_{\max})] = \rho_P(1 - 2P/P_{\max})P[PP(1P/Pmax)]=P(12P/Pmax).

Also P[KK(E)P]=KK(E)\partial_P[-\kappa_K K(E) P] = -\kappa_K K(E)P[KK(E)P]=KK(E) (here we treated K\kappa_KK constant; if K\kappa_KK depends on P via K(E) only E dependent so ok). And E[ESE(E)KK(E)P]=ESE(E)KPK(E)\partial_E [\eta_E S_E(E) - \kappa_K K(E) P] = \eta_E S_E'(E) - \kappa_K P \,K'(E)E[ESE(E)KK(E)P]=ESE(E)KPK(E).

Hence

TP=T,EP=ESE(E)KPK(E),PP=P(12PPmax)KK(E).\begin{aligned} \partial_T \dot P &= -\delta_T,\\[4pt] \partial_E \dot P &= \eta_E S_E'(E^*) \;-\; \kappa_K P^* K'(E^*),\\[4pt] \partial_P \dot P &= \rho_P\Big(1 - \dfrac{2P^*}{P_{\max}}\Big) \;-\; \kappa_K K(E^*). \end{aligned}TPEPPP=T,=ESE(E)KPK(E),=P(1Pmax2P)KK(E).

Collecting all:

J(X)=(TSE(E)BLTBL(1T)SE(E)BPTEPTESE(E)KPK(E)P(12P/Pmax)KK(E)).(11)\boxed{% J(X^*) = \begin{pmatrix} -\alpha_T S_E(E^*) B_L & \alpha_T B_L (1-T^*) S_E'(E^*) & -B_P\\[6pt] -\phi_T & -\lambda_E & \phi_P\\[6pt] -\delta_T & \eta_E S_E'(E^*) - \kappa_K P^* K'(E^*) & \rho_P(1 - 2P^*/P_{\max}) - \kappa_K K(E^*) \end{pmatrix}. } \tag{11}J(X)=TSE(E)BLTTTBL(1T)SE(E)EESE(E)KPK(E)BPPP(12P/Pmax)KK(E).(11)

Here SE(E)=dSEdE=E1/(E+)2S_E'(E) = \dfrac{dS_E}{dE} = \kappa\,\theta^\kappa\,E^{\kappa-1}/(E^\kappa+\theta^\kappa)^2SE(E)=dEdSE=E1/(E+)2 (or the equivalent derivative of the Hill function used). K(E)K'(E)K(E) is derivative of the cooptation function w.r.t. EEE (zero for our original KKK that depends on leadership parameters only --- if KKK depends only on \Theta and not on EEE, then K(E)=0K'(E)=0K(E)=0, simplifying the expression).

6. Bifurcation conditions and nondegeneracy

Fold (saddle-node): algebraically encoded by (9a--9b) (double root in PPP for some EE^*E). Equivalently, at (Xc,c)(X_c,\mu_c)(Xc,c) the Jacobian has a simple zero eigenvalue: detJ(Xc)=0\det J(X_c)=0detJ(Xc)=0 and rankJJJ=2. Nondegeneracy requires the quadratic nonlinearity coefficient a=12wD2F[v,v]0a = \tfrac12 w^\top D^2F[v,v] \neq 0a=21wD2F[v,v]=0 (see Section 2.D) and b=wF0b = w^\top \partial_{\mu}F \neq 0b=wF=0. For our model F=(0,1,0)\partial_{\mu}F = (0,1,0)^\topF=(0,1,0), so b=w2b= w_2b=w2 (the second component of the left nullvector).
Hopf (if it appears): occurs when the Jacobian has a simple complex conjugate pair crossing the imaginary axis. Routh--Hurwitz criteria on the cubic characteristic polynomial (or direct eigenvalue computation) determine the Hopf locus. The first Lyapunov coefficient 1\ell_11 computed from multilinear forms determines criticality.
7. Recipe to compute crit()\mu_{\mathrm{crit}}(\Theta)crit() analytically / numerically

a. Compute Ai(E)A_i(E)Ai(E) symbolic expressions from definitions in Section 3 (requires SE(E)S_E(E)SE(E), K(E)K(E)K(E), and mapping of leadership parameters to coefficients).
b. Solve discriminant equation (10): A1(E)24A2(E)A0(E)=0A_1(E)^2 - 4 A_2(E) A_0(E) = 0A1(E)24A2(E)A0(E)=0 for EE^*E in (0,1)(0,1)(0,1). Each real root EE^*E is a candidate fold height.
c. For each candidate EE^*E, compute P=A1(E)/(2A2(E))P^* = -A_1(E^*)/(2 A_2(E^*))P=A1(E)/(2A2(E)) (must be 0\ge00 and Pmax\le P_{\max}Pmax). Then compute TT^*T from (5) and finally compute E\mu_EE via (6). That E\mu_EE is crit()\mu_{\mathrm{crit}}(\Theta)crit() associated with that fold.
d. Check Jacobian at (T,E,P)(T^*,E^*,P^*)(T,E,P): compute JJJ in (11). Verify zero eigenvalue (within tolerance) and check that nullspace has dimension 1 (simple eigenvalue). Compute the nondegeneracy coefficients aaa and bbb as in Section 2.D to ensure classical saddle-node.
e. If symbolic closed forms are desired: compute resultant eliminating PP^*P between (8) and its derivative (9b) to obtain a polynomial in EE^*E alone (this is exactly the discriminant), then solve symbolically if possible (SymPy). Otherwise solve numerically with reliable root solvers and continuation.
8. Additional remarks (interpretation & practicalities)

The reduction to a quadratic (8) is central: it demonstrates why fold bifurcations are generic here --- the protest equation includes a logistic (quadratic) nonlinearity and linear couplings to TTT and EEE, so multiplicity of roots (coexisting low- and high-PPP equilibria) is expected. The leadership parameters enter A1,A0A_1,A_0A1,A0 through KKK, SES_ESE, T\delta_TT, etc., and thereby modulate the discriminant (E)\Delta(E)(E). In particular, increasing legitimacy LLL generally increases BLB_LBL, decreasing the TBP/(TBLSE) \delta_T B_P /(\alpha_T B_L S_E)TBP/(TBLSE) term and thereby shifting the discriminant to favor single-root (more stable) regimes.
When KKK does not depend on EEE (pure leadership control), A1A_1A1 simplifies and the discriminant becomes easier; if KKK depends on EEE (e.g., elite cooptation weakens with economic stress), that coupling can create richer algebraic structure (multiple fold points).
For the full 4D augmented model (including GGG), a similar elimination procedure can be attempted but will produce cubic or quartic algebraic relations requiring higher-order resultants. Practically, use continuation tools to track folds and Hopf curves in parameter space.

B Jacobian derivation and eigenvalue conditions

1. Notation and setup

We consider the reduced dynamical system (restating for convenience)

T=F1(T,E,P;,E),E=F2(T,E,P;,E),P=F3(T,E,P;,E),\begin{aligned} \dot T &= F_1(T,E,P;\Theta,\mu_E),\\ \dot E &= F_2(T,E,P;\Theta,\mu_E),\\ \dot P &= F_3(T,E,P;\Theta,\mu_E), \end{aligned}TEP=F1(T,E,P;,E),=F2(T,E,P;,E),=F3(T,E,P;,E),

where \Theta denotes the vector of seven leadership parameters (Consensus CCC, Legitimacy LLL, Crisis management MMM, Narrative NNN, Economic stability ESESES, Elite coordination EMEMEM, Repression vs consensus RRR) and E\mu_EE is the exogenous economic stress (control parameter). Let X=(T,E,P)X=(T,E,P)^\topX=(T,E,P) and write F(X;,E)F(X;\Theta,\mu_E)F(X;,E) for the vector field.

An equilibrium satisfies F(X;,E)=0F(X^*;\Theta,\mu_E)=0F(X;,E)=0. Local linear stability is governed by the Jacobian matrix evaluated at the equilibrium:

J(X;,E)=DXF(X;,E)(the 33 matrix of partial derivatives).J(X^*;\Theta,\mu_E) \;=\; D_X F(X^*;\Theta,\mu_E) \qquad\text{(the 33 matrix of partial derivatives).}J(X;,E)=DXF(X;,E)(the 33 matrix of partial derivatives).

Below I give explicit expressions for JJJ in terms of the model functions used in Section 2/4.A, then derive eigenvalue conditions and the bifurcation criteria.

2. Jacobian --- explicit partial derivatives

From the model in Section 4.A (notation consistent with that section),

F1(T,E,P)=TSE(E)BL(1T)BPP,F2(T,E,P)=EEE+PPTT,F3(T,E,P)=PP(1PPmax)+ESE(E)TTKK(E,)P.\begin{aligned} F_1(T,E,P) &= \alpha_T S_E(E)\,B_L\,(1-T) - B_P P,\\ F_2(T,E,P) &= \mu_E - \lambda_E E + \phi_P P - \phi_T T,\\ F_3(T,E,P) &= \rho_P P\Big(1-\tfrac{P}{P_{\max}}\Big) + \eta_E S_E(E) - \delta_T T - \kappa_K K(E,\Theta) P. \end{aligned}F1(T,E,P)F2(T,E,P)F3(T,E,P)=TSE(E)BL(1T)BPP,=EEE+PPTT,=PP(1PmaxP)+ESE(E)TTKK(E,)P.

(If your exact expressions for some coefficients differ, replace them symbolically; the structure below is generic.)

Differentiate to obtain the Jacobian entries Jij=Fi/xjJ_{ij} = \partial F_i/\partial x_jJij=Fi/xj evaluated at X=(T,E,P)X^*=(T^*,E^*,P^*)X=(T,E,P).

Calculate the (analytical) partials:

SE(E)S_E(E)SE(E) derivative:
SE(E)=ddEEE+=E1(E+)2.S_E'(E) \;=\; \frac{d}{dE}\frac{E^\kappa}{E^\kappa + \theta^\kappa} \;=\; \frac{\kappa\,\theta^\kappa\,E^{\kappa-1}}{(E^\kappa+\theta^\kappa)^2}.SE(E)=dEdE+E=(E+)2E1.

Entries:
J11=TF1=TSE(E)BL,J12=EF1=TBL(1T)SE(E),J13=PF1=BP,\begin{aligned} J_{11} &= \partial_T F_1 = -\alpha_T\,S_E(E^*)\,B_L,\\[4pt] J_{12} &= \partial_E F_1 = \alpha_T\,B_L\,(1-T^*)\,S_E'(E^*),\\[4pt] J_{13} &= \partial_P F_1 = -B_P, \end{aligned}J11J12J13=TF1=TSE(E)BL,=EF1=TBL(1T)SE(E),=PF1=BP, J21=TF2=T,J22=EF2=E,J23=PF2=P,\begin{aligned} J_{21} &= \partial_T F_2 = -\phi_T,\\[4pt] J_{22} &= \partial_E F_2 = -\lambda_E,\\[4pt] J_{23} &= \partial_P F_2 = \phi_P, \end{aligned}J21J22J23=TF2=T,=EF2=E,=PF2=P, J31=TF3=T,J32=EF3=ESE(E)KPK(E,),J33=PF3=P(12PPmax)KK(E,).\begin{aligned} J_{31} &= \partial_T F_3 = -\delta_T,\\[4pt] J_{32} &= \partial_E F_3 = \eta_E S_E'(E^*) - \kappa_K P^* \, K'(E^*,\Theta),\\[4pt] J_{33} &= \partial_P F_3 = \rho_P\Big(1 - 2\frac{P^*}{P_{\max}}\Big) - \kappa_K K(E^*,\Theta). \end{aligned}J31J32J33=TF3=T,=EF3=ESE(E)KPK(E,),=PF3=P(12PmaxP)KK(E,).

Notes:

If KKK (elite cooptation) does not depend on EEE in your parameterization, then K(E)=0K'(E^*)=0K(E)=0 and J32=ESE(E)J_{32}=\eta_E S_E'(E^*)J32=ESE(E).
All coefficient symbols (e.g. T,E,P,E,T,K,BL,BP\alpha_T,\lambda_E,\rho_P,\eta_E,\delta_T,\kappa_K,B_L,B_PT,E,P,E,T,K,BL,BP) are functions of the leadership vector \Theta via the mappings in Section 2.B--2.C; therefore JJJ is an explicit function of \Theta and XX^*X.
Compactly:

J(X)=(J11J12J13J21J22J23J31J32J33).J(X^*) = \begin{pmatrix} J_{11} & J_{12} & J_{13}\\[4pt] J_{21} & J_{22} & J_{23}\\[4pt] J_{31} & J_{32} & J_{33} \end{pmatrix}.J(X)=J11J21J31J12J22J32J13J23J33.

3. Characteristic polynomial and stability (Routh--Hurwitz)

The eigenvalues \lambda satisfy the characteristic equation

det(IJ)=0.\det(\lambda I - J) \;=\; 0.det(IJ)=0.

Expanding gives a cubic polynomial:

3+c12+c2+c3=0,\lambda^3 + c_1 \lambda^2 + c_2 \lambda + c_3 = 0,3+c12+c2+c3=0,

with the standard coefficients expressed in terms of traces and minors of JJJ:

c1=tr(J)=(J11+J22+J33)c_1 = -\operatorname{tr}(J) = -(J_{11}+J_{22}+J_{33})c1=tr(J)=(J11+J22+J33),
c2=c_2 =c2= sum of principal 222\times222 minors
 c2=J11J22+J11J33+J22J33J12J21J13J31J23J32,c_2 = J_{11}J_{22} + J_{11}J_{33} + J_{22}J_{33} - J_{12}J_{21} - J_{13}J_{31} - J_{23}J_{32},c2=J11J22+J11J33+J22J33J12J21J13J31J23J32,

c3=det(J).c_3 = -\det(J).c3=det(J).
(Conventions sign: here polynomial written 3+c12+c2+c3\lambda^3 + c_1\lambda^2 + c_2\lambda + c_33+c12+c2+c3.)

Linear stability (all eigenvalues have negative real parts) is characterized by the Routh--Hurwitz conditions for a cubic:

c1>0,c3>0,c1c2>c3.\boxed{c_1 > 0,\qquad c_3 > 0,\qquad c_1 c_2 > c_3.}c1>0,c3>0,c1c2>c3.

If any condition fails, the equilibrium loses stability.

4. Bifurcation conditions: saddle-node and Hopf

(a) Saddle-node (fold)

A saddle-node occurs when an eigenvalue crosses through zero: =0\lambda=0=0 is a simple root of the characteristic polynomial. Algebraically this means

detJ(Xc)=0(i.e. c3=0),\det J(X_c) = 0 \quad\text{(i.e. } c_3 = 0\text{),}detJ(Xc)=0(i.e. c3=0),

with the zero eigenvalue being simple (algebraic multiplicity 1, geometric multiplicity 1). Equivalently:

c3=0c_3 = 0c3=0 and c20c_2 \neq 0c2=0 (generically).
Rank condition: rankJ=2\mathrm{rank}\,J = 2rankJ=2 (one dimensional nullspace).
Nondegeneracy / transversality for a saddle-node requires:

Let vvv span the right nullspace: Jv=0J v = 0Jv=0.
Let www span the left nullspace: wJ=0w^\top J = 0wJ=0.
Normalization wv=1w^\top v = 1wv=1.
Define
 a=12wDX2F(Xc)[v,v](quadratic coefficient),b=wEF(Xc),a \;=\; \tfrac{1}{2}\, w^\top D_X^2 F(X_c)[v,v] \quad\text{(quadratic coefficient)},\qquad b \;=\; w^\top \partial_{\mu_E} F(X_c),a=21wDX2F(Xc)[v,v](quadratic coefficient),b=wEF(Xc),
where DX2F[X][u,v]D_X^2 F[X][u,v]DX2F[X][u,v] is the bilinear form of second derivatives (Hessian applied to u,vu,vu,v).
A generic saddle-node requires a0a\neq 0a=0 and b0b\neq 0b=0. In our model EF=(0,1,0)\partial_{\mu_E}F=(0,1,0)^\topEF=(0,1,0) so b=w2b = w_2b=w2 (the second component of the left nullvector).
Interpretation: if a0a\neq0a=0 and b0b\neq0b=0 then locally near the fold the reduced scalar normal form is

y=ay2+b(Ec)+higher orders,\dot y = a y^2 + b (\mu_E-\mu_c) + \text{higher orders},y=ay2+b(Ec)+higher orders,

so the local scaling of the two merging equilibria and their behavior follow classical saddle-node normal form.

(b) Hopf

A Hopf bifurcation requires a pair of complex conjugate eigenvalues to cross the imaginary axis:

At Hopf, 1,2=i0\lambda_{1,2} = \pm i\omega_01,2=i0 (simple pair) and the third eigenvalue has (3)0\Re(\lambda_3) \neq 0(3)=0.
Algebraically at the Hopf point the cubic coefficients satisfy (Hurwitz borderline)
c1c2=c3,c1>0,c2>0,c3>0,c_1 c_2 = c_3,\qquad c_1>0,\quad c_2>0,\quad c_3>0,c1c2=c3,c1>0,c2>0,c3>0,
and the crossing is transversal:
ddE(1,2)E=h0.\left.\frac{d}{d\mu_E}\Re(\lambda_{1,2})\right|_{\mu_E=\mu_h} \neq 0.dEd(1,2)E=h=0.

Nondegeneracy (to determine criticality) requires the first Lyapunov coefficient 1\ell_11 (computed from second and third derivatives) to be nonzero. If 1<0\ell_1<01<0 the Hopf is supercritical (stable small amplitude limit cycles born); if 1>0\ell_1>01>0 it is subcritical (unstable cycles and often hysteresis).
We give the standard formula for 1\ell_11 below.

5. Formulae for normal-form coefficients (fold a,ba,ba,b and Hopf 1\ell_11)

(a) Saddle-node coefficients a,ba,ba,b

Right nullvector vvv and left nullvector www normalized so wv=1w^\top v = 1wv=1.
Quadratic form DX2F(Xc)[u,v]D_X^2 F(X_c)[u,v]DX2F(Xc)[u,v] is defined componentwise:
(DX2F(Xc)[u,v])i=j,k2Fixjxk(Xc)ujvk.\big(D_X^2 F(X_c)[u,v]\big)_i = \sum_{j,k} \frac{\partial^2 F_i}{\partial x_j \partial x_k}(X_c)\,u_j v_k.(DX2F(Xc)[u,v])i=j,kxjxk2Fi(Xc)ujvk.

Then
a=12wDX2F(Xc)[v,v],b=wF(Xc),a = \tfrac12\, w^\top D_X^2 F(X_c)[v,v],\qquad b = w^\top F_{\mu}(X_c),a=21wDX2F(Xc)[v,v],b=wF(Xc),
with F=F/EF_{\mu} = \partial F/\partial \mu_EF=F/E.
(These are exactly the coefficients that appear in the reduced scalar normal form y=ay2+b(c)+ \dot y = a y^2 + b (\mu-\mu_c) + \cdotsy=ay2+b(c)+.)

(b) Hopf first Lyapunov coefficient 1\ell_11 (Kuznetsov formula)

Let J=DXF(Xc)J = D_X F(X_c)J=DXF(Xc) have a simple pair of pure imaginary eigenvalues i0\pm i\omega_0i0 with right eigenvector vvv and left eigenvector www normalized as w,v=wv=1\langle w, v\rangle = w^\dagger v = 1w,v=wv=1. Define the multilinear forms BBB and CCC:

B(u,v)=DX2F(Xc)[u,v]B(u,v) \;=\; D_X^2 F(X_c)[u,v]B(u,v)=DX2F(Xc)[u,v] (vector valued bilinear form),
C(u,v,w)=DX3F(Xc)[u,v,w]C(u,v,w) \;=\; D_X^3 F(X_c)[u,v,w]C(u,v,w)=DX3F(Xc)[u,v,w] (vector valued trilinear form).
Then (see Kuznetsov, Elements of Applied Bifurcation Theory, 3.4) the first Lyapunov coefficient is

1=120{w,C(v,v,v)2w,B(v,(J2i0I)1B(v,v))+w,B(v,(J)1B(v,v))}.\boxed{% \ell_1 \;=\; \frac{1}{2\omega_0}\,\Re\Big\{\,\langle w,\, C(v,v,\bar v)\rangle - 2\langle w,\, B\big(v,\, (J-2i\omega_0 I)^{-1} B(v,v)\big)\rangle + \langle w,\, B\big(\bar v,\, (J)^{-1} B(v,\bar v)\big)\rangle \Big\}.}1=201{w,C(v,v,v)2w,B(v,(J2i0I)1B(v,v))+w,B(v,(J)1B(v,v))}.

Here:

v\bar vv is complex conjugate of vvv.
Inverse operators (J2i0I)1(J-2i\omega_0 I)^{-1}(J2i0I)1 and J1J^{-1}J1 act on the right-hand side vectors; these inverses exist under the generic Hopf nonresonance assumptions.
If 1<0\ell_1<01<0 the Hopf is supercritical.
This formula can be expanded componentwise using the explicit second and third derivatives of FFF. For a 3-D system it is straightforward but algebraically lengthy --- using symbolic algebra packages or finite-difference multilinear approximations is common in practice.

6. Eigenvalue sensitivity and transversality derivatives

For transversality checks you often need the derivative of a particular eigenvalue ()\lambda(\mu)() with respect to the parameter =E\mu=\mu_E=E. For a simple eigenvalue with right eigenvector vvv and left www normalized wv=1w^\top v=1wv=1,

dd=wF=wF(X)\boxed{ \; \frac{d\lambda}{d\mu}\;=\; w^\top \frac{\partial F}{\partial \mu}\;=\; w^\top F_{\mu}(X^*) \; }dd=wF=wF(X)

evaluated at the bifurcation point (since F\partial_\mu FF is the parameter derivative of the vector field). In our case F=(0,1,0)F_\mu=(0,1,0)F=(0,1,0), so dd=w2\dfrac{d\lambda}{d\mu} = w_2dd=w2.

For a Hopf, the transversality condition is dd=h0\left.\dfrac{d}{d\mu}\Re\lambda\right|_{\mu=\mu_h} \neq 0dd=h=0. Numerically you can compute the slope by finite differences of \Re\lambda across \mu or use the eigenvector formula above when www and vvv are available.

7. Practical numerical recipe (steps for reproducible computation)

a. Find equilibrium XX^*X for target \mu and \Theta (Newton + continuation is recommended).
b. Assemble Jacobian JJJ using the analytic partials above (preferred) or high-accuracy finite differences (if analytic derivatives are tedious). Use double precision.
c. Compute eigenpairs of JJJ: eigenvalues i\lambda_ii and eigenvectors viv_ivi; compute left eigenvectors wiw_iwi of JJ^\topJ. Normalize each pair so wivi=1w_i^\dagger v_i = 1wivi=1.
d. Stability test: evaluate c1,c2,c3c_1,c_2,c_3c1,c2,c3 and Routh--Hurwitz inequalities.
e. Fold detection: check detJ0\det J \approx 0detJ0 (or min real eigenvalue 0). If a simple zero eigenvalue is found, compute right/left nullvectors v,wv,wv,w. Compute a=12wD2F[v,v]a = \tfrac12 w^\top D^2F[v,v]a=21wD2F[v,v] and b=wFb = w^\top F_\mub=wF. If aaa and bbb are not too small, you have a generic fold; record c\mu_cc, XcX_cXc, a,ba,ba,b.
Compute D2F[v,v]D^2F[v,v]D2F[v,v] analytically or by central finite differences: D2F[v,v]F(X+v)+F(Xv)2F(X)2D^2F[v,v] \approx \dfrac{F(X+\varepsilon v) + F(X-\varepsilon v) - 2F(X)}{\varepsilon^2}D2F[v,v]2F(X+v)+F(Xv)2F(X) for small \varepsilon. Use \varepsilon around 10610^{-6}106--10410^{-4}104 and test convergence.
f. Hopf detection: check for a complex conjugate pair with 0\Re\lambda \approx 00 and 0=>0\omega_0 = |\Im\lambda|>00=>0. If found, verify transversality via d/dd\Re\lambda/d\mud/d (use eigenvector formula or finite differences). Compute 1\ell_11 using the Kuznetsov formula: compute BBB and CCC multilinear forms (via analytic derivatives or finite differences). Use careful finite-difference parameters (e.g., central differences eps ~ 10510^{-5}105--10410^{-4}104).
g. Sanity checks & sensitivity: perturb \Theta slightly to verify smooth dependences and check that computed a,b,1a,b,\ell_1a,b,1 signs are robust.
8. Recommended implementation details and pitfalls

Prefer analytic derivatives whenever possible. The Jacobian and second/third derivatives can be derived symbolically (SymPy, Mathematica) from the model equations and exported to code. This eliminates finite-difference error accumulation.
If using finite differences, use central differences and test convergence by varying step sizes (e.g., 10610^{-6}106 10410^{-4}104 10210^{-2}102) and check that results converge. For higher derivatives (third order) increase step sizes slightly to avoid round-off noise.
Eigenvector normalization: use complex conjugate left eigenvectors properly; compute left eigenvector www as eigenvector of JJ^\topJ corresponding to the conjugate eigenvalue and normalize wv=1w^\dagger v = 1wv=1.
Near-degenerate cases: if the zero eigenvalue at fold is not simple (rare), the normal form degenerates; numerical indicators are very small denominators when normalizing wvw^\top vwv --- treat cautiously.
Continuation: use pseudo-arclength or AUTO/MATCONT for reliable branch tracing. Continuation significantly reduces numerical sensitivity in locating bifurcation points and enables computation of curves c()\mu_c(\Theta)c().
Validation: once you identify a fold or Hopf, perturb \mu slightly and integrate the full nonlinear ODE to confirm the qualitative behavior (e.g., two equilibria coexisting or small amplitude oscillations). For Hopf, check growth/decay of small perturbations to verify super/subcritical nature predicted by 1\ell_11.
9. How leadership parameters \Theta enter the eigenvalue conditions

Leadership parameters \Theta enter JJJ through the coefficient maps (e.g., BL(),BP(),T(),B_L(\Theta), B_P(\Theta), \delta_T(\Theta),BL(),BP(),T(), etc.). Therefore:

All entries Jij=Jij(X,)J_{ij} = J_{ij}(X^*,\Theta)Jij=Jij(X,). Changes in a leadership coordinate (say LLL or RRR) modify the entries and shift the cubic coefficients c1,c2,c3c_1,c_2,c_3c1,c2,c3, thereby moving eigenvalues in the complex plane.
Sensitivity derivatives of a critical eigenvalue with respect to a leadership parameter pp\in\Thetap are given by
p=wFp,\frac{\partial \lambda}{\partial p} = w^\top \frac{\partial F}{\partial p},p=wpF,
evaluated at the equilibrium (for a simple eigenvalue), where pF\partial_p FpF is the parameter derivative of the vector field (this follows the same eigenvalue perturbation formula as for \mu). This is useful for ranking which leadership attributes most strongly shift the stability boundary c()\mu_c(\Theta)c().
Using these formulas you can compute gradients c()\nabla_\Theta \mu_c(\Theta)c() (sensitivity of the fold location to leadership) via implicit differentiation of the fold normal-form equation ay2+b(c)=0a y^2 + b(\mu-\mu_c) =0ay2+b(c)=0 or via differentiating the algebraic discriminant as in Section 4.A.
Short algorithm 

Solve F(X;,)=0F(X;\Theta,\mu)=0F(X;,)=0 for XX^*X.
Compute analytic Jacobian JJJ from formulas above.
Compute eigenvalues/eigenvectors of JJJ. Check Routh--Hurwitz.
If min()0\min \Re(\lambda)\approx 0min()0 refine \mu by bisection/continuation to locate the critical c\mu_cc.
At critical point compute nullvectors v,wv,wv,w and normal-form coefficients a,ba,ba,b for fold or 1\ell_11 for Hopf.
Validate by direct time simulation (perturbation and integration) to observe expected dynamics.

C. Sensitivity analysis: how legitimacy, narrative, and repression shift _crit

To assess how specific leadership dimensions affect systemic resilience, we performed a parameter sensitivity analysis centered on the critical economic stress threshold crit\mu_{\text{crit}}crit. This threshold marks the onset of instability, where equilibrium stability is lost via a saddle--node or Hopf bifurcation. By perturbing selected leadership parameters around their baseline values, we can quantify both the direction and magnitude of their impact on crit\mu_{\text{crit}}crit.

1. Legitimacy

Legitimacy acts primarily through the stabilizing coefficient T\delta_TT, which links trust erosion to the underlying political mandate. Increasing legitimacy raises T\delta_TT, thereby dampening the erosion of trust in the presence of protests. Numerical sweeps show that higher legitimacy systematically shifts crit\mu_{\text{crit}}crit upward, meaning that a system with strong legitimacy can absorb greater economic shocks before reaching instability. Conversely, low legitimacy sharply lowers crit\mu_{\text{crit}}crit, making unrest more likely even under modest stress.

2. Narrative control

Narrative influences the trust-building pathway BLB_LBL and the damping of protest growth. Strong narrative capacity amplifies the conversion of economic stress into resilience by maintaining coherence in public perception. Sensitivity curves demonstrate that improving narrative control increases the slope of the stability region in the (E,P)(E,P)(E,P) plane, thereby delaying bifurcation. However, unlike legitimacy, narrative has a nonlinear impact: once narrative control surpasses a moderate threshold, marginal improvements yield diminishing returns. Weak narrative, on the other hand, can collapse stability abruptly, lowering crit\mu_{\text{crit}}crit to values near zero.

3. Repression vs. consensus

The repression--consensus balance modifies both the direct protest damping term (c4c_4c4) and the backfire channel feeding protest intensity. Sensitivity experiments reveal a U-shaped effect:

At moderate levels, repression stabilizes the system by containing protests, thus pushing crit\mu_{\text{crit}}crit upward.
At excessive levels, repression generates backlash, reducing consensus and indirectly amplifying protest dynamics. This lowers crit\mu_{\text{crit}}crit, making the system fragile to even small shocks.
Conversely, a tilt toward consensus-building (low repression) broadens the stability basin, but only if supported by sufficient legitimacy. Otherwise, low repression combined with low legitimacy leads to rapid erosion of trust.
4. Comparative influence

When normalized, legitimacy exerts the strongest linear effect on crit\mu_{\text{crit}}crit, followed by repression (nonlinear and regime-dependent), and finally narrative (saturating impact). Taken together, these three parameters provide interpretable levers of systemic resilience. Our formalism thus operationalizes the intuitive claim that leaders with higher legitimacy, credible narratives, and calibrated use of repression extend the stress threshold before unrest becomes inevitable.

D Simulation of black horse emergence (HHH)

An essential feature of our extended model is the variable H(t)H(t)H(t), representing the black horse potential---the latent probability that a new challenger or outsider figure rises when systemic instability creates political openings. Unlike trust (TTT), economic stress (EEE), or protest intensity (PPP), the black horse potential is not directly observable until crystallization occurs. Instead, it can be modeled as an emergent state variable driven by two forces: systemic instability and leadership weakness.

1. Governing mechanism

We extend the reduced system with a fourth equation:

H=PP+EELL(1T)HH,\dot H \;=\; \sigma_P \, P \;+\; \sigma_E \, E \;-\; \sigma_L \, L \,(1-T) \;-\; \gamma_H H,H=PP+EELL(1T)HH,

where:

P\sigma_PP weights protest intensity as a driver of alternative figure emergence;
E\sigma_EE links black horse probability to economic hardship;
L\sigma_LL measures legitimacy as a suppressor of outsider viability;
H\gamma_HH is a decay term reflecting natural dissipation of outsider momentum in stable contexts.
In this representation, HHH rises when protests intensify and the economy weakens, especially under low legitimacy. Conversely, high legitimacy dampens HHH, while stable conditions reduce it over time.

2. Simulation insights

Numerical simulations illustrate how H(t)H(t)H(t) behaves under different leadership profiles:

Soeharto. High legitimacy and repression kept HHH suppressed for decades. Only when E\mu_EE surged during the 1997--98 Asian Financial Crisis did HHH rise abruptly, enabling outsider figures to emerge (e.g., Reformasi leaders).
SBY. Sustained legitimacy and moderate narrative control kept HHH low, even amid shocks. Outsider momentum was weak, and no black horse displaced him.
Jokowi. Narrative saturation lowered protest intensity, but when legitimacy wavered, HHH rose intermittently---manifest in recurrent outsider challenges in elections, though none unseated him during stable phases.
Prabowo. Low legitimacy combined with fragile repression balance elevates baseline HHH. Simulations show that even moderate protests can sustain black horse potential, making the regime vulnerable to outsider competition (e.g., figures like Anies Baswedan or populist celebrities).
3. Interpretation

The emergence of a black horse is thus not random but structurally predictable: it occurs when systemic stress pushes the equilibrium past crit\mu_{\text{crit}}crit, creating a positive feedback loop in HHH. Leadership quality, particularly legitimacy, acts as the main suppressor. Therefore, our simulations suggest that outsider breakthroughs are most likely under regimes that combine low legitimacy with rising economic stress and poorly calibrated repression.

V. Case Study: Indonesia

A. Soeharto: high consolidation, low legitimacy abrupt tipping

Under Soeharto's New Order regime (1966--1998), the system displayed a paradoxical stability. Elite management, institutional consensus, and narrative control were consolidated to near-maximum levels, creating a strong coefficient in our formalism. This yielded high apparent resilience: protests were suppressed effectively, and the black horse potential () remained close to zero for decades. However, legitimacy was chronically low, derived from authoritarian control rather than democratic mandate. In our sensitivity analysis, this situates Soeharto near the upper peak of the repression curve but on the low end of the legitimacy curve.

The consequence was a regime that could withstand routine disturbances but was extremely brittle under extreme stress. When the Asian Financial Crisis of 1997--98 sharply increased economic stress (), the system crossed its bifurcation threshold abruptly. Trust () collapsed, protest intensity () surged beyond repression capacity, and the black horse potential () rapidly emerged, enabling Reformasi leaders to challenge the regime. Mathematically, the trajectory resembled a saddle--node tipping: stability persisted until the system was pushed beyond the fold, at which point equilibrium vanished and regime collapse became irreversible.

Thus, Soeharto's case illustrates how high consolidation without legitimacy can create long-lived but brittle equilibria, where stability is an illusion that masks underlying vulnerability. In our formalism, Soeharto exemplifies an abrupt tipping regime: resilient under moderate shocks, but catastrophically unstable once critical stress is exceeded.

B. SBY: high legitimacy, moderate coalition resilient

Susilo Bambang Yudhoyono's presidency (2004--2014) represents a contrasting configuration to Soeharto's New Order. Whereas Soeharto relied on consolidation and coercion, SBY governed with comparatively high legitimacy derived from democratic elections and sustained public trust. In our formalism, this translates into a high value of the legitimacy parameter, which strongly raises the stability coefficient . As a result, the system's critical threshold was shifted upward, allowing the polity to absorb larger shocks without losing stability.

Coalition management under SBY was moderate rather than absolute. While he did not command the same total consensus as Soeharto or Jokowi, his reliance on institutionalized coalitions and consultative styles created a balanced equilibrium. This configuration sits near the rising linear portion of the legitimacy curve in our sensitivity analysis: resilience grew steadily with legitimacy, even in the absence of extreme repression or elite saturation.

Simulations of the model under SBY-like parameters show that trust () decays slowly under stress, while protest intensity () remains bounded. Importantly, the black horse potential () seldom rises above a negligible baseline, reflecting the lack of viable outsider challengers during his tenure. Economic turbulence, such as the global financial crisis of 2008, did not destabilize the system: protests were limited, and trust levels recovered relatively quickly once external stress eased. This dynamic resembles a system near but safely inside the stability basin, with no bifurcation triggered.

In short, SBY's regime exemplifies resilient equilibrium: high legitimacy provided a strong buffer against economic stress, moderate coalition management preserved systemic flexibility, and the calibrated use of narrative and repression kept black horse emergence suppressed. Within our framework, SBY demonstrates how legitimacy-centered leadership can prevent tipping events even amid significant external shocks.

C. Jokowi: strong narrative, elite cooptation stable until sudden shift

Joko Widodo's presidency (2014--2024) highlights the power and limits of narrative saturation combined with elite cooptation. Unlike Soeharto's coercion-heavy model or SBY's legitimacy-centered approach, Jokowi emphasized broad narrative framing---"infrastructure acceleration" and "down-to-earth populism"---that boosted trust () even during periods of economic stress. In our model, this translates into a high narrative coefficient in , raising the conversion of stress into resilience. Simultaneously, Jokowi's skillful elite cooptation elevated the consensus parameter, pushing the stabilizing term near its maximum.

The result was a regime with an unusually wide basin of attraction: protest intensity () remained suppressed, and black horse potential () stayed low despite structural inequalities and intermittent crises. Numerical simulations show that for moderate increases in , trust erosion is offset by narrative reinforcement, and equilibrium stability is preserved.

However, the sensitivity analysis also reveals diminishing returns. Narrative's effect on stability is saturating: beyond a certain point, additional narrative control no longer significantly raises . If legitimacy declines---due to perceptions of democratic backsliding or overcentralization---then the system becomes vulnerable to an abrupt shift. In simulations, once trust () falls below a threshold, repression and narrative can no longer compensate, and the system tips into instability. At this point, black horse potential () rises rapidly, reflecting outsider figures gaining traction in the political landscape.

Thus, Jokowi's configuration exemplifies a stable-until-sudden-shift regime: stability is prolonged by strong narrative and elite consensus, but fragility accumulates silently when legitimacy wanes. In formal terms, the system hovers near the plateau of the narrative curve, where resilience no longer grows proportionally, leaving it susceptible to sharp bifurcations once critical thresholds are crossed.

D. Prabowo: low legitimacy, high repression low _crit, unstable under shocks

Prabowo Subianto's current configuration illustrates the vulnerabilities of combining fragile legitimacy with strong elite cooptation and coercive tendencies. Unlike Soeharto, who paired repression with long-term institutional embedding, Prabowo's legitimacy remains contested, deriving more from elite bargains than broad-based democratic mandate. In our formalism, this maps to a low legitimacy parameter, which depresses the stabilizing coefficient and lowers the critical stress threshold .

At the same time, Prabowo exhibits strong elite management and a proclivity toward repression. This places him on the descending slope of the repression curve: moderate repression may contain protests in the short term, but excessive reliance on it generates backlash, amplifying protest intensity () and feeding the black horse potential (). Simulations under Prabowo-like parameters show that even moderate increases in economic stress () can destabilize the system, causing trust () to erode rapidly. Once falls below a critical threshold, protests surge and outsider challengers gain immediate viability.

The consequence is a regime prone to early instability: lies much lower than in the cases of Soeharto, SBY, or Jokowi. Where those leaders could absorb considerable stress before tipping, Prabowo's system risks bifurcation under relatively small shocks---whether economic downturns, policy failures, or sudden legitimacy crises. In this sense, the black horse potential is already elevated ex ante, making the emergence of alternative figures not a rare event but an almost structural feature of the political landscape.

In summary, Prabowo's profile exemplifies a fragile equilibrium: elite control and repression appear to sustain stability, but without legitimacy as ballast, the equilibrium is shallow. Our formalism predicts that this configuration is highly sensitive to perturbations, with black horse dynamics activated at lower stress thresholds than in any of the preceding cases.

VI. Discussion

A. Theoretical implications: leadership as mathematical modulator

Our results demonstrate that leadership qualities are not merely normative or descriptive constructs; they can be formalized as mathematical modulators of systemic stability. By embedding seven leadership parameters into a nonlinear dynamical system, we show that attributes such as legitimacy, narrative capacity, coalition management, and repression function as control coefficients that shift bifurcation thresholds. In this sense, leadership acts as a form of "political elasticity," directly modulating how much exogenous stress a system can absorb before tipping into unrest.

This formalism offers three key theoretical implications. First, it bridges the gap between political science and applied mathematics by treating leadership variables as endogenous levers that alter system dynamics. Second, it demonstrates that stability or instability cannot be reduced solely to structural or economic conditions; leadership qualities fundamentally determine whether a society approaches, resists, or surpasses bifurcation points. Third, the model underscores that leadership effects are nonlinear: legitimacy exhibits near-linear stabilization, narrative saturates, and repression has a U-shaped trajectory that can flip from stabilizing to destabilizing depending on calibration.

In short, leadership is not an exogenous backdrop but a mathematical modulator of systemic resilience. This perspective reframes classical debates in comparative politics by embedding agency into dynamical systems theory: leaders are not only decision-makers but also parametric forces that reshape the geometry of stability landscapes.

B. Policy implications: stability trade-offs between repression and legitimacy

Our formalism highlights that repression and legitimacy constitute two distinct yet interacting levers of systemic stability. Both can raise the critical stress threshold , but they operate through fundamentally different mechanisms and with opposite long-term consequences.

Repression provides a short-term stabilizer. By suppressing protest intensity (), it can temporarily delay bifurcation, giving the impression of resilience. Yet our U-shaped sensitivity curve shows that beyond a moderate level, repression generates backlash, lowering consensus and accelerating trust erosion (). Excessive reliance on coercion therefore creates a paradox of stability: the more repression is applied, the shallower the equilibrium becomes, making the system brittle under even modest shocks.

Legitimacy, by contrast, is a long-term stabilizer. It directly strengthens the trust feedback loop, raising and allowing the system to absorb greater stress before tipping. Unlike repression, legitimacy does not exhibit strong diminishing returns within the observed range. Higher legitimacy consistently pushes upward, broadening the stability basin and suppressing the black horse potential ().

The trade-off is thus clear: repression offers speed but erodes durability, while legitimacy offers durability but requires sustained investment in democratic credibility and inclusiveness. For policymakers, this means that the pursuit of stability through coercion is ultimately self-defeating, while strategies that cultivate legitimacy---transparent governance, electoral credibility, equitable policies---extend resilience in a structurally robust way.

In practical terms, our model suggests that regimes prioritizing legitimacy over repression are better insulated against systemic shocks. Conversely, those that substitute repression for legitimacy may appear stable but are in fact perched near bifurcation thresholds, where instability can be triggered abruptly by relatively small exogenous disturbances.

C. Integration with bifurcation theory: from abstract dynamics to political diagnostics

The incorporation of leadership parameters into a bifurcation framework allows us to move beyond abstract dynamical analysis toward concrete political diagnostics. Traditional bifurcation theory identifies when equilibria lose stability as a control parameter crosses a critical value. In our formalism, exogenous stress () functions as the primary bifurcation driver, while leadership parameters determine the location of the critical threshold and the qualitative nature of the transition.

This integration produces two key advances. First, it makes bifurcation analysis politically interpretable: a saddle--node bifurcation corresponds to abrupt regime collapse (as in Soeharto's case), while a Hopf bifurcation corresponds to recurrent protest cycles and oscillatory instability (a possibility in fragile democratic contexts). Second, it allows the calibration of early-warning indicators. By monitoring trust levels (), protest intensity (), and black horse potential (), policymakers and analysts can detect when a system is approaching a bifurcation boundary.

In this sense, our contribution extends bifurcation theory from abstract nonlinear systems into the domain of political diagnostics. Just as in physics or biology, where bifurcations indicate tipping points in ecosystems or chemical reactions, in politics they represent thresholds of stability where small shocks can precipitate disproportionate change. Leadership qualities, modeled as parametric modulators, thus become not only normative variables but diagnostic instruments that determine how close a polity is to its tipping point.

D. Limitations and directions for future work

While our formalism demonstrates that leadership parameters can be mathematically integrated into bifurcation models of political stability, several limitations must be acknowledged.

First, the model relies on stylized functional forms and simplified mappings of complex political realities. Leadership qualities such as legitimacy or narrative capacity were normalized into single coefficients, which inevitably compresses multidimensional phenomena into scalar terms. This abstraction is necessary for tractability, but it obscures contextual nuances such as cultural variation, institutional legacies, and the role of media ecosystems.

Second, the simulations used heuristic parameter values calibrated from case knowledge rather than empirically estimated coefficients. Although the qualitative dynamics align with historical trajectories of Soeharto, SBY, Jokowi, and Prabowo, a more rigorous approach would involve empirical fitting to longitudinal political and economic data. Integrating survey-based trust indices, protest event datasets, and macroeconomic indicators could strengthen the model's predictive validity.

Third, the current framework treats leadership as exogenous modulators. In reality, leadership strategies evolve endogenously in response to crisis and opposition. Future work should extend the formalism into adaptive or game-theoretic settings, where leaders adjust repression, narrative, or coalition strategies dynamically, and citizens adapt their responses accordingly.

Finally, the model has not yet incorporated international and transnational variables, such as foreign capital flows, geopolitical shocks, or social media contagion. These exogenous forces may interact with domestic leadership parameters in nonlinear ways, potentially shifting bifurcation thresholds in unexpected directions.

Future research should therefore aim at three directions: (i) empirical calibration of the leadership--stability mapping; (ii) integration of adaptive leadership strategies into the dynamical system; and (iii) expansion to multi-level systems that capture both domestic and global shocks. Such developments would enhance both the explanatory depth and the predictive utility of the framework, advancing the project of a mathematically grounded science of political stability.

VII. Conclusion

A. Leadership quality is quantifiable as parameters in nonlinear dynamics

This study demonstrates that leadership quality, often treated as a qualitative or normative construct, can be operationalized within the rigorous language of nonlinear dynamics. By defining seven leadership parameters---consensus, legitimacy, crisis management, narrative control, economic stability, elite management, and repression balance---we embed leadership into the coefficients of a dynamical system governing trust, economic stress, protest intensity, and black horse potential.

The result is a tractable formalism in which leadership directly modulates bifurcation thresholds and stability landscapes. Our simulations and case studies show that high legitimacy consistently raises the critical stress threshold, narrative control provides resilience with diminishing returns, and repression offers only transient stability before backfire effects emerge. Leadership therefore operates as a mathematical modulator: its qualities determine not only how long stability can persist but also the manner in which instability unfolds---abrupt collapse, oscillatory cycles, or gradual erosion.

In this way, leadership becomes quantifiable: it is no longer an abstract cultural variable but a set of parameters that shift systemic dynamics in predictable ways. This formalization offers a foundation for building a comparative political science of tipping points, where leadership is understood as both agency and structure, both narrative and number.

B. Integration of leadership parameters with bifurcation theory provides novel insights into political unrest and regime transitions

By embedding leadership parameters within a bifurcation-theoretic framework, we gain a novel vantage point on the mechanisms that govern political unrest and regime transitions. Traditional accounts of instability emphasize exogenous shocks---economic crises, geopolitical disruptions, or demographic pressures. Our model shows that while such shocks are critical, their systemic consequences are filtered through the internal architecture of leadership.

This integration yields three insights. First, regime stability is not determined solely by the magnitude of external stress but by how leadership parameters shift the critical thresholds at which bifurcations occur. Second, the qualitative type of transition---whether abrupt collapse, cyclical protest waves, or gradual erosion---depends on the nonlinear configuration of legitimacy, narrative, repression, and coalition management. Third, outsider emergence is not an anomaly but a mathematically predictable consequence of instability: the black horse potential () rises systematically when endogenous leadership parameters fail to stabilize trust and suppress protest growth.

The synthesis of leadership analysis with bifurcation theory therefore transforms political instability from a largely descriptive domain into one with predictive and diagnostic potential. It allows political scientists to move beyond post hoc explanations of regime change toward anticipatory modeling, in which early-warning indicators are tied explicitly to parametric weaknesses in leadership. In doing so, it opens the door to a comparative science of regime transitions that is both conceptually rich and mathematically precise.

C. Toward a predictive science of political stability

The integration of leadership parameters into nonlinear dynamical models opens a pathway toward a predictive science of political stability. By quantifying leadership attributes as parameters that modulate bifurcation thresholds, it becomes possible to construct models that generate early-warning indicators of unrest, identify conditions under which outsider figures are likely to emerge, and anticipate the qualitative form of regime transitions.

This predictive capacity rests on three pillars. First, the mathematical formalism ensures internal coherence: the same system of equations can be applied across cases, enabling comparative analysis. Second, the sensitivity of critical thresholds to leadership parameters provides a structured means to diagnose vulnerability---whether a regime is brittle like Soeharto's, resilient like SBY's, or stable-until-sudden-shift like Jokowi's. Third, the model highlights the conditions under which legitimacy deficits and repression surpluses converge to lower the tipping threshold, as illustrated by the Prabowo case.

Moving toward predictive application will require empirical calibration, integration of dynamic adaptation, and incorporation of global interdependencies. Yet even in its stylized form, the model demonstrates that political stability is neither mysterious nor wholly contingent: it follows discernible mathematical patterns shaped by leadership. By uniting bifurcation theory with leadership analysis, we take a step toward transforming political stability studies into a predictive, diagnostic, and comparative science.

References

Core Books & Theoretical Foundations

Acemoglu, D., & Robinson, J. A. (2012). Why Nations Fail: The Origins of Power, Prosperity, and Poverty. New York: Crown Business.

Goldstone, J. A. (1991). Revolution and Rebellion in the Early Modern World. Berkeley: University of California Press.

Hirschman, A. O. (1970). Exit, Voice, and Loyalty: Responses to Decline in Firms, Organizations, and States. Cambridge, MA: Harvard University Press.

Kuznetsov, Y. A. (2004). Elements of Applied Bifurcation Theory (3rd ed.). New York: Springer.

Prigogine, I., & Stengers, I. (1984). Order Out of Chaos: Man's New Dialogue with Nature. New York: Bantam Books.

Scheffer, M. (2009). Critical Transitions in Nature and Society. Princeton: Princeton University Press.

Strogatz, S. H. (2018). Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (2nd ed.). Boca Raton, FL: CRC Press.

Politics, Instability, and Indonesia Case Literature

Aspinall, E., & Mietzner, M. (Eds.). (2010). Problems of Democratisation in Indonesia: Elections, Institutions and Society. Singapore: ISEAS Publishing.

Mietzner, M. (2024). Elite Collusion in Indonesia: How It Has Both Enabled and Limited Executive Aggrandizement. The Annals of the American Academy of Political and Social Science, 75(1).

Robison, R., & Hadiz, V. R. (2004). Reorganising Power in Indonesia: The Politics of Oligarchy in an Age of Markets. London: RoutledgeCurzon.

Recent Journal Articles (2023--2025)

Doces, J. A., & Lupton, D. L. (2023). Leader Turnover, Regime Type, and Growth Volatility. Global Studies Quarterly, 2.

Kuzenbayev, N., & Pelizzo, R. (2023). Political Stability, Confidence in the Future, and Values. Social Sciences, 12(2), Article 82.

Setiawan, K. M. P., & Tomsa, D. (2023). Defending a Vulnerable yet Resilient Democracy: Civil Society Activism in Jokowi's Indonesia. Asian Journal of Comparative Politics.

Putri, F. A. (2024). Political Opportunity, Media Framing and Resource Mobilization: Social Protest Against Indonesia's Regional Election Law Revisions. Journal of Kuasa, 1(2), 64--71.

Savirani, A. (2023). University Curricula and Student Movement in Indonesia: "Gejayan Memanggil" Protest in Yogyakarta. Journal of Governance, 8(3).

Misran, M., Sutan, A. J., & Efendi, D. (2023). Social Media Protest and Resistance: Case of Indonesia Job Creation Bill Policy. Jurnal Komunikasi, 15(1), 1--17.

Regional Stability & Societal Systems (Comparative Perspectives)

Salma, U., & Khan, M. F. H. (2023). The Connection Between Political Stability and Inflation: Insights from Four South Asian Nations. arXiv preprint.

Tilman, A. R., Krueger, E. H., McManus, L. C., & Watson, J. R. (2023). Maintaining human wellbeing as socio-environmental systems undergo regime shifts. arXiv preprint.