I'm trying to prove that the maximal number of linearly independent Killing vector fields on a Riemannian manifold $(M,g)$ of dimension $N$ is $N(N+1)/2$.
Here's what I got so far. By definition a Killing field $\xi$ satisfies $\mathcal L_\xi g=0$, which can be written as $\nabla_i\xi_j + \nabla_j\xi_i=0$ in a chart. This then implies that $\nabla_i\nabla_j\xi_k = -R^l{}_{ijk}\xi_l$ in terms of the curvature tensor.
I get the impression that the remaining part of the proof comes down to the fact that this second order PDE, which is essentially of the form $$ \partial_i\partial_j\xi_m = \phi^{ml}_{ij}\xi_k + \phi^{mkl}_{ij}\partial_k\xi_l$$ with $\phi^{mk}_{ij},\phi^{mkl}_{ij}$ smooth functions on an open subset or $\mathbb R^N$, has at most one solution when initial conditions $\xi_i(p),\partial_i\xi_j(p)$ are specified at some point $p$. If so, then the proof is complete, because there are $N$ linearly independent possibilities for $\xi_j(p)$ and $N(N-1)/2$ for $\partial_i\xi_j(p)$ (as $\partial_i\xi_j(p)$ is related in a certain way to $\partial_j\xi_i(p)$ by Killing's equation), so that the linear space of initial conditions has dimension $N + N(N-1)/2 = N(N+1)/2$.
Assuming for the moment that this approach is correct, how do we know that this PDE has at most one solution given such initial conditions? And is the degree of smoothness essential here? I'm mainly interested in $C^\infty$ manifolds, but if the result generalizes for instance two $C^2$ manifolds, which I can imagine, that would be useful.
EDIT:
Maybe it is helpful to notice that if we assume our Killing vectors to be real-analytic, then uniqueness of the solution to the PDE is not hard to prove. In this case, if $\xi(p)$ and $\partial \xi(p)$ are known (I'm omitting indices for the moment) then the second order PDE tells us the value of $\partial\partial\xi (p)$. Differentiating the PDE now yields a formula for $\partial\partial\partial\xi$ in terms of lower order derivatives, which we already know. Continuing like this we obtain all order derivatives of $\xi$ at $p$, which uniquely fixes the Tayler series of $\xi$ around $p$ and by real-analiticity this implies that $\xi$, locally being equal to its Taylor series, is uniquely determined.
Unfortunately this argument is based heavily on the supposition that $\xi$ is real-analytic. A proof for $\xi$ merely being smooth would probably need to take en entirely different route.
I found the proof of uniqueness in Wald's General Relativity (appendix C).
The system of PDEs $\nabla_i\nabla_j\xi_k = -R^l_{ijk}\xi_l$ is equivalent to the system of PDEs
$$ \nabla_iL_{jk} = -R^l_{ijk}\xi_l,\qquad \nabla_j\xi_k = L_{jk}, $$
which again is equivalent to the statement that for any vector field $v$
$$ v^i\nabla_iL_{jk} = -R^l_{ijk}v^i\xi_l,\qquad v^j\nabla_j\xi_k = L_{jk}v^j. $$
If we choose $v$ to be a tangent field along a curve $\gamma$ starting at $p\in M$ (at which we provide initial conditions), then along this curve the equations reduce to the system of ODEs
$$\frac{\text{d}}{\text{d} t} L_{jk} - \Gamma^l_{ij}L_{lk}v^i - \Gamma^l_{ik}L_{jl}v^i= -R^l_{ijk}v^i\xi_l,\qquad \frac{\text{d}}{\text{d} t}\xi_k - \Gamma^l_{jk}v^j\xi_l = L_{jk}v^j.$$
By the well-known existence and uniqueness theorem for ODEs, this system admits at most one solution when initial conditions for $\xi$ and $L$ are provided at $p$. From this it follows that our original system of PDEs also can have at most one solution.
Concering smoothness assumptions, in order for the proposition stated in the question to hold it is sufficient that $M$ is a $C^2$ manifold and that Killing fields are also required to be $C^2$. But note that if the Killing field is $C^2$ then the system of second order PDEs implies that all second derivatives are in fact $C^1$, so that $\xi$ is $C^3$. But then the PDEs imply that the second derivatives are $C^2$, so $\xi$ is $C^4$, ad infinitum, i.e., any $C^2$ Killing field is necessarily $C^\infty$.