[proofplan] We reduce to a canonical form by splitting coordinates along the invariant subspaces $E^s$, $E^u$, and $E^c$. The centre manifold is then sought as the graph of a $C^r$ function $h: E^c \to E^s \oplus E^u$ satisfying $h(\mathbf{0}) = \mathbf{0}$ and $Dh_{\mathbf{0}} = 0$ (the tangency condition). Substituting the graph ansatz into the ODE yields a quasilinear PDE for $h$ — the **invariance equation**. We solve this equation by reformulating it as a fixed-point problem in a Banach space of $C^r$ functions with small Lipschitz constant, and invoke the Banach fixed-point theorem to obtain a unique solution in that space. The stable and unstable manifolds are constructed identically (and follow already from the [Stable Manifold Theorem](/theorems/2778) applied to the hyperbolic subsystem). Non-uniqueness of the centre manifold is demonstrated by an explicit counterexample. [/proofplan] [step:Decompose the system along the invariant subspaces $E^s$, $E^u$, $E^c$] The Jacobian matrix $A := Jf_{\mathbf{0}} \in \mathbb{R}^{n \times n}$ has eigenvalues that partition into three groups: - **Stable eigenvalues:** $\operatorname{Re}(\lambda) < 0$, spanning the stable subspace $E^s$ of dimension $n_s$. - **Unstable eigenvalues:** $\operatorname{Re}(\lambda) > 0$, spanning the unstable subspace $E^u$ of dimension $n_u$. - **Centre eigenvalues:** $\operatorname{Re}(\lambda) = 0$, spanning the centre subspace $E^c$ of dimension $n_c$. Since $\mathbf{0}$ is non-hyperbolic, $n_c \geq 1$. The space $\mathbb{R}^n$ decomposes as a direct sum $\mathbb{R}^n = E^s \oplus E^u \oplus E^c$, and each of these subspaces is $A$-invariant. Choose a basis adapted to this decomposition, so that writing $x = (x_s, x_u, x_c)$ with $x_s \in E^s \cong \mathbb{R}^{n_s}$, $x_u \in E^u \cong \mathbb{R}^{n_u}$, $x_c \in E^c \cong \mathbb{R}^{n_c}$, the system $\dot{x} = f(x)$ takes the form \begin{align*} \dot{x}_s &= A_s x_s + g_s(x_s, x_u, x_c), \\ \dot{x}_u &= A_u x_u + g_u(x_s, x_u, x_c), \\ \dot{x}_c &= A_c x_c + g_c(x_s, x_u, x_c), \end{align*} where $A_s \in \mathbb{R}^{n_s \times n_s}$ has all eigenvalues with $\operatorname{Re}(\lambda) < 0$, $A_u \in \mathbb{R}^{n_u \times n_u}$ has all eigenvalues with $\operatorname{Re}(\lambda) > 0$, $A_c \in \mathbb{R}^{n_c \times n_c}$ has all eigenvalues with $\operatorname{Re}(\lambda) = 0$, and the nonlinear terms $g_s, g_u, g_c$ are $C^r$ functions satisfying \begin{align*} g_\sigma(\mathbf{0}) = \mathbf{0}, \qquad D(g_\sigma)_{\mathbf{0}} = 0, \qquad \sigma \in \{s, u, c\}. \end{align*} The vanishing of $g_\sigma$ and its derivative at the origin follows from $f(\mathbf{0}) = \mathbf{0}$ and the definition of $A = Jf_{\mathbf{0}}$. [guided] The first step is purely algebraic: we decompose $\mathbb{R}^n$ according to the spectral properties of the Jacobian $A = Jf_{\mathbf{0}} \in \mathbb{R}^{n \times n}$. The eigenvalues of $A$ partition into three groups based on their real parts: - **Stable eigenvalues** with $\operatorname{Re}(\lambda) < 0$, whose generalised eigenspaces span $E^s$ of dimension $n_s$. - **Unstable eigenvalues** with $\operatorname{Re}(\lambda) > 0$, spanning $E^u$ of dimension $n_u$. - **Centre eigenvalues** with $\operatorname{Re}(\lambda) = 0$, spanning $E^c$ of dimension $n_c$. Since $\mathbf{0}$ is non-hyperbolic, $n_c \geq 1$. The three subspaces are each $A$-invariant (generalised eigenspaces are invariant under the matrix that defines them) and together span $\mathbb{R}^n$ because these three cases exhaust all possibilities. Hence $\mathbb{R}^n = E^s \oplus E^u \oplus E^c$. Choosing a basis adapted to this decomposition puts $A$ in block-diagonal form $A = \operatorname{diag}(A_s, A_u, A_c)$. Writing $x = (x_s, x_u, x_c)$, the system $\dot{x} = f(x)$ splits as \begin{align*} \dot{x}_s &= A_s x_s + g_s(x_s, x_u, x_c), \\ \dot{x}_u &= A_u x_u + g_u(x_s, x_u, x_c), \\ \dot{x}_c &= A_c x_c + g_c(x_s, x_u, x_c), \end{align*} where the nonlinear terms $g_\sigma$ satisfy $g_\sigma(\mathbf{0}) = \mathbf{0}$ and $D(g_\sigma)_{\mathbf{0}} = 0$. These conditions say that $g_\sigma(x) = O(|x|^2)$ near the origin -- the nonlinearities vanish to at least second order. This is the standard consequence of having already extracted the linear part $Ax$ from $f(x) = Ax + (f(x) - Ax)$, using $f(\mathbf{0}) = \mathbf{0}$ and $Jf_{\mathbf{0}} = A$. Why does this decomposition matter for the rest of the proof? The three linear blocks have qualitatively different behaviours: $e^{A_s t}$ decays exponentially as $t \to +\infty$, $e^{A_u t}$ grows exponentially as $t \to +\infty$ (but decays as $t \to -\infty$), and $e^{A_c t}$ neither grows nor decays -- it oscillates or remains polynomially bounded. The centre manifold will be a surface tangent to $E^c$ on which the dynamics is governed by the centre block $A_c$ and the nonlinear terms restricted to that surface. The different exponential rates of $A_s$ and $A_u$ are what allow the fixed-point argument (in later steps) to succeed: they provide the contraction needed for the integral operators to converge. [/guided] [/step] [step:Formulate the centre manifold as the graph of a function $h$ and derive the invariance equation] We seek the centre manifold $W^c$ locally as the graph of a $C^r$ function \begin{align*} h: B_\rho^c \subset E^c &\to E^s \oplus E^u, \\ x_c &\mapsto h(x_c) = (h_s(x_c),\, h_u(x_c)), \end{align*} where $B_\rho^c = \{x_c \in E^c : |x_c| < \rho\}$ for some $\rho > 0$, subject to the tangency conditions \begin{align*} h(\mathbf{0}) = \mathbf{0}, \qquad Dh_{\mathbf{0}} = 0. \end{align*} The candidate centre manifold is then $W^c_{\mathrm{loc}} = \{(x_s, x_u, x_c) : (x_s, x_u) = h(x_c),\; x_c \in B_\rho^c\}$. For $W^c_{\mathrm{loc}}$ to be invariant under the flow, the trajectory starting at any point $(x_s, x_u, x_c) = (h_s(x_c), h_u(x_c), x_c)$ on the graph must remain on the graph. Differentiating the relation $x_s(t) = h_s(x_c(t))$ with respect to $t$ and using the chain rule yields \begin{align*} \dot{x}_s = Dh_{s,x_c}\, \dot{x}_c. \end{align*} Substituting the ODE for $\dot{x}_s$ and $\dot{x}_c$ from the decomposed system, and evaluating the nonlinear terms on the graph $(x_s, x_u) = h(x_c)$, we obtain the **invariance equation** for $h_s$: \begin{align*} A_s\, h_s(x_c) + g_s(h(x_c), x_c) = Dh_{s,x_c}\!\left[A_c\, x_c + g_c(h(x_c), x_c)\right]. \end{align*} An analogous equation holds for $h_u$. Together, these form a system of quasilinear PDEs for $h = (h_s, h_u)$. [guided] We seek the centre manifold $W^c$ locally as the graph of a $C^r$ function $h: B_\rho^c \subset E^c \to E^s \oplus E^u$ satisfying $h(\mathbf{0}) = \mathbf{0}$ and $Dh_{\mathbf{0}} = 0$. The graph ansatz $(x_s, x_u) = h(x_c)$ is the standard way to represent a manifold passing through the origin and tangent to $E^c$: the condition $h(\mathbf{0}) = \mathbf{0}$ places the graph at the origin, and $Dh_{\mathbf{0}} = 0$ ensures that the tangent plane at $\mathbf{0}$ is exactly $E^c$ (since the graph of a function with vanishing derivative at a point is tangent to the domain at that point). The candidate centre manifold is $W^c_{\mathrm{loc}} = \{(x_s, x_u, x_c) : (x_s, x_u) = h(x_c),\; x_c \in B_\rho^c\}$. How do we derive the invariance equation? For $W^c_{\mathrm{loc}}$ to be invariant under the flow, a trajectory starting on the graph must remain on the graph. This means $x_s(t) = h_s(x_c(t))$ must hold for all $t$. Differentiating both sides with respect to $t$ gives, by the chain rule: \begin{align*} \dot{x}_s = Dh_{s,x_c}\, \dot{x}_c. \end{align*} The left-hand side is the ODE for $x_s$: $\dot{x}_s = A_s x_s + g_s(x_s, x_u, x_c)$. The right-hand side uses the ODE for $x_c$: $\dot{x}_c = A_c x_c + g_c(x_s, x_u, x_c)$. Evaluating everything on the graph (replacing $(x_s, x_u)$ by $h(x_c)$) and equating gives the **invariance equation** for $h_s$: \begin{align*} A_s\, h_s(x_c) + g_s(h(x_c), x_c) = Dh_{s,x_c}\!\left[A_c\, x_c + g_c(h(x_c), x_c)\right]. \end{align*} An analogous equation holds for $h_u$. This is a **quasilinear PDE** for $h$: the highest derivative $Dh_{s,x_c}$ appears multiplied by terms that depend on $h$ itself (through $g_c(h(x_c), x_c)$). Solving it directly by PDE methods is difficult, which motivates the reformulation as a fixed-point problem in a Banach space in the subsequent steps. [/guided] [/step] [step:Cut off the nonlinearity to obtain a globally Lipschitz system] The nonlinear terms $g_\sigma$ are only defined locally and may have large Lipschitz constants away from the origin. To work in a Banach space setting and apply the contraction mapping principle, we introduce a smooth cut-off. Let $\chi: \mathbb{R} \to [0,1]$ be a $C^\infty$ bump function with $\chi(s) = 1$ for $|s| \leq 1$ and $\chi(s) = 0$ for $|s| \geq 2$. For $\rho > 0$, define the modified nonlinearities \begin{align*} \tilde{g}_\sigma(x_s, x_u, x_c) := \chi\!\left(\frac{|(x_s, x_u, x_c)|}{\rho}\right) g_\sigma(x_s, x_u, x_c), \qquad \sigma \in \{s, u, c\}. \end{align*} The modified system $\dot{x} = Ax + \tilde{g}(x)$ agrees with the original system inside the ball $B(\mathbf{0}, \rho)$ and is identically linear outside $B(\mathbf{0}, 2\rho)$. Since $D(g_\sigma)_{\mathbf{0}} = 0$, the Lipschitz constant of $\tilde{g}_\sigma$ can be made arbitrarily small by taking $\rho$ sufficiently small: \begin{align*} \operatorname{Lip}(\tilde{g}_\sigma) \leq C\rho \end{align*} for a constant $C > 0$ depending on the $C^2$ norm of $g_\sigma$ near the origin. We fix $\rho$ small enough that $\operatorname{Lip}(\tilde{g}_\sigma) < \varepsilon$ for a threshold $\varepsilon > 0$ to be determined below. [guided] Why do we need a cut-off? The contraction mapping argument in the next step requires working in a Banach space of functions defined on all of $E^c$, with operators that are global contractions. This requires global control on the nonlinearity -- specifically, a small global Lipschitz constant. But the original $g_\sigma$ is only defined on a neighbourhood $E$ of the origin and may have large Lipschitz constants away from $\mathbf{0}$. The solution is to multiply $g_\sigma$ by a smooth bump function $\chi$ that equals $1$ near $\mathbf{0}$ and $0$ far from $\mathbf{0}$. Let $\chi: \mathbb{R} \to [0,1]$ be $C^\infty$ with $\chi(s) = 1$ for $|s| \leq 1$ and $\chi(s) = 0$ for $|s| \geq 2$. The modified nonlinearity \begin{align*} \tilde{g}_\sigma(x) := \chi\!\left(\frac{|x|}{\rho}\right) g_\sigma(x) \end{align*} agrees with $g_\sigma$ inside $B(\mathbf{0}, \rho)$ (where $\chi = 1$) and vanishes identically outside $B(\mathbf{0}, 2\rho)$ (where $\chi = 0$). The key estimate is $\operatorname{Lip}(\tilde{g}_\sigma) \leq C\rho$, which can be made arbitrarily small by choosing $\rho$ small. Why does this hold? Since $D(g_\sigma)_{\mathbf{0}} = 0$ and $g_\sigma$ is $C^r$ ($r \geq 2$), the mean value inequality gives $|D(g_\sigma)_{x}| \leq C_1 |x|$ for $|x| \leq 2\rho$. When we differentiate the product $\tilde{g}_\sigma = \chi \cdot g_\sigma$, the product rule gives $D\tilde{g}_\sigma = D\chi \cdot g_\sigma + \chi \cdot Dg_\sigma$. The first term is bounded by $|D\chi| \cdot |g_\sigma| = O(1/\rho) \cdot O(|x|^2) = O(\rho)$ (since $|x| \leq 2\rho$ on the support). The second term is bounded by $|\chi| \cdot |Dg_\sigma| = 1 \cdot O(|x|) = O(\rho)$. Combining: $|D\tilde{g}_\sigma| = O(\rho)$ everywhere, so $\operatorname{Lip}(\tilde{g}_\sigma) \leq C\rho$. The modified system $\dot{x} = Ax + \tilde{g}(x)$ is identical to the original inside $B(\mathbf{0}, \rho)$ and is purely linear ($\dot{x} = Ax$) outside $B(\mathbf{0}, 2\rho)$. Any invariant manifold we construct for the modified system restricts to a local invariant manifold of the original system within $B(\mathbf{0}, \rho)$. We fix $\rho$ small enough that $\operatorname{Lip}(\tilde{g}_\sigma) < \varepsilon$ for a threshold $\varepsilon > 0$ determined by the contraction condition in the next step. [/guided] [/step] [step:Reformulate as a fixed-point problem in a Banach space and apply the contraction mapping principle] We reformulate the invariance equation as a fixed-point problem. For the stable component $h_s$, integrate the $x_s$-equation forward in time along trajectories on the manifold. On the centre manifold, the dynamics is governed by \begin{align*} \dot{x}_c = A_c\, x_c + \tilde{g}_c(h(x_c), x_c). \end{align*} Denote by $x_c(t; y)$ the solution of this reduced equation with $x_c(0; y) = y$. The stable component must satisfy \begin{align*} h_s(y) = \int_{-\infty}^{0} e^{-A_s \tau}\, \tilde{g}_s\bigl(h(x_c(\tau; y)),\, x_c(\tau; y)\bigr)\, d\mathcal{L}^1(\tau). \end{align*} This integral converges because $e^{-A_s \tau}$ decays exponentially as $\tau \to -\infty$ (since all eigenvalues of $A_s$ have strictly negative real part, $e^{A_s t}$ decays as $t \to +\infty$, so $e^{-A_s \tau} = e^{A_s(-\tau)}$ decays as $\tau \to -\infty$) and the integrand is bounded by $\operatorname{Lip}(\tilde{g}_s) \cdot \sup |x_c(\tau; y)|$. Similarly, for the unstable component, integrate the $x_u$-equation backward: \begin{align*} h_u(y) = -\int_{0}^{\infty} e^{-A_u \tau}\, \tilde{g}_u\bigl(h(x_c(\tau; y)),\, x_c(\tau; y)\bigr)\, d\mathcal{L}^1(\tau). \end{align*} This integral converges because $e^{-A_u \tau}$ decays exponentially as $\tau \to +\infty$ (all eigenvalues of $A_u$ have strictly positive real part). Define the Banach space $\mathcal{X} = \{h \in C^0(E^c;\, E^s \oplus E^u) : h(\mathbf{0}) = \mathbf{0},\; \|h\|_{\mathrm{Lip}} < \infty\}$ equipped with the Lipschitz norm $\|h\|_{\mathrm{Lip}} = \sup_{y \neq y'} |h(y) - h(y')| / |y - y'|$, and let $\mathcal{X}_\delta = \{h \in \mathcal{X} : \|h\|_{\mathrm{Lip}} \leq \delta\}$ for a small $\delta > 0$. The integral equations above define an operator \begin{align*} \mathcal{T}: \mathcal{X}_\delta &\to \mathcal{X}, \\ h &\mapsto \mathcal{T}[h], \end{align*} where $\mathcal{T}[h] = (\mathcal{T}_s[h], \mathcal{T}_u[h])$ is given by the two integral formulas. Since $\operatorname{Lip}(\tilde{g}_\sigma) < \varepsilon$ and $\|e^{A_s t}\| \leq K e^{-\alpha t}$ for $t \geq 0$ (with $\alpha > 0$ depending on the spectral gap of $A_s$) and $\|e^{-A_u t}\| \leq K e^{-\beta t}$ for $t \geq 0$ (with $\beta > 0$ depending on the spectral gap of $A_u$), a standard estimate gives \begin{align*} \|\mathcal{T}[h_1] - \mathcal{T}[h_2]\|_{\mathrm{Lip}} \leq \frac{K\varepsilon}{\min(\alpha, \beta)} \cdot (1 + \delta) \cdot \|h_1 - h_2\|_{\mathrm{Lip}}. \end{align*} Choosing $\varepsilon$ (i.e., $\rho$) small enough that $K\varepsilon(1 + \delta)/\min(\alpha, \beta) < 1$, the operator $\mathcal{T}$ is a contraction on $\mathcal{X}_\delta$. By the [Banach Fixed-Point Theorem](/theorems/???), there exists a unique $h^* \in \mathcal{X}_\delta$ with $\mathcal{T}[h^*] = h^*$. The graph of $h^*$ is the local centre manifold $W^c_{\mathrm{loc}}$. [guided] The key idea is to convert the invariance PDE into an integral equation using the variation-of-constants formula. Rather than solving the PDE for $h$ directly, we reformulate the problem: given a candidate function $h$, the dynamics on the centre manifold is governed by the reduced ODE $\dot{x}_c = A_c x_c + \tilde{g}_c(h(x_c), x_c)$. Let $x_c(t; y)$ denote the solution with $x_c(0; y) = y$. **Deriving the stable integral formula.** For the trajectory on the manifold to be well-behaved as $t \to -\infty$ in the $x_s$-direction, we use the variation-of-constants formula for $x_s$: $x_s(t) = e^{A_s t} h_s(y) + \int_0^t e^{A_s(t-\tau)} \tilde{g}_s(\cdots) \, d\mathcal{L}^1(\tau)$. As $t \to -\infty$, the term $e^{A_s t} h_s(y)$ grows (since $e^{A_s t}$ grows in reverse time) unless $h_s(y)$ is chosen to cancel the growing mode. The unique choice that keeps $x_s$ bounded as $t \to -\infty$ is: \begin{align*} h_s(y) = \int_{-\infty}^{0} e^{-A_s \tau}\, \tilde{g}_s\bigl(h(x_c(\tau; y)),\, x_c(\tau; y)\bigr)\, d\mathcal{L}^1(\tau). \end{align*} This integral converges because $e^{-A_s \tau} = e^{A_s(-\tau)}$ decays exponentially as $\tau \to -\infty$ (since $A_s$ has eigenvalues with negative real part). Similarly, the unstable component $h_u$ is determined by integrating forward to $+\infty$, exploiting the decay of $e^{-A_u \tau}$ as $\tau \to +\infty$. **The contraction estimate.** Define the operator $\mathcal{T}$ by these integral formulas. For two candidates $h_1, h_2 \in \mathcal{X}_\delta$, the difference $\mathcal{T}[h_1] - \mathcal{T}[h_2]$ involves differences $\tilde{g}_s(h_1(\cdot), \cdot) - \tilde{g}_s(h_2(\cdot), \cdot)$, which are bounded by $\operatorname{Lip}(\tilde{g}_s) \cdot \|h_1 - h_2\|$. The exponential decay of $e^{A_s t}$ (with rate $\alpha$) means the integral gains a factor $K/\alpha$. Combining: \begin{align*} \|\mathcal{T}[h_1] - \mathcal{T}[h_2]\|_{\mathrm{Lip}} \leq \frac{K\varepsilon}{\min(\alpha, \beta)} \cdot (1 + \delta) \cdot \|h_1 - h_2\|_{\mathrm{Lip}}. \end{align*} The factor $(1 + \delta)$ accounts for the dependence of the flow $x_c(t; y)$ on $h$ through the reduced equation. Choosing $\varepsilon$ (equivalently, $\rho$) small enough that $K\varepsilon(1 + \delta)/\min(\alpha, \beta) < 1$ makes $\mathcal{T}$ a contraction. **The Banach space.** Define $\mathcal{X} = \{h \in C^0(E^c; E^s \oplus E^u) : h(\mathbf{0}) = \mathbf{0},\; \|h\|_{\mathrm{Lip}} < \infty\}$ with the Lipschitz norm, and let $\mathcal{X}_\delta = \{h \in \mathcal{X} : \|h\|_{\mathrm{Lip}} \leq \delta\}$ for small $\delta > 0$. The integral formulas above define the operator $\mathcal{T}: \mathcal{X}_\delta \to \mathcal{X}$ by $\mathcal{T}[h] = (\mathcal{T}_s[h], \mathcal{T}_u[h])$. **Why Lipschitz functions?** The space $\mathcal{X}_\delta$ consists of Lipschitz functions (not just continuous ones) because the integral formula involves composing $h$ with the flow $x_c(t; y)$. By the chain rule for Lipschitz functions, $\operatorname{Lip}(h \circ x_c(t; \cdot)) \leq \operatorname{Lip}(h) \cdot \operatorname{Lip}(x_c(t; \cdot))$, so controlling the Lipschitz constant of $\mathcal{T}[h]$ requires knowing $\operatorname{Lip}(h)$, which is the norm in $\mathcal{X}_\delta$. The [Banach Fixed-Point Theorem](/theorems/???) now applies: $\mathcal{T}$ maps $\mathcal{X}_\delta$ to itself (verified by a norm estimate similar to the contraction estimate) and is a contraction with constant $K\varepsilon(1+\delta)/\min(\alpha,\beta) < 1$. There exists a unique $h^* \in \mathcal{X}_\delta$ with $\mathcal{T}[h^*] = h^*$. The graph of $h^*$ is the local centre manifold $W^c_{\mathrm{loc}}$. [/guided] [/step] [step:Verify the tangency conditions and $C^r$ regularity] The fixed point $h^*$ satisfies $h^*(\mathbf{0}) = \mathbf{0}$ by construction (since $\mathbf{0}$ is in $\mathcal{X}_\delta$, which requires $h(\mathbf{0}) = \mathbf{0}$). For the tangency condition $Dh^*_{\mathbf{0}} = 0$, differentiate the fixed-point equation $\mathcal{T}[h^*] = h^*$ at $y = \mathbf{0}$. Since $\tilde{g}_\sigma(\mathbf{0}) = \mathbf{0}$ and $D(\tilde{g}_\sigma)_{\mathbf{0}} = 0$, the linearisation of $\mathcal{T}$ at $h = 0$ is the zero operator, which gives $Dh^*_{\mathbf{0}} = 0$. This confirms $T_{\mathbf{0}} W^c = E^c$. For the $C^r$ regularity, one bootstraps: the contraction argument above produces a Lipschitz $h^*$. Differentiating the fixed-point equation and applying the implicit function theorem (or repeating the contraction argument in the space of $C^1$ functions with the $C^1$ norm), one obtains $h^* \in C^1$. Iterating, each additional derivative of $\tilde{g}_\sigma$ (available because $f$ is $C^r$) yields one more derivative of $h^*$, up to $C^r$. The details of this regularity bootstrap are standard (see, e.g., the fibre contraction theorem of Hirsch, Pugh, and Shub): at each stage, the derivative $D^k h^*$ satisfies a linear integral equation whose coefficients depend on the lower-order derivatives $D^j h^*$ ($j < k$), already known to be continuous. Contractivity of the linearised operator then gives continuity of $D^k h^*$. [guided] **Tangency: $Dh^*_{\mathbf{0}} = 0$.** The tangency condition guarantees that the centre manifold $W^c$ is tangent to $E^c$ at the origin, not merely passing through it. To verify it, we differentiate the fixed-point equation $\mathcal{T}[h^*] = h^*$ at $y = \mathbf{0}$. The integral formula for $(\mathcal{T}_s[h])(y)$ involves $\tilde{g}_s(h(x_c(\tau; y)), x_c(\tau; y))$, which vanishes to second order at the origin since $\tilde{g}_s(\mathbf{0}) = \mathbf{0}$ and $D(\tilde{g}_s)_{\mathbf{0}} = 0$. When we differentiate with respect to $y$ and evaluate at $y = \mathbf{0}$, every term in the result contains a factor of either $D(\tilde{g}_s)_{\mathbf{0}} = 0$ or $Dh^*_{\mathbf{0}}$ (which appears linearly through the chain rule). The resulting equation for $Dh^*_{\mathbf{0}}$ has the form $Dh^*_{\mathbf{0}} = L(Dh^*_{\mathbf{0}})$ where $L$ is a linear operator. Since $\mathcal{T}$ is a contraction with contraction constant $< 1$, the linearised operator $L$ has norm $< 1$. The only solution to $M = L(M)$ with $\|L\| < 1$ is $M = 0$. Therefore $Dh^*_{\mathbf{0}} = 0$, confirming $T_{\mathbf{0}} W^c = E^c$. **$C^r$ regularity via the fibre contraction principle.** The Banach fixed-point argument above produces $h^* \in C^{0,1}$ (Lipschitz continuous). To obtain higher regularity, one differentiates the fixed-point equation $\mathcal{T}[h^*] = h^*$ to derive an equation for $Dh^*$: \begin{align*} Dh^* = D\mathcal{T}_{h^*}[Dh^*], \end{align*} which is itself a fixed-point problem, now in a space of continuous matrix-valued functions $E^c \to \mathrm{Hom}(E^c, E^s \oplus E^u)$. The contraction property of this linearised operator carries over from the original contraction because the spectral gap $\alpha$ (from the exponential decay of $e^{A_s t}$) still provides the necessary convergence of the integral operator. The Banach fixed-point theorem gives a unique continuous solution $Dh^*$, proving $h^* \in C^1$. This process -- called the **fibre contraction principle** (due to Hirsch, Pugh, and Shub) -- iterates: at the $k$-th stage, $D^k h^*$ satisfies a linear integral equation whose coefficients depend on $D^j h^*$ for $j < k$ (already known to be continuous), and contractivity of the linearised operator gives continuity of $D^k h^*$. Since $f$ is $C^r$, the modified nonlinearities $\tilde{g}_\sigma$ are also $C^r$, and the bootstrap can be iterated $r$ times, yielding $h^* \in C^r$. [/guided] [/step] [step:Construct $W^s$ and $W^u$ via the Stable Manifold Theorem, and verify all claimed properties] The stable and unstable manifolds $W^s$ and $W^u$ are constructed by applying the [Stable Manifold Theorem](/theorems/2778) to the system. The fixed point $\mathbf{0}$ may not be hyperbolic for the full system, but the subsystem restricted to $E^s \oplus E^u$ (with the centre directions projected out) is hyperbolic. The Stable Manifold Theorem produces unique local manifolds $W^s_{\mathrm{loc}}$ and $W^u_{\mathrm{loc}}$ tangent to $E^s$ and $E^u$ respectively, each of class $C^r$, and each locally invariant under the flow. Their dimensions are $n_s$ and $n_u$ by construction. Alternatively, $W^s$ and $W^u$ can be constructed by the same graph-transform/fixed-point argument used for $W^c$: seek $W^s$ as the graph of a function $E^s \to E^u \oplus E^c$ and solve the analogous integral equation. The exponential decay of $e^{A_s t}$ as $t \to +\infty$ and the exponential decay of $e^{-A_u \tau}$ as $\tau \to +\infty$ provide the necessary contractivity, and the construction proceeds identically to the centre manifold case. We now verify all four claimed properties: 1. **Dimensions:** $\dim W^s = n_s$, $\dim W^u = n_u$, and $\dim W^c = n_c$ by construction — each manifold is the graph of a function defined on the corresponding linear subspace. 2. **Tangency:** $T_{\mathbf{0}} W^s = E^s$, $T_{\mathbf{0}} W^u = E^u$, $T_{\mathbf{0}} W^c = E^c$ because $Dh^*_{\mathbf{0}} = 0$ in each case. 3. **Uniqueness of $W^s$ and $W^u$; non-uniqueness of $W^c$:** The stable and unstable manifolds are unique in a sufficiently small neighbourhood because the contraction mapping principle produces a unique fixed point. For the centre manifold, uniqueness can fail because the integral equation over $(-\infty, 0]$ for $h_s$ depends on the behaviour of trajectories as $t \to -\infty$ on the centre manifold, and different solutions of the modified system (with different cut-off radii $\rho$) can yield different global invariant manifolds that agree only to finite order at the origin. A classical counterexample is $\dot{x} = -x^3$, $\dot{y} = -y$ on $\mathbb{R}^2$, where $E^c = \{y = 0\}$ and $E^s = \{x = 0\}$. The centre manifold is the graph of a function $h: \mathbb{R} \to \mathbb{R}$ with $h(0) = 0$ and $h'(0) = 0$. Any function $h(x) = c\, e^{-1/(2x^2)}$ for $c \in \mathbb{R}$ satisfies the invariance equation and the tangency conditions, yielding a one-parameter family of distinct centre manifolds. All are flat at the origin (all derivatives vanish), which is why they agree to infinite order but are not equal. 4. **Invariance under the flow:** By construction, each manifold is the graph of a function satisfying the invariance equation: the ODE restricted to the graph stays on the graph. Hence if $x(0) \in W^\sigma$, then $x(t) \in W^\sigma$ for all $t$ in the maximal interval of existence of the trajectory. This completes the proof of all four parts of the Centre Manifold Theorem. [guided] The stable and unstable manifolds $W^s$ and $W^u$ are constructed by applying the [Stable Manifold Theorem](/theorems/2778) to the full system. Although $\mathbf{0}$ is non-hyperbolic, the subsystem on $E^s \oplus E^u$ is hyperbolic, and the theorem produces unique local manifolds $W^s_{\mathrm{loc}}$ tangent to $E^s$ and $W^u_{\mathrm{loc}}$ tangent to $E^u$, each $C^r$ and locally invariant. We verify the four claimed properties: (1) dimensions match the respective subspaces by construction; (2) tangency follows from $Dh^*_{\mathbf{0}} = 0$ in each case; (3) $W^s$ and $W^u$ are unique while $W^c$ may not be; (4) invariance holds because each manifold is the graph of a function satisfying the invariance equation. Let us address the most surprising aspect: **why is $W^c$ non-unique while $W^s$ and $W^u$ are unique?** The fundamental difference is the exponential rate in the integral formula. For the stable manifold, the function $h_s$ is determined by integrating forward to $+\infty$: the exponential decay $\|e^{A_s t}\| \leq Ke^{-\alpha t}$ as $t \to +\infty$ kills all ambiguity. The integral $\int_0^\infty e^{-A_s \tau} (\cdots) \, d\mathcal{L}^1(\tau)$ converges to a unique value regardless of the global behaviour of the solution, because the exponential decay suppresses any differences. Similarly, $W^u$ involves integration to $-\infty$ with $\|e^{-A_u \tau}\| \leq Ke^{-\beta \tau}$ as $\tau \to +\infty$. For the centre manifold, the eigenvalues of $A_c$ have zero real part, so $e^{A_c t}$ does not decay in either time direction -- it merely oscillates or grows polynomially (for non-trivial Jordan blocks). The contraction in the fixed-point argument comes entirely from the smallness of $\varepsilon = \operatorname{Lip}(\tilde{g}_\sigma)$, not from exponential decay. Different choices of cut-off radius $\rho$ produce different modified systems, and the global fixed points of the corresponding operators $\mathcal{T}$ may differ. These different global solutions restrict to different manifolds near the origin. The classical counterexample illustrates this concretely. Consider $\dot{x} = -x^3$, $\dot{y} = -y$ on $\mathbb{R}^2$, with $E^c = \{y = 0\}$ and $E^s = \{x = 0\}$. Along the $x$-axis, solutions $\dot{x} = -x^3$ give $x(t) = x_0 / \sqrt{1 + 2x_0^2 t}$, which decays only algebraically (not exponentially). The $y$-component decays exponentially: $y(t) = y_0 e^{-t}$. The function $h(x) = c\,e^{-1/(2x^2)}$ for any $c \in \mathbb{R}$ satisfies the invariance equation and the tangency conditions $h(0) = 0$, $h'(0) = 0$ (since $e^{-1/(2x^2)}$ and all its derivatives vanish at $x = 0$). This gives a one-parameter family of distinct centre manifolds -- all tangent to the $x$-axis at the origin, all agreeing to infinite order, but genuinely different as sets. [/guided] [/step]