$\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 its flow.
Is there a reasonable way to express then $n$-th derivative $\frac{\partial^n }{\partial t^n}\left|_{t=0}\right.\tr_g\big(\psi_t^*g\big)=\operatorname{tr}_g(\mathcal{L}_X \dots \mathcal{L}_X(\mathcal{L}_Xg))$ for $n \ge 3$?
I am particularly interested in the case where $X$ is divergence-free.
For $n=1$, we have $\frac{\partial }{\partial t}\left|_{t=0}\right.\tr_g\big(\psi_t^*g\big)=\tr_g(L_Xg)=2\div (X)$.
For $n=2$, this answer here gives $$\frac{\partial^2 }{\partial t^2}\left|_{t=0}\right.\tr_g\big(\psi_t^*g\big)=\operatorname{tr}_g(\mathcal{L}_X(\mathcal{L}_Xg))=2X(\operatorname{div}X)+2\|\nabla X\|^2+2\operatorname{tr}(\nabla X\circ\nabla X)$$
I am not even sure how to approach the case of $n=3$.
Here's a (hopefully) more straightforward description of the computation from the previous answer, and how it applies to higher derivatives.
Let $T$ be a $(0,2)$ tensor field. The Lie derivative distributes over tensor evaluations, so for $U,V\in\mathfrak{X}M$, we have $$ (\mathcal{L}_XT)(U,V)=\mathcal{L}_X(T(U,V))-T(\mathcal{L}_XU,V)-T(U,\mathcal{L}_XV) $$ The same is true of the Levy-Civita connection, which is to say we may replace $\mathcal{L}_X$ with $\nabla_X$ in the above expression. Combining these two identities with the torsion free condition $\mathcal{L}_UV=\nabla_UV-\nabla_VU$, we obtain $$ (\mathcal{L}_XT)(U,V)=(\nabla_XT)(U,V)+T(\nabla_UX,V)+T(U,\nabla_VX) $$ or, in abstract index notation $$ (\mathcal{L}_XT)_{ab}=X^cT_{ab;c}+X^c_{;a}T_{cb}+T_{ac}X^c_{;b} $$ To compute the higher order Lie derivatives of the metric in terms of $X$ and its covariant derivatives, one can apply this formula repeatedly to $g$. Since the results become quite complicated with successive applications, there is some choice in how to organize this computation.
One approach involves writing everything in terms of $(1,1)$ tensors (fiberwise linear maps $TM\to TM$) and using matrix notation. To this end, define $$ N^a{}_b=X^a_{;b}\ \ \ \ \ \ \ \bar{N}^a{}_b=g^{ac}g_{cd}X^d_{;c}\ \ \ \ \ \ \ (L^{(n)})^a{}_b=g^{ac}(\mathcal{L}_X^ng)_{cb} $$ That is, $N$ is the (abstract) matrix representation of $\nabla X$, $\bar{N}$ is its adjoint w.r.t. the metric, and $L^{(n)}$ are the matrices whose traces we would like to compute. In terms of these matrices, the recurrence relation becomes $$ L^{(0)}=I,\ \ \ \ \ \ \ L^{(n+1)}=\nabla_X(L^{(n)})+\bar{N}L^{(n)}+L^{(n)}N $$ Computing the first few cases, we have $$ L^{(1)}=\bar{N}+N \\ L^{(2)}=\nabla_X(\bar{N}+N)+\bar{N}^2+2\bar{N}N+N^2 \\ L^{(3)}=\nabla_X^2(\bar{N}+N)+\nabla_X(\bar{N}^2+2\bar{N}N+N^2) \\ +\bar{N}\nabla_X(\bar{N}+N)+\nabla_X(\bar{N}+N)N \\ +\bar{N}^3+3\bar{N}^2N+3\bar{N}N^2+N^3 $$ and so on. Simplifying the traces of these matrices is still not so straightforward, and the combinatorics to describe the general term are rather involved. Nonetheless, they can be dealt with individually using the fact that the trace is invariant under adjoints and cyclic permutations and commutes with $\nabla_X$. For instance, the first three cases simplify to $$ \operatorname{tr}_g(\mathcal{L}_Xg)=2\operatorname{tr} N \\ \operatorname{tr}_g(\mathcal{L}_X^2g)=2\nabla_X(\operatorname{tr}N)+2\operatorname{tr}(\bar{N}N)+2\operatorname{tr}(N^2) \\ \operatorname{tr}_g(\mathcal{L}_X^3g)=2\nabla_X^2(\operatorname{tr}N)+3\nabla_X(\operatorname{tr}(N^2)+\operatorname{tr}(\bar{N}N))+2\operatorname{tr}(N^3)+6\operatorname{tr}(\bar{N}N^2) $$ Subsequent entries can likewise be written as traces of products of covariant derivatives of $\bar{N}$ and $N$.