[proofplan]
We compute the second-order part of the variation of the Ricci tensor in local coordinates, because the principal symbol only records the coefficients of the highest-order derivatives of the variation tensor $h$. The variation of the Christoffel symbols gives the second derivatives of $h$ appearing in $D\operatorname{Ric}_g(h)$; all curvature, Christoffel, and first-derivative terms are lower order and therefore do not enter the symbol. Replacing each covariant derivative $\nabla_i$ in the second-order part by the covector component $\xi_i$ gives the displayed principal symbol. Finally, substituting $h=\xi\otimes\eta+\eta\otimes\xi$ shows that these infinitesimal-diffeomorphism directions are annihilated by the symbol.
[/proofplan]
[step:Compute the second-order part of the linearized Ricci tensor]
Fix $p \in M$, and let $(U,x)$ be a smooth coordinate chart around $p$ with coordinate vector fields $\partial_{x_1},\dots,\partial_{x_n}$. Write $\Gamma_{ij,k}$ for the coefficient of $\partial_{x_k}$ in $\nabla_{\partial_{x_i}}\partial_{x_j}$. Let $h \in \Gamma(S^2T^*M)$ be a smooth symmetric $2$-tensor.
We use the Levi-Civita connection coefficient formula in the chart $(U,x)$ and differentiate it with respect to the metric parameter. Equivalently, differentiating the Koszul formula for the Levi-Civita connection gives a tensorial first-variation formula for the connection; after expressing that formula in the chart $(U,x)$, the only first-order terms in $h$ are the three covariant derivatives displayed below. The terms produced by differentiating $(g_s^{-1})_{k\ell}$ contain no derivatives of $h$, and after the principal symbol freezes coefficients at $p$, all such coefficient-variation terms are lower order.
For a smooth one-parameter family of metrics $g_s$ with $g_0=g$ and $\frac{d}{ds}\big|_{s=0}g_s=h$, the first variation of the Christoffel coefficients is
\begin{align*}
(D\Gamma_{ij,k})_g(h)
=
\frac{1}{2}\sum_{\ell=1}^n (g^{-1})_{k\ell}
\left(
\nabla_i h_{j\ell}
+
\nabla_j h_{i\ell}
-
\nabla_\ell h_{ij}
\right)
+
\text{terms of order zero in } h.
\end{align*}
Here $\nabla_i h_{j\ell}$ denotes $(\nabla_{\partial_{x_i}}h)(\partial_{x_j},\partial_{x_\ell})$.
By tracing the coordinate expression for the Riemann curvature tensor in the definition of the Ricci tensor, the Ricci tensor has local coordinate expression
\begin{align*}
\operatorname{Ric}(g)_{ij}
=
\sum_{k=1}^n \partial_{x_k}\Gamma_{ij,k}
-
\sum_{k=1}^n \partial_{x_j}\Gamma_{ik,k}
+
\text{terms quadratic in the Christoffel coefficients}.
\end{align*}
When this expression is linearized, the quadratic Christoffel terms contribute only zeroth-order and first-order terms in $h$. Replacing the exterior partial derivatives by covariant derivatives changes only lower-order terms, because the difference between $\partial_{x_i}$ and $\nabla_i$ on tensor components is multiplication by Christoffel coefficients. Define the two second-order tensor expressions $A(h)$ and $B(h)$ by
\begin{align*}
A(h)_{ij}
&:=
\frac{1}{2}\sum_{k,\ell=1}^n (g^{-1})_{k\ell}
\nabla_k
\left(
\nabla_i h_{j\ell}
+
\nabla_j h_{i\ell}
-
\nabla_\ell h_{ij}
\right),
\end{align*}
and
\begin{align*}
B(h)_{ij}
&:=
\frac{1}{2}\sum_{k,\ell=1}^n (g^{-1})_{k\ell}
\nabla_j
\left(
\nabla_i h_{k\ell}
+
\nabla_k h_{i\ell}
-
\nabla_\ell h_{ik}
\right).
\end{align*}
Hence the second-order part of $D\operatorname{Ric}_g(h)$ is
\begin{align*}
(D\operatorname{Ric})_g(h)_{ij}
&\equiv A(h)_{ij}-B(h)_{ij},
\end{align*}
where $\equiv$ denotes equality modulo terms of differential order at most one in $h$.
[guided]
We only need the second-order part because the principal symbol ignores every term involving at most one derivative of $h$. The Christoffel symbols depend on first derivatives of the metric, so their variation contains first derivatives of $h$. Differentiating the varied Christoffel symbols inside the Ricci tensor is therefore the unique source of second derivatives of $h$.
Fix $p \in M$ and a smooth coordinate chart $(U,x)$ around $p$. The notation $\Gamma_{ij,k}$ means the scalar function on $U$ defined by
\begin{align*}
\nabla_{\partial_{x_i}}\partial_{x_j}=\sum_{k=1}^n \Gamma_{ij,k}\partial_{x_k}.
\end{align*}
Let $g_s$ be a smooth family of metrics with $g_0=g$ and variation $h=\frac{d}{ds}\big|_{s=0}g_s$. Differentiating the Koszul formula for the Levi-Civita connection gives the tensorial variation formula for $(D\nabla)_g(h)$; in the chart $(U,x)$, its first-order part has coefficients
\begin{align*}
(D\Gamma_{ij,k})_g(h)
=
\frac{1}{2}\sum_{\ell=1}^n (g^{-1})_{k\ell}
\left(
\nabla_i h_{j\ell}
+
\nabla_j h_{i\ell}
-
\nabla_\ell h_{ij}
\right)
+
\text{terms of order zero in } h.
\end{align*}
The displayed part is first order in $h$. The omitted part comes from differentiating inverse-metric and coefficient factors, contains no derivative of $h$, and after freezing coefficients at $p$ contributes only lower-order terms to the principal symbol.
Tracing the coordinate expression for the curvature tensor gives the Ricci tensor formula
\begin{align*}
\operatorname{Ric}(g)_{ij}
=
\sum_{k=1}^n \partial_{x_k}\Gamma_{ij,k}
-
\sum_{k=1}^n \partial_{x_j}\Gamma_{ik,k}
+
\text{terms quadratic in the Christoffel coefficients}.
\end{align*}
The quadratic Christoffel terms contain only the metric and its first derivatives. Their linearization therefore contains at most first derivatives of $h$, so they are lower order for a second-order operator. Also, replacing $\partial_{x_i}$ by $\nabla_i$ changes tensor-component formulas by terms involving Christoffel coefficients multiplying lower derivatives of $h$; this again changes only lower-order terms. Thus the second-order part is obtained by differentiating only the first-order part of $D\Gamma$:
\begin{align*}
(D\operatorname{Ric})_g(h)_{ij}
&\equiv
\sum_{k=1}^n \nabla_k (D\Gamma_{ij,k})_g(h)
-
\sum_{k=1}^n \nabla_j (D\Gamma_{ik,k})_g(h).
\end{align*}
Substituting the first-variation formula for $(D\Gamma)_g(h)$ into this expression gives
\begin{align*}
(D\operatorname{Ric})_g(h)_{ij}
&\equiv
\frac{1}{2}\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\nabla_k\left(\nabla_i h_{j\ell}+\nabla_j h_{i\ell}-\nabla_\ell h_{ij}\right)
-
\frac{1}{2}\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\nabla_j\left(\nabla_i h_{k\ell}+\nabla_k h_{i\ell}-\nabla_\ell h_{ik}\right).
\end{align*}
The symbol $\equiv$ means equality after discarding terms with at most one derivative of $h$. This is exactly the [equivalence relation](/page/Equivalence%20Relation) relevant to the principal symbol.
[/guided]
[/step]
[step:Replace the second covariant derivatives by covector components]
The operator $L_g=D(-2\operatorname{Ric})_g$ is $-2$ times the preceding linearization. The second-order part of $L_g(h)$ is therefore
\begin{align*}
L_g h_{ij}
&\equiv
-\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\nabla_k\left(\nabla_i h_{j\ell}+\nabla_j h_{i\ell}-\nabla_\ell h_{ij}\right)
+
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\nabla_j\left(\nabla_i h_{k\ell}+\nabla_k h_{i\ell}-\nabla_\ell h_{ik}\right).
\end{align*}
The principal symbol is obtained by replacing each occurrence of $\nabla_a\nabla_b h$ by multiplication by $\xi_a\xi_b h$ and freezing the coefficients at $p$. Hence
\begin{align*}
\sigma_\xi(L_g)(h)_{ij}
&=
-\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_k\left(\xi_i h_{j\ell}+\xi_j h_{i\ell}-\xi_\ell h_{ij}\right)
+
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_j\left(\xi_i h_{k\ell}+\xi_k h_{i\ell}-\xi_\ell h_{ik}\right).
\end{align*}
Using the symmetry $h_{ab}=h_{ba}$ and the definitions of $|\xi|_g^2$ and $\operatorname{tr}_g h$, this becomes
\begin{align*}
\sigma_\xi(L_g)(h)_{ij}=|\xi|_g^2 h_{ij}+\xi_i\xi_j\operatorname{tr}_g h-\sum_{k,\ell=1}^n \xi_i(g^{-1})_{k\ell}\xi_\ell h_{kj}-\sum_{k,\ell=1}^n \xi_j(g^{-1})_{k\ell}\xi_\ell h_{ki}.
\end{align*}
This is the asserted formula.
[guided]
The operator under consideration is $L_g=D(-2\operatorname{Ric})_g$, so its second-order part is obtained by multiplying the second-order part of $D\operatorname{Ric}_g$ by $-2$. Thus
\begin{align*}
L_g h_{ij}
&\equiv
-\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\nabla_k\left(\nabla_i h_{j\ell}+\nabla_j h_{i\ell}-\nabla_\ell h_{ij}\right)
+
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\nabla_j\left(\nabla_i h_{k\ell}+\nabla_k h_{i\ell}-\nabla_\ell h_{ik}\right).
\end{align*}
To form the principal symbol at $p$ in the covector direction $\xi\in T_p^*M$, we freeze the coefficient functions $(g^{-1})_{k\ell}$ at $p$ and replace each second covariant derivative $\nabla_a\nabla_b h$ by multiplication by $\xi_a\xi_b h$. Lower-order terms do not contribute by definition of the principal symbol. This gives
\begin{align*}
\sigma_\xi(L_g)(h)_{ij}
&=
-\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_k\left(\xi_i h_{j\ell}+\xi_j h_{i\ell}-\xi_\ell h_{ij}\right)
+
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_j\left(\xi_i h_{k\ell}+\xi_k h_{i\ell}-\xi_\ell h_{ik}\right).
\end{align*}
Now use the symmetry $h_{ab}=h_{ba}$. The positive term in the first sum is
\begin{align*}
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_k\xi_\ell h_{ij}=|\xi|_g^2h_{ij},
\end{align*}
and the trace term in the second sum is
\begin{align*}
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_j\xi_i h_{k\ell}=\xi_i\xi_j\operatorname{tr}_g h.
\end{align*}
Relabelling dummy indices in the two remaining contraction terms yields
\begin{align*}
\sigma_\xi(L_g)(h)_{ij}=|\xi|_g^2 h_{ij}+\xi_i\xi_j\operatorname{tr}_g h-\sum_{k,\ell=1}^n \xi_i(g^{-1})_{k\ell}\xi_\ell h_{kj}-\sum_{k,\ell=1}^n \xi_j(g^{-1})_{k\ell}\xi_\ell h_{ki}.
\end{align*}
This is precisely the displayed formula in the theorem statement.
[/guided]
[/step]
[step:Exhibit kernel directions generated by symmetric products with $\xi$]
Assume now that $\xi \ne 0$. For any $\eta \in T_p^*M$, define $h_{\xi,\eta}\in S^2T_p^*M$ by
\begin{align*}
(h_{\xi,\eta})_{ij}:=\xi_i\eta_j+\xi_j\eta_i.
\end{align*}
Let
\begin{align*}
a:=\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_k\eta_\ell.
\end{align*}
Then
\begin{align*}
\operatorname{tr}_g h_{\xi,\eta}
=
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}(\xi_k\eta_\ell+\xi_\ell\eta_k)
=
2a,
\end{align*}
and
\begin{align*}
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_\ell (h_{\xi,\eta})_{kj}=\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_\ell(\xi_k\eta_j+\xi_j\eta_k)=|\xi|_g^2\eta_j+a\xi_j.
\end{align*}
The analogous contraction with the index $i$ gives
\begin{align*}
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_\ell (h_{\xi,\eta})_{ki}
=
|\xi|_g^2\eta_i+a\xi_i.
\end{align*}
Substituting these three identities into the symbol formula gives
\begin{align*}
\sigma_\xi(L_g)(h_{\xi,\eta})_{ij}=|\xi|_g^2(\xi_i\eta_j+\xi_j\eta_i)+2a\xi_i\xi_j-\xi_i(|\xi|_g^2\eta_j+a\xi_j)-\xi_j(|\xi|_g^2\eta_i+a\xi_i)=0.
\end{align*}
Thus $h_{\xi,\eta}\in \ker\sigma_\xi(L_g)$ for every $\eta\in T_p^*M$.
Since $\xi\ne 0$, choosing $\eta=\xi$ gives
\begin{align*}
h_{\xi,\xi}=2\xi\otimes\xi\ne 0.
\end{align*}
Therefore $\ker\sigma_\xi(L_g)$ contains a non-zero element, so $\sigma_\xi(L_g)$ has non-trivial kernel for every non-zero $\xi\in T_p^*M$.
[guided]
Fix a non-zero covector $\xi\in T_p^*M$. The kernel directions we test are the symmetric products of $\xi$ with another covector. For any $\eta\in T_p^*M$, define the symmetric tensor $h_{\xi,\eta}\in S^2T_p^*M$ by
\begin{align*}
(h_{\xi,\eta})_{ij}:=\xi_i\eta_j+\xi_j\eta_i.
\end{align*}
Define the scalar $a\in\mathbb{R}$ by
\begin{align*}
a:=\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_k\eta_\ell.
\end{align*}
This scalar is the $g$-[inner product](/page/Inner%20Product) of the covectors $\xi$ and $\eta$ after raising one index with $g^{-1}$.
First compute the trace of $h_{\xi,\eta}$. Since $g^{-1}$ is symmetric and $h_{\xi,\eta}$ is symmetric,
\begin{align*}
\operatorname{tr}_g h_{\xi,\eta}
&=
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}(\xi_k\eta_\ell+\xi_\ell\eta_k)
=
2a.
\end{align*}
Next compute the contraction appearing in the third term of the symbol. Expanding the definition of $h_{\xi,\eta}$ gives
\begin{align*}
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_\ell (h_{\xi,\eta})_{kj}
=
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_\ell(\xi_k\eta_j+\xi_j\eta_k).
\end{align*}
Separating the two summands and using the definition of $a$ gives
\begin{align*}
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_\ell (h_{\xi,\eta})_{kj}
=
|\xi|_g^2\eta_j+a\xi_j.
\end{align*}
The same computation with $j$ replaced by $i$ gives
\begin{align*}
\sum_{k,\ell=1}^n (g^{-1})_{k\ell}\xi_\ell (h_{\xi,\eta})_{ki}
&=
|\xi|_g^2\eta_i+a\xi_i.
\end{align*}
Substituting these three identities into the principal-symbol formula yields
\begin{align*}
\sigma_\xi(L_g)(h_{\xi,\eta})_{ij}
=
|\xi|_g^2(\xi_i\eta_j+\xi_j\eta_i)+2a\xi_i\xi_j
-\xi_i(|\xi|_g^2\eta_j+a\xi_j)
-\xi_j(|\xi|_g^2\eta_i+a\xi_i).
\end{align*}
The two $|\xi|_g^2$ terms cancel the corresponding products $|\xi|_g^2\xi_i\eta_j$ and $|\xi|_g^2\xi_j\eta_i$, while the two $a$-terms cancel $2a\xi_i\xi_j$. Hence
\begin{align*}
\sigma_\xi(L_g)(h_{\xi,\eta})_{ij}=0.
\end{align*}
Thus every tensor $h_{\xi,\eta}$ lies in $\ker\sigma_\xi(L_g)$. To make this kernel element non-zero, choose $\eta=\xi$. Since $\xi\ne 0$, the tensor
\begin{align*}
h_{\xi,\xi}=2\xi\otimes\xi
\end{align*}
is non-zero in $S^2T_p^*M$. Hence $\ker\sigma_\xi(L_g)$ contains a non-zero element for every non-zero $\xi\in T_p^*M$.
[/guided]
[/step]
