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:
