[proofplan]
We construct the local stable manifold $W^s_{\mathrm{loc}}(0)$ as the graph of a $C^r$ function from a neighbourhood of $0$ in the stable subspace $E^s$ to the unstable subspace $E^u$, using the Banach fixed-point theorem in a space of Lipschitz graph maps. The key steps are: (1) write the ODE in coordinates adapted to the $E^s \oplus E^u$ decomposition, (2) reformulate the condition "$\varphi(t, x) \to 0$ as $t \to +\infty$" as an integral equation on the space of exponentially decaying trajectories, (3) solve this integral equation via contraction, and (4) verify that the solution set is a $C^r$ manifold tangent to $E^s$ at the origin. The unstable manifold $W^u_{\mathrm{loc}}(0)$ is obtained by reversing time.
[/proofplan]
[step:Set up coordinates adapted to the stable-unstable decomposition]
Without loss of generality, $f(0) = 0$ (the fixed point is at the origin). Since $f$ is $C^r$ with $r \geq 1$, write
\begin{align*}
f(x) = Ax + g(x),
\end{align*}
where $A := Jf_0$ and $g(x) = O(|x|^2)$ satisfies $g(0) = 0$, $Jg_0 = 0$.
Since $0$ is hyperbolic, the eigenvalues of $A$ split into $n^s$ eigenvalues with $\operatorname{Re}(\lambda) < 0$ and $n^u$ eigenvalues with $\operatorname{Re}(\lambda) > 0$, with $n^s + n^u = n$. The corresponding generalised eigenspaces give the direct-sum decomposition $\mathbb{R}^n = E^s \oplus E^u$. In coordinates $x = (x_s, x_u)$ with $x_s \in E^s \cong \mathbb{R}^{n^s}$ and $x_u \in E^u \cong \mathbb{R}^{n^u}$, the matrix $A$ is block-diagonal:
\begin{align*}
A = \begin{pmatrix} A_s & 0 \\ 0 & A_u \end{pmatrix},
\end{align*}
where $A_s \in \mathbb{R}^{n^s \times n^s}$ has all eigenvalues with negative real part and $A_u \in \mathbb{R}^{n^u \times n^u}$ has all eigenvalues with positive real part. Choose constants $\alpha, \beta > 0$ with $\operatorname{Re}(\lambda) \leq -\alpha$ for eigenvalues of $A_s$ and $\operatorname{Re}(\lambda) \geq \beta$ for eigenvalues of $A_u$. Using a norm adapted to the Jordan structure, there exists $C > 0$ such that
\begin{align*}
|e^{A_s t}| &\leq C e^{-\alpha t} \quad \text{for } t \geq 0, \\
|e^{A_u t}| &\leq C e^{\beta t} \quad \text{for } t \leq 0 \quad (\text{equivalently } |e^{-A_u t}| \leq C e^{-\beta t} \text{ for } t \geq 0).
\end{align*}
[/step]
[step:Characterise stable manifold trajectories via an integral equation]
Write $g = (g_s, g_u)$ in the stable-unstable coordinates. The system $\dot{x} = Ax + g(x)$ becomes
\begin{align*}
\dot{x}_s &= A_s x_s + g_s(x_s, x_u), \\
\dot{x}_u &= A_u x_u + g_u(x_s, x_u).
\end{align*}
A trajectory $x(t) = (x_s(t), x_u(t))$ with $x(0) = (\xi, \eta)$ satisfies $\varphi(t, (\xi, \eta)) \to 0$ as $t \to +\infty$ if and only if both components decay. By the variation-of-constants formula:
\begin{align*}
x_s(t) &= e^{A_s t} \xi + \int_0^t e^{A_s(t-\tau)} g_s(x(\tau)) \, d\mathcal{L}^1(\tau), \\
x_u(t) &= e^{A_u t} \eta + \int_0^t e^{A_u(t-\tau)} g_u(x(\tau)) \, d\mathcal{L}^1(\tau).
\end{align*}
For the unstable component to decay as $t \to +\infty$, the growing term $e^{A_u t} \eta$ must be cancelled. Rewriting the unstable equation using the condition $x_u(t) \to 0$:
\begin{align*}
x_u(t) = -\int_t^{\infty} e^{A_u(t-\tau)} g_u(x(\tau)) \, d\mathcal{L}^1(\tau),
\end{align*}
where the integral converges because $e^{A_u(t - \tau)} = e^{-A_u(\tau - t)}$ decays exponentially for $\tau > t$ (since $A_u$ has eigenvalues with positive real part) and $g_u(x(\tau)) \to 0$ as $\tau \to +\infty$ (since $x(\tau) \to 0$ and $g(0) = 0$).
Setting $t = 0$ determines $\eta$ uniquely in terms of the trajectory:
\begin{align*}
\eta = -\int_0^{\infty} e^{-A_u \tau} g_u(x(\tau)) \, d\mathcal{L}^1(\tau).
\end{align*}
Thus the initial condition $\eta$ on $E^u$ is not free — it is determined by $\xi$ and the requirement of decay.
[guided]
The integral equation reformulation is the heart of the proof. We write the decomposed system $\dot{x}_s = A_s x_s + g_s(x)$, $\dot{x}_u = A_u x_u + g_u(x)$ and apply the variation-of-constants formula to each component separately. For a trajectory $x(t) = (x_s(t), x_u(t))$ with initial condition $x(0) = (\xi, \eta)$:
\begin{align*}
x_s(t) &= e^{A_s t} \xi + \int_0^t e^{A_s(t-\tau)} g_s(x(\tau)) \, d\mathcal{L}^1(\tau), \\
x_u(t) &= e^{A_u t} \eta + \int_0^t e^{A_u(t-\tau)} g_u(x(\tau)) \, d\mathcal{L}^1(\tau).
\end{align*}
The key observation is that the stable part $e^{A_s t}$ automatically decays as $t \to +\infty$ (since all eigenvalues of $A_s$ have negative real part), but the unstable part $e^{A_u t}$ grows exponentially. For the full trajectory to decay to $0$, both components must decay. The stable component decays automatically (the linear term $e^{A_s t} \xi$ decays, and the integral term is controlled by the smallness of $g$). The unstable component, however, has a growing linear term $e^{A_u t} \eta$ that must be cancelled.
How do we force the unstable component to decay? We impose a "boundary condition at $t = +\infty$": require $x_u(t) \to 0$ as $t \to +\infty$. Rearranging the variation-of-constants formula and sending $t \to +\infty$, the growing mode $e^{A_u t}$ forces a specific choice of $\eta$. Equivalently, we can write the unstable component by integrating backward from $+\infty$:
\begin{align*}
x_u(t) = -\int_t^{\infty} e^{A_u(t-\tau)} g_u(x(\tau)) \, d\mathcal{L}^1(\tau).
\end{align*}
This integral converges because $e^{A_u(t - \tau)} = e^{-A_u(\tau - t)}$ decays exponentially for $\tau > t$ (since $A_u$ has eigenvalues with positive real part, $e^{-A_u s}$ decays as $s \to +\infty$), and $g_u(x(\tau)) \to 0$ as $\tau \to +\infty$ (since $x(\tau) \to 0$ and $g(0) = 0$). Setting $t = 0$ determines $\eta$ uniquely:
\begin{align*}
\eta = -\int_0^{\infty} e^{-A_u \tau} g_u(x(\tau)) \, d\mathcal{L}^1(\tau).
\end{align*}
This means the stable manifold is parameterised by $\xi \in E^s$ alone: for each small $\xi$, the requirement that the trajectory decays to $0$ uniquely determines the $E^u$-component $\eta = \eta(\xi)$ of the initial condition. The graph $\{(\xi, \eta(\xi)) : \xi \in B_{E^s}(0, \delta)\}$ is the local stable manifold $W^s_{\mathrm{loc}}(0)$. The parameter $\eta$ is not free -- it is a function of $\xi$ and the entire future trajectory, computed as an integral to $+\infty$.
[/guided]
[/step]
[step:Solve the integral equation by contraction on a space of decaying paths]
Define the Banach space of exponentially decaying continuous paths:
\begin{align*}
\mathcal{Y}_\gamma := \left\{x \in C([0, \infty), \mathbb{R}^n) : \|x\|_\gamma := \sup_{t \geq 0} e^{\gamma t} |x(t)| < \infty\right\},
\end{align*}
where $\gamma > 0$ is chosen with $0 < \gamma < \alpha$ (so that $e^{-\gamma t}$ decays slower than $e^{-\alpha t}$ but still decays). For $\xi \in E^s$ with $|\xi|$ small, define the operator $\mathcal{T}_\xi: \mathcal{Y}_\gamma \to \mathcal{Y}_\gamma$ by
\begin{align*}
(\mathcal{T}_\xi x)_s(t) &:= e^{A_s t} \xi + \int_0^t e^{A_s(t - \tau)} g_s(x(\tau)) \, d\mathcal{L}^1(\tau), \\
(\mathcal{T}_\xi x)_u(t) &:= -\int_t^{\infty} e^{A_u(t - \tau)} g_u(x(\tau)) \, d\mathcal{L}^1(\tau).
\end{align*}
We verify that $\mathcal{T}_\xi$ maps a small ball in $\mathcal{Y}_\gamma$ to itself. Using the exponential bounds on $e^{A_s t}$ and $e^{-A_u t}$ together with the estimate $|g(x)| \leq \varepsilon |x|$ (valid for $|x| \leq \delta$ with $\varepsilon \to 0$ as $\delta \to 0$), the norms of the integral terms are bounded by
\begin{align*}
\sup_{t \geq 0} e^{\gamma t} \left|\int_0^t C e^{-\alpha(t-\tau)} \varepsilon |x(\tau)| \, d\mathcal{L}^1(\tau)\right| \leq \frac{C\varepsilon}{\alpha - \gamma} \|x\|_\gamma,
\end{align*}
and similarly for the unstable component with $\beta$ replacing $\alpha$. The contraction estimate follows from the Lipschitz property of $g$: for $x, y \in \mathcal{Y}_\gamma$,
\begin{align*}
\|\mathcal{T}_\xi x - \mathcal{T}_\xi y\|_\gamma \leq C\varepsilon \left(\frac{1}{\alpha - \gamma} + \frac{1}{\beta - \gamma}\right) \|x - y\|_\gamma.
\end{align*}
Choosing $\delta$ (and hence $\varepsilon$) small enough that $C\varepsilon((\alpha - \gamma)^{-1} + (\beta - \gamma)^{-1}) < 1$, the operator $\mathcal{T}_\xi$ is a contraction. By the [Banach Fixed-Point Theorem](/theorems/???), there exists a unique fixed point $x^\xi \in \mathcal{Y}_\gamma$.
[/step]
[step:Extract the manifold as a graph and verify tangency]
For each small $\xi \in E^s$, the unique fixed point $x^\xi$ of $\mathcal{T}_\xi$ gives a trajectory with $x^\xi_s(0) = \xi$ and $x^\xi_u(0) = -\int_0^\infty e^{-A_u \tau} g_u(x^\xi(\tau)) \, d\mathcal{L}^1(\tau) =: \sigma(\xi)$. Define
\begin{align*}
\sigma: B_{E^s}(0, \delta') &\to E^u \\
\xi &\mapsto -\int_0^{\infty} e^{-A_u \tau} g_u(x^\xi(\tau)) \, d\mathcal{L}^1(\tau).
\end{align*}
The local stable manifold is
\begin{align*}
W^s_{\mathrm{loc}}(0) = \{(\xi, \sigma(\xi)) : \xi \in B_{E^s}(0, \delta')\}.
\end{align*}
**Dimension:** $\dim W^s_{\mathrm{loc}}(0) = \dim E^s = n^s$, since $W^s_{\mathrm{loc}}(0)$ is the graph of a function from an $n^s$-dimensional domain. This verifies (i).
**Tangency:** Since $g(0) = 0$ and $Jg_0 = 0$, the map $\xi \mapsto x^\xi$ depends differentiably on $\xi$ (by the smooth dependence of fixed points on parameters in the Banach fixed-point theorem), and $J\sigma_0 = 0$. Therefore $T_0 W^s_{\mathrm{loc}}(0) = E^s \oplus \{0\} = E^s$, verifying (ii).
**Regularity:** Since $f$ is $C^r$ and the fixed-point map $\xi \mapsto x^\xi$ depends $C^r$-smoothly on the parameter $\xi$ (by the implicit function theorem in Banach spaces applied to the contraction), $\sigma$ is $C^r$.
**Asymptotic characterisation (iii):** By construction, $(\xi, \sigma(\xi)) \in W^s_{\mathrm{loc}}(0)$ if and only if $\varphi(t, (\xi, \sigma(\xi))) \to 0$ as $t \to +\infty$, with the decay rate $|\varphi(t, (\xi, \sigma(\xi)))| \leq \text{const} \cdot e^{-\gamma t}$.
**Invariance (iv):** If $x(0) \in W^s_{\mathrm{loc}}(0)$, then $x(t) = \varphi(t, x(0)) \to 0$ as $t \to +\infty$. For $t \geq 0$, the trajectory $\tau \mapsto \varphi(\tau, \varphi(t, x(0))) = \varphi(\tau + t, x(0))$ also decays to $0$, so $\varphi(t, x(0)) \in W^s_{\mathrm{loc}}(0)$ (provided it remains in the neighbourhood where the local manifold is defined). Hence $\varphi(t, W^s_{\mathrm{loc}}(0)) \subset W^s_{\mathrm{loc}}(0)$ for $t \geq 0$.
[/step]
[step:Obtain the unstable manifold by time reversal]
The system $\dot{x} = -f(x)$ has the same fixed point $0$, with Jacobian $-A$. The eigenvalues of $-A$ are $-\lambda_1, \ldots, -\lambda_n$, so the stable subspace of $-A$ is $E^u$ and the unstable subspace is $E^s$. Applying the stable manifold construction to $\dot{x} = -f(x)$ yields a $C^r$ manifold tangent to $E^u$ at $0$, consisting of points whose forward orbits under $-f$ converge to $0$ — equivalently, points whose backward orbits under $f$ converge to $0$. This is precisely $W^u_{\mathrm{loc}}(0)$. Properties (i)–(iv) for the unstable manifold follow from those of the stable manifold of the time-reversed system.
[/step]