When the Jacobian $Jf_p$ of a smooth system $\dot{X} = F(X)$ at an equilibrium $p$ has a pair of purely imaginary eigenvalues $\pm i\omega_0$ with $\omega_0 > 0$ — and all other eigenvalues have nonzero real part — the linearisation produces exact rotation at frequency $\omega_0$. It neither contracts nor expands, making the linear theory inconclusive about stability. Resolving this ambiguity requires examining the nonlinear terms through **complexification**, which isolates the single scalar quantity — the First Lyapunov Coefficient — that determines stability. When a parameter is introduced and the eigenvalues cross the imaginary axis, the same coefficient governs the **Hopf bifurcation**.

## Definition

[definition: Equilibrium With Purely Imaginary Eigenvalues] Let $F: U \subseteq \mathbb{R}^n \to \mathbb{R}^n$ be a smooth vector field with an equilibrium $p$ (i.e., $F(p) = 0$). The equilibrium has **purely imaginary eigenvalues** if the Jacobian $A := Jf_p$ has a simple pair of eigenvalues $\lambda_{1,2} = \pm i\omega_0$ with $\omega_0 > 0$, and every other eigenvalue $\lambda_j$ satisfies $\operatorname{Re}(\lambda_j) \neq 0$. [/definition]

The condition $\operatorname{Re}(\lambda_{1,2}) = 0$ places the equilibrium exactly on the stability [boundary](/page/Boundary): the imaginary axis for [continuous](/page/Continuity) systems plays the role of the unit circle for maps (see [Fixed Points of Maps](/page/Fixed%20Points%20of%20Maps)). The linear flow on the center eigenspace is rigid rotation — all orbits are circles of constant radius. Whether the full nonlinear orbits spiral inward (stable) or outward (unstable) depends entirely on higher-order terms.

When $n = 2$, the entire phase plane is the center eigenspace and no reduction is needed. When $n > 2$, the [Center Manifold Theorem](/page/Center%20Manifold%20Theorem) provides a 2D invariant surface on which the stability is decided, and the remaining directions are exponentially contracting or expanding.

## Complexification

The nonlinear analysis begins by writing the system in complex coordinates that diagonalise the linear part. In real coordinates, the linearisation is a $2 \times 2$ rotation matrix that mixes the two variables at every step. In complex coordinates, rotation becomes multiplication by $e^{i\omega_0 t}$, cleanly separating the rotational dynamics from the amplitude dynamics.

[definition: Complexification Of A Planar System] Consider a smooth system on $\mathbb{R}^2$ with an equilibrium at the origin whose linearisation has eigenvalues $\pm i\omega_0$. After a real linear change of coordinates (if needed) to put the linear part in rotation form $\begin{pmatrix} 0 & -\omega_0 \\ \omega_0 & 0 \end{pmatrix}$, define the **complex coordinate**: \begin{align*} z := x_1 - ix_2, \qquad \bar{z} = x_1 + ix_2, \end{align*} with inverse $x_1 = \frac{z + \bar{z}}{2}$, $x_2 = \frac{z - \bar{z}}{-2i}$. In these coordinates the system becomes: \begin{align*} \dot{z} = i\omega_0 z + \sum_{j+k \ge 2} \frac{1}{j!\,k!}\,g_{jk}\,z^j\bar{z}^k, \end{align*} where the Taylor coefficients $g_{jk} \in \mathbb{C}$ are determined by the nonlinear terms of the original system. [/definition]

The computation proceeds mechanically: [set](/page/Set) $z = x_1 - ix_2$, compute $\dot{z} = \dot{x}_1 - i\dot{x}_2$ by substituting the original ODEs, then re-express everything in terms of $z$ and $\bar{z}$ using the inverse relations. The result is matched to the expansion $\dot{z} = i\omega_0 z + \frac{1}{2}g_{20}z^2 + g_{11}z\bar{z} + \frac{1}{2}g_{02}\bar{z}^2 + \frac{1}{2}g_{21}z^2\bar{z} + \cdots$ to read off the coefficients $g_{20}, g_{11}, g_{02}, g_{21}$.

