$\newcommand{\tr}{\operatorname{tr}}$ $\renewcommand{\div}{\operatorname{div}}$
Let $(M,g)$ be a smooth Riemannian manifold. Given a vector field $X$ on $M$, Let $\psi_t:M \to M$ be the flow of $X$.
Is there a reasonable way to express $\frac{\partial^2 }{\partial t^2}\left|_{t=0}\right.\tr_g\big((\psi_t^*g)_p\big)$?
If we denote $f(t,p)=\tr_g\big((\psi_t^*g)_p\big)$ then $\frac{\partial }{\partial t}f(t,p)=\tr_g(\psi_t^*L_Xg)$. (In particular, $\frac{\partial }{\partial t}f(0,p)=\tr_g(L_Xg)=2\div (X)$.)
So, the question is how to obtain an expression for $\frac{\partial }{\partial t} \tr_g(\psi_t^*L_Xg)$.
The expression you wtite down is the trace of the second Lie derivative of the metric $\operatorname{tr}_g(\mathcal{L}_X(\mathcal{L}_Xg))$. This can be rewritten a bit using properties of the Lie derivative. Let $\hat{g}$ be the inverse metric. The trace $\operatorname{tr}_g$ can be seen as a tensor product with $\hat{g}$ followed by a contraction $\delta:T^2_2 M\to\mathbb{R}$. We can use this to expand the expression, since $\mathcal{L}_X$ commutes with contractions and distributes over tensor products. $$ \operatorname{Tr}_g(\mathcal{L}_X(\mathcal{L}_Xg))=\delta[\hat{g}\otimes\mathcal{L}_X(\mathcal{L}_X(g))] \\ =\delta[\mathcal{L}_X(\mathcal{L}_X(\hat{g}\otimes g))-\mathcal{L}_X(\mathcal{L}_X\hat{g}\otimes g)-\mathcal{L}_X\hat{g}\otimes\mathcal{L}_Xg] \\ =\mathcal{L}_X(\mathcal{L}_X(\delta(\hat{g}\otimes g)))-\mathcal{L}_X(\delta(\mathcal{L}_X\hat{g}\otimes g))-\delta(\mathcal{L}_X\hat{g}\otimes\mathcal{L}_Xg) \\ $$ The first term vanishes, since it is the lie derivative of the constant function $\operatorname{tr}_gg$. The other two can be written as $$ =-\mathcal{L}_X((\mathcal{L}_X\hat{g})^{ij}g_{ij})-(\mathcal{L}_X\hat{g})^{ij}(\mathcal{L}_Xg)_{ij} $$ Then using the fact that $(\mathcal{L}_X\hat{g})^{ij}=-\hat{g}^{ik}\hat{g}^{jl}(\mathcal{L}_Xg)_{kl}$ (which can be seen by differentiating the expression $\hat{g}^{ij}g_{jk}=\delta^i_k$), we get $$ =\mathcal{L}_X(\hat{g}^{ij}(\mathcal{L}_Xg)_{ij})+\hat{g}^{ik}\hat{g}^{jl}(\mathcal{L}_Xg)_{ij}(\mathcal{L}_Xg)_{kl} $$ Finally, using Killing's equation $(\mathcal{L}_X g)_{ij}=X_{i;j}+X_{j;i}$ $$ =2X^jX^i_{;ij}+2g_{ij}\hat{g}^{kl}X^i_{;k}X^j_{;l}+2X^i_{;j}X^j_{;i} $$ Or, in more coordinate free notation: $$ =2X(\operatorname{div}X)+2\|\nabla X\|^2+2\operatorname{tr}(\nabla X\circ\nabla X) $$ Where $\|\ \ \|$ is the norm on $(1,1)$ tensors induced by $g$ and $\circ$ is composition of $(1,1)$ tensors, thought of as fiberwise linear maps $TM\to TM$.