First of all, this question may be more suited for the Math stack exchange site.  If anyone finds this question irrelevant here, please transfer to the relevant site.

Recall that in terms of Weyl and Ricci tensor, the Bach tensor in $4$-dimension can be expressed as the following:

\begin{equation}
B_{\alpha \beta} \equiv \nabla^\mu \nabla^\nu C_{\mu \alpha \beta \nu}-\frac{1}{2} C_{\mu \alpha \beta \nu} R^{\mu \nu}.
\end{equation}

Problem Setup:

Let $(M^4, g)$ be a space-time with local coordinates $x^\alpha(\alpha=0,1,2,3)$ and $\Sigma^3 \subset M^4$ a spacelike hypersurface with intrinsic coordinates $t^i(i=1,2,3)$, unit normal $n=n^\alpha \partial_\alpha$, and connecting quantities $x_i^\alpha=\partial x^\alpha / \partial t^i$. The first and second fundamental forms of $\Sigma$

$$
\begin{aligned}
\mathrm{I} & =g_{i j} \,d t^i\, dt^j = g_{\alpha \beta} x_i^\alpha \, dt^i \, x_j^\beta \, d t^j, \\
\mathrm{II} & =P_{i j} \, dt^i \, dt^j = \left(\nabla_\alpha n_\beta\right) x_i^\alpha \, d t^i \, x_j^\beta \, dt^i
\end{aligned}
$$

are defined above. Let's define two more fundamental forms on $\Sigma$:

$$
\begin{aligned}
& \mathrm{III}=Q_{i j} \, dt^i \, dt^j = n^\mu x_i^\alpha x_j^\beta n^\nu C_{\mu \alpha \beta \nu} \, dt^i \, dt^j \\
& \mathrm{IV}=S_{i j} \, dt^i \, dt^j = n^\mu x_i^\alpha x_j^\beta \, \nabla^\nu C_{\mu \alpha \nu \beta} \, d t^i \, d t^j
\end{aligned}
$$

We also use shorthand notations:

$$
\begin{gathered}
B_{n n}=n^\alpha n^\beta B_{\alpha \beta}, \quad B_{n i}=n^\alpha x_i^\beta B_{\alpha \beta}, \quad B_{i j} = x_i^\alpha x_j^\beta B_{\alpha \beta} \\
R_{i j k l}=x_i^\alpha x_j^\nu x_k^\mu x_l^{\beta} R_{\alpha \beta \mu \nu} \\
R_{i j k n}=x_i^\alpha x_j^\beta x_k^\mu n^\nu R_{\alpha \beta \mu \nu}
\end{gathered}
$$

etc. Gauss-Codazzi equations tell us that $R_{i j k l}, R_{i j k n}, C_{i j k l}, C_{i j k n}$ are $\Sigma$-intrinsic differential expressions in I, II. This can be established using van der Waerden's $D_i$-symbols: The $D_i$ act as $x_i^\alpha \nabla_\alpha$ on $M$-tensors and as $i$-covariant derivatives on $\Sigma$-tensors.

Problem/Question Statement:

In an old paper titled "Cauchy's problem for Bach's equations of General Relativity" 
R. Schimming stated the following formula for Bach's normal-normal component without proof:

$$
B_{n n}=D^i\left(D^j Q_{i j}+P^{j k} C_{i j k n}\right)-Q^{i j}\left(Q_{i j}+R^{k}_{i k j}\right)+P^{i j} S_{i j}.
$$

Can anyone please recommend any source from which I might find its derivation? I spent a few days trying to derive this, but my calculations are getting nastier, particularly $ Q_{i j} Q^{i j}$ becomes very wild. Again, I apologize if this question is silly. I would be grateful if you could suggest a nice (or any) way of doing this calculation. Thank you very much.

Edit: I corrected the formula after Professor Case's comments.