[example: Complexification Of A Planar System] Consider the system: \begin{align*} \dot{x} &= y, \\ \dot{y} &= -x - 6x^2 + 2xy. \end{align*} **Step 1: Eigenvalues.** The Jacobian at the origin is $A = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$ with eigenvalues $\pm i$, so $\omega_0 = 1$. The linear part is already in rotation form (with the convention that gives $\dot{z} = iz$ for $z = x - iy$). **Step 2: Complex coordinate.** Set $z = x - iy$, so $x = \frac{z+\bar{z}}{2}$ and $y = i\frac{z - \bar{z}}{2}$. Compute $\dot{z} = \dot{x} - i\dot{y}$: \begin{align*} \dot{z} = y + ix + 6ix^2 - 2ixy. \end{align*} The linear part gives $y + ix = i(x - iy) = iz$. For the nonlinear terms, substitute the inverse relations: \begin{align*} 6ix^2 &= \frac{6i}{4}(z + \bar{z})^2 = \frac{3i}{2}z^2 + 3iz\bar{z} + \frac{3i}{2}\bar{z}^2, \\ -2ixy &= -2i \cdot \frac{z+\bar{z}}{2} \cdot i\frac{z - \bar{z}}{2} = \frac{1}{2}z^2 - \frac{1}{2}\bar{z}^2. \end{align*} Combining all terms: \begin{align*} \dot{z} = iz + \left(\tfrac{1}{2} + \tfrac{3i}{2}\right)z^2 + 3iz\bar{z} + \left(-\tfrac{1}{2} + \tfrac{3i}{2}\right)\bar{z}^2. \end{align*} **Step 3: Read off coefficients.** Matching to $\dot{z} = iz + \frac{1}{2}g_{20}z^2 + g_{11}z\bar{z} + \frac{1}{2}g_{02}\bar{z}^2 + \cdots$: \begin{align*} g_{20} &= 1 + 3i, \qquad g_{11} = 3i, \qquad g_{02} = -1 + 3i, \qquad g_{21} = 0. \end{align*} There is no $z^2\bar{z}$ term in the expansion, so $g_{21} = 0$. [/example]

## The First Lyapunov Coefficient

With the system in complex form, the stability reduces to a single number: the coefficient of the leading nonlinear term that systematically changes the amplitude of oscillation. The following motivation explains why this is the $z^2\bar{z}$ term and how the formula arises.

