[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 $\mathbf{x} = (\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c)$ with $\mathbf{x}_s \in E^s \cong \mathbb{R}^{n_s}$, $\mathbf{x}_u \in E^u \cong \mathbb{R}^{n_u}$, $\mathbf{x}_c \in E^c \cong \mathbb{R}^{n_c}$, the system $\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x})$ takes the form \begin{align*} \dot{\mathbf{x}}_s &= A_s \mathbf{x}_s + \mathbf{g}_s(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c), \\ \dot{\mathbf{x}}_u &= A_u \mathbf{x}_u + \mathbf{g}_u(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c), \\ \dot{\mathbf{x}}_c &= A_c \mathbf{x}_c + \mathbf{g}_c(\mathbf{x}_s, \mathbf{x}_u, \mathbf{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 $\mathbf{g}_s, \mathbf{g}_u, \mathbf{g}_c$ are $C^r$ functions satisfying \begin{align*} \mathbf{g}_\sigma(\mathbf{0}) = \mathbf{0}, \qquad D(\mathbf{g}_\sigma)_{\mathbf{0}} = 0, \qquad \sigma \in \{s, u, c\}. \end{align*} The vanishing of $\mathbf{g}_\sigma$ and its derivative at the origin follows from $\mathbf{f}(\mathbf{0}) = \mathbf{0}$ and the definition of $A = Jf_{\mathbf{0}}$.[/step]

[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 $\mathbf{x} = (\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c)$, the system $\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x})$ splits as \begin{align*} \dot{\mathbf{x}}_s &= A_s \mathbf{x}_s + \mathbf{g}_s(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c), \\ \dot{\mathbf{x}}_u &= A_u \mathbf{x}_u + \mathbf{g}_u(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c), \\ \dot{\mathbf{x}}_c &= A_c \mathbf{x}_c + \mathbf{g}_c(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c), \end{align*} where the nonlinear terms $\mathbf{g}_\sigma$ satisfy $\mathbf{g}_\sigma(\mathbf{0}) = \mathbf{0}$ and $D(\mathbf{g}_\sigma)_{\mathbf{0}} = 0$. These conditions say that $\mathbf{g}_\sigma(\mathbf{x}) = O(|\mathbf{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 $A\mathbf{x}$ from $\mathbf{f}(\mathbf{x}) = A\mathbf{x} + (\mathbf{f}(\mathbf{x}) - A\mathbf{x})$, using $\mathbf{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: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, \\ \mathbf{x}_c &\mapsto h(\mathbf{x}_c) = (h_s(\mathbf{x}_c),\, h_u(\mathbf{x}_c)), \end{align*} where $B_\rho^c = \{\mathbf{x}_c \in E^c : |\mathbf{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}} = \{(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c) : (\mathbf{x}_s, \mathbf{x}_u) = h(\mathbf{x}_c),\; \mathbf{x}_c \in B_\rho^c\}$. For $W^c_{\mathrm{loc}}$ to be invariant under the flow, the trajectory starting at any point $(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c) = (h_s(\mathbf{x}_c), h_u(\mathbf{x}_c), \mathbf{x}_c)$ on the graph must remain on the graph. Differentiating the relation $\mathbf{x}_s(t) = h_s(\mathbf{x}_c(t))$ with respect to $t$ and using the chain rule yields \begin{align*} \dot{\mathbf{x}}_s = Dh_{s,\mathbf{x}_c}\, \dot{\mathbf{x}}_c. \end{align*} Substituting the ODE for $\dot{\mathbf{x}}_s$ and $\dot{\mathbf{x}}_c$ from the decomposed system, and evaluating the nonlinear terms on the graph $(\mathbf{x}_s, \mathbf{x}_u) = h(\mathbf{x}_c)$, we obtain the **invariance equation** for $h_s$: \begin{align*} A_s\, h_s(\mathbf{x}_c) + \mathbf{g}_s(h(\mathbf{x}_c), \mathbf{x}_c) = Dh_{s,\mathbf{x}_c}\!\left[A_c\, \mathbf{x}_c + \mathbf{g}_c(h(\mathbf{x}_c), \mathbf{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)$.[/step]

[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 $(\mathbf{x}_s, \mathbf{x}_u) = h(\mathbf{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}} = \{(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c) : (\mathbf{x}_s, \mathbf{x}_u) = h(\mathbf{x}_c),\; \mathbf{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 $\mathbf{x}_s(t) = h_s(\mathbf{x}_c(t))$ must hold for all $t$. Differentiating both sides with respect to $t$ gives, by the chain rule: \begin{align*} \dot{\mathbf{x}}_s = Dh_{s,\mathbf{x}_c}\, \dot{\mathbf{x}}_c. \end{align*} The left-hand side is the ODE for $\mathbf{x}_s$: $\dot{\mathbf{x}}_s = A_s \mathbf{x}_s + \mathbf{g}_s(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c)$. The right-hand side uses the ODE for $\mathbf{x}_c$: $\dot{\mathbf{x}}_c = A_c \mathbf{x}_c + \mathbf{g}_c(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c)$. Evaluating everything on the graph (replacing $(\mathbf{x}_s, \mathbf{x}_u)$ by $h(\mathbf{x}_c)$) and equating gives the **invariance equation** for $h_s$: \begin{align*} A_s\, h_s(\mathbf{x}_c) + \mathbf{g}_s(h(\mathbf{x}_c), \mathbf{x}_c) = Dh_{s,\mathbf{x}_c}\!\left[A_c\, \mathbf{x}_c + \mathbf{g}_c(h(\mathbf{x}_c), \mathbf{x}_c)\right]. \end{align*} An analogous equation holds for $h_u$. This is a **quasilinear PDE** for $h$: the highest derivative $Dh_{s,\mathbf{x}_c}$ appears multiplied by terms that depend on $h$ itself (through $\mathbf{g}_c(h(\mathbf{x}_c), \mathbf{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:Cut off the nonlinearity to obtain a globally Lipschitz system]The nonlinear terms $\mathbf{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{\mathbf{g}}_\sigma(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c) := \chi\!\left(\frac{|(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c)|}{\rho}\right) \mathbf{g}_\sigma(\mathbf{x}_s, \mathbf{x}_u, \mathbf{x}_c), \qquad \sigma \in \{s, u, c\}. \end{align*} The modified system $\dot{\mathbf{x}} = A\mathbf{x} + \tilde{\mathbf{g}}(\mathbf{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(\mathbf{g}_\sigma)_{\mathbf{0}} = 0$, the Lipschitz constant of $\tilde{\mathbf{g}}_\sigma$ can be made arbitrarily small by taking $\rho$ sufficiently small: \begin{align*} \operatorname{Lip}(\tilde{\mathbf{g}}_\sigma) \leq C\rho \end{align*} for a constant $C > 0$ depending on the $C^2$ norm of $\mathbf{g}_\sigma$ near the origin. We fix $\rho$ small enough that $\operatorname{Lip}(\tilde{\mathbf{g}}_\sigma) < \varepsilon$ for a threshold $\varepsilon > 0$ to be determined below.[/step]

[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 $\mathbf{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 $\mathbf{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{\mathbf{g}}_\sigma(\mathbf{x}) := \chi\!\left(\frac{|\mathbf{x}|}{\rho}\right) \mathbf{g}_\sigma(\mathbf{x}) \end{align*} agrees with $\mathbf{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{\mathbf{g}}_\sigma) \leq C\rho$, which can be made arbitrarily small by choosing $\rho$ small. Why does this hold? Since $D(\mathbf{g}_\sigma)_{\mathbf{0}} = 0$ and $\mathbf{g}_\sigma$ is $C^r$ ($r \geq 2$), the mean value inequality gives $|D(\mathbf{g}_\sigma)_{\mathbf{x}}| \leq C_1 |\mathbf{x}|$ for $|\mathbf{x}| \leq 2\rho$. When we differentiate the product $\tilde{\mathbf{g}}_\sigma = \chi \cdot \mathbf{g}_\sigma$, the product rule gives $D\tilde{\mathbf{g}}_\sigma = D\chi \cdot \mathbf{g}_\sigma + \chi \cdot D\mathbf{g}_\sigma$. The first term is bounded by $|D\chi| \cdot |\mathbf{g}_\sigma| = O(1/\rho) \cdot O(|\mathbf{x}|^2) = O(\rho)$ (since $|\mathbf{x}| \leq 2\rho$ on the support). The second term is bounded by $|\chi| \cdot |D\mathbf{g}_\sigma| = 1 \cdot O(|\mathbf{x}|) = O(\rho)$. Combining: $|D\tilde{\mathbf{g}}_\sigma| = O(\rho)$ everywhere, so $\operatorname{Lip}(\tilde{\mathbf{g}}_\sigma) \leq C\rho$. The modified system $\dot{\mathbf{x}} = A\mathbf{x} + \tilde{\mathbf{g}}(\mathbf{x})$ is identical to the original inside $B(\mathbf{0}, \rho)$ and is purely linear ($\dot{\mathbf{x}} = A\mathbf{x}$) 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{\mathbf{g}}_\sigma) < \varepsilon$ for a threshold $\varepsilon > 0$ determined by the contraction condition in the next step.[/guided]

[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 $\mathbf{x}_s$-equation forward in time along trajectories on the manifold. On the centre manifold, the dynamics is governed by \begin{align*} \dot{\mathbf{x}}_c = A_c\, \mathbf{x}_c + \tilde{\mathbf{g}}_c(h(\mathbf{x}_c), \mathbf{x}_c). \end{align*} Denote by $\mathbf{x}_c(t; \mathbf{y})$ the solution of this reduced equation with $\mathbf{x}_c(0; \mathbf{y}) = \mathbf{y}$. The stable component must satisfy \begin{align*} h_s(\mathbf{y}) = \int_{-\infty}^{0} e^{-A_s \tau}\, \tilde{\mathbf{g}}_s\bigl(h(\mathbf{x}_c(\tau; \mathbf{y})),\, \mathbf{x}_c(\tau; \mathbf{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{\mathbf{g}}_s) \cdot \sup |\mathbf{x}_c(\tau; \mathbf{y})|$. Similarly, for the unstable component, integrate the $\mathbf{x}_u$-equation backward: \begin{align*} h_u(\mathbf{y}) = -\int_{0}^{\infty} e^{-A_u \tau}\, \tilde{\mathbf{g}}_u\bigl(h(\mathbf{x}_c(\tau; \mathbf{y})),\, \mathbf{x}_c(\tau; \mathbf{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_{\mathbf{y} \neq \mathbf{y}'} |h(\mathbf{y}) - h(\mathbf{y}')| / |\mathbf{y} - \mathbf{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{\mathbf{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}}$.[/step]

[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{\mathbf{x}}_c = A_c \mathbf{x}_c + \tilde{\mathbf{g}}_c(h(\mathbf{x}_c), \mathbf{x}_c)$. Let $\mathbf{x}_c(t; \mathbf{y})$ denote the solution with $\mathbf{x}_c(0; \mathbf{y}) = \mathbf{y}$. **Deriving the stable integral formula.** For the trajectory on the manifold to be well-behaved as $t \to -\infty$ in the $\mathbf{x}_s$-direction, we use the variation-of-constants formula for $\mathbf{x}_s$: $\mathbf{x}_s(t) = e^{A_s t} h_s(\mathbf{y}) + \int_0^t e^{A_s(t-\tau)} \tilde{\mathbf{g}}_s(\cdots) \, d\mathcal{L}^1(\tau)$. As $t \to -\infty$, the term $e^{A_s t} h_s(\mathbf{y})$ grows (since $e^{A_s t}$ grows in reverse time) unless $h_s(\mathbf{y})$ is chosen to cancel the growing mode. The unique choice that keeps $\mathbf{x}_s$ bounded as $t \to -\infty$ is: \begin{align*} h_s(\mathbf{y}) = \int_{-\infty}^{0} e^{-A_s \tau}\, \tilde{\mathbf{g}}_s\bigl(h(\mathbf{x}_c(\tau; \mathbf{y})),\, \mathbf{x}_c(\tau; \mathbf{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{\mathbf{g}}_s(h_1(\cdot), \cdot) - \tilde{\mathbf{g}}_s(h_2(\cdot), \cdot)$, which are bounded by $\operatorname{Lip}(\tilde{\mathbf{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 $\mathbf{x}_c(t; \mathbf{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 $\mathbf{x}_c(t; \mathbf{y})$. By the chain rule for Lipschitz functions, $\operatorname{Lip}(h \circ \mathbf{x}_c(t; \cdot)) \leq \operatorname{Lip}(h) \cdot \operatorname{Lip}(\mathbf{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: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 $\mathbf{y} = \mathbf{0}$. Since $\tilde{\mathbf{g}}_\sigma(\mathbf{0}) = \mathbf{0}$ and $D(\tilde{\mathbf{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{\mathbf{g}}_\sigma$ (available because $\mathbf{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^*$.[/step]

[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 $\mathbf{y} = \mathbf{0}$. The integral formula for $(\mathcal{T}_s[h])(\mathbf{y})$ involves $\tilde{\mathbf{g}}_s(h(\mathbf{x}_c(\tau; \mathbf{y})), \mathbf{x}_c(\tau; \mathbf{y}))$, which vanishes to second order at the origin since $\tilde{\mathbf{g}}_s(\mathbf{0}) = \mathbf{0}$ and $D(\tilde{\mathbf{g}}_s)_{\mathbf{0}} = 0$. When we differentiate with respect to $\mathbf{y}$ and evaluate at $\mathbf{y} = \mathbf{0}$, every term in the result contains a factor of either $D(\tilde{\mathbf{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 $\mathbf{f}$ is $C^r$, the modified nonlinearities $\tilde{\mathbf{g}}_\sigma$ are also $C^r$, and the bootstrap can be iterated $r$ times, yielding $h^* \in C^r$.[/guided]

[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 $\mathbf{x}(0) \in W^\sigma$, then $\mathbf{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.[/step]

[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{\mathbf{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]