[motivation] ### Polar Coordinates and the Amplitude Equation Substituting $z = re^{i\theta}$ into $\dot{z} = i\omega_0 z + \sum c_{jk}z^j\bar{z}^k$ and using $z^j\bar{z}^k = r^{j+k}e^{i(j-k)\theta}$, the amplitude equation becomes: \begin{align*} \dot{r} = \sum_{j+k \ge 2} r^{j+k}\,\operatorname{Re}\!\left(c_{jk}\,e^{i(j-k-1)\theta}\right). \end{align*} Each monomial contributes a term to $\dot{r}$ that oscillates in $\theta$ at frequency $j - k - 1$. To determine whether $r$ systematically grows or shrinks, we compute the net change in $r$ over one full rotation by integrating over $\theta$ from $0$ to $2\pi$: \begin{align*} \Delta r = \frac{1}{\omega_0}\int_0^{2\pi} \dot{r}\,d\theta + O(r^5). \end{align*} The change of integration variable from $t$ to $\theta$ uses $d\theta = \dot{\theta}\,dt \approx \omega_0\,dt$, with corrections of order $r^2$ that only affect higher-order terms. ### Resonant vs Non-Resonant Terms Each [integral](/page/Integral) evaluates to: \begin{align*} \int_0^{2\pi} e^{in\theta}\,d\theta = \begin{cases} 2\pi & \text{if } n = 0, \\ 0 & \text{if } n \neq 0. \end{cases} \end{align*} A monomial $z^j\bar{z}^k$ produces $n = j - k - 1$ in the amplitude equation. When $n \neq 0$, the contribution oscillates and integrates to exactly zero — the term pushes $r$ up for part of the rotation and pulls it down equally for the rest, with no net effect. When $n = 0$ (i.e., $j = k + 1$), the contribution is constant in $\theta$ and produces a nonzero cumulative change. These are the **resonant terms**: the only monomials that systematically affect the amplitude. ### Why Cubic Is the Leading Order At quadratic order ($j + k = 2$), the resonance condition $j = k + 1$ gives $j = 3/2$, which is not an integer — no quadratic monomial is resonant. At cubic order ($j + k = 3$), it gives $j = 2, k = 1$: the monomial $z^2\bar{z}$. This is the leading resonant term, and its coefficient determines: \begin{align*} \Delta r = \frac{2\pi}{\omega_0}\operatorname{Re}(c_{21})\,r_0^3 + O(r_0^5). \end{align*} ### The Indirect Contribution of Quadratic Terms Although quadratic terms have zero direct contribution to $\Delta r$, they distort the trajectory within each rotation. A near-identity coordinate change $z \mapsto z + h_{20}z^2 + h_{11}z\bar{z} + h_{02}\bar{z}^2$ eliminates them, but generates new cubic terms through interaction with the linear part. The $L_1$ formula accounts for both the direct resonant cubic term ($g_{21}$) and the indirect cubic contributions generated by eliminating the quadratics ($ig_{20}g_{11}$). [/motivation]

[definition: Resonant Monomial] A monomial $z^j\bar{z}^k$ in the complex expansion of a planar system with linear part $\dot{z} = i\omega_0 z$ is **resonant** if $j = k + 1$. Equivalently, its contribution to the amplitude equation $\dot{r}$ is independent of $\theta$ and does not integrate to zero over one rotation. The resonant monomials at successive orders are $z^2\bar{z}$ (cubic), $z^3\bar{z}^2$ (quintic), $z^4\bar{z}^3$ (septic), with corresponding Lyapunov coefficients $L_1, L_2, L_3$. [/definition]

After performing successive near-identity coordinate changes to eliminate all non-resonant monomials, the system reduces to its **normal form**, in which only resonant terms survive:

\begin{align*} \dot{z} = i\omega_0 z + c_1 z^2\bar{z} + c_2 z^3\bar{z}^2 + c_3 z^4\bar{z}^3 + \cdots \end{align*}

where $c_k \in \mathbb{C}$. In polar coordinates $z = re^{i\theta}$, each resonant monomial $z^{k+1}\bar{z}^k = r^{2k+1}e^{i\theta}$ contributes purely to the amplitude, giving:

\begin{align*} \dot{r} = \operatorname{Re}(c_1)\,r^3 + \operatorname{Re}(c_2)\,r^5 + \operatorname{Re}(c_3)\,r^7 + \cdots \end{align*}

[definition: Lyapunov Coefficients] The **$k$-th Lyapunov coefficient** is \begin{align*} L_k := \operatorname{Re}(c_k), \end{align*} where $c_k$ is the coefficient of the $k$-th resonant monomial $z^{k+1}\bar{z}^k$ in the normal form. The amplitude equation becomes: \begin{align*} \dot{r} = L_1 r^3 + L_2 r^5 + L_3 r^7 + \cdots \end{align*} The stability of the equilibrium is determined by the **first nonzero** $L_k$: for $r$ sufficiently small, that term dominates all higher powers of $r$, so $\operatorname{sgn}(\dot{r}) = \operatorname{sgn}(L_k)$. [/definition]

The sign of the leading nonzero Lyapunov coefficient determines stability. The following result makes this precise and provides the formula for $L_1$.