Let $(M, g)$ be a Riemannian manifold and $\Psi: M \rightarrow M$ a diffeomorphism. To simplify things we can take $g$ to be the Euclidean metric.
From this question it seems that if $\Psi$ is an isometry, then for any function $f: M \rightarrow \mathbb{R}$ we can write the covariant Hessian as $$ \left[\nabla d\right]f(\Psi({\bf x})) = \left[d\Psi\right]^\top\big \vert_{\bf x} \left[\nabla d f\right]\vert_{\Psi({\bf x})} \left[d\Psi\right]\big \vert_{\bf x},$$ where $\left[d\Psi\right]^\top$ is the adjoint of $d\Psi$.
Different from the prior question, I want to understand how the covariant Hessian transforms when $\Psi$ is not an isometry. Specifically, I'm most interested in the case when $$\textrm{det}\left(\left[d\Psi\right]\vert_{\bf x}\right) = 1, \ \forall {\bf x} \in M.$$
Clearly, applying the chain rule once in terms of the exterior derivative gives $$ \left[\nabla d\right]f(\Psi({\bf x})) = \nabla \left( \left[d\Psi\right]^\top\big \vert_{\bf x} \left[d f\right]\vert_{\Psi({\bf x})} \right),$$ though this is the point where I get confused. Specifically, I'm not sure how the coviariant derivative $\nabla$ acts on $ \left[d\Psi\right]^\top\big \vert_{\bf x}$.
Naively plugging away at this by applying the chain rule again, it seems like we should get something like $$\nabla \left( \left[d\Psi\right]^\top\big \vert_{\bf x} \left[d f\right]\vert_{\Psi({\bf x})} \right) = \left[\nabla [d \Psi]^\top\right]\big \vert_{\bf x} \left[d f\right]\vert_{\Psi({\bf x})} \ + \left[d\Psi\right]^\top\big \vert_{\bf x} \left[\nabla d f\right]\vert_{\Psi({\bf x})} \left[d\Psi\right]\big \vert_{\bf x}.$$
In the case where $\Psi$ is an isometry, it appears we have $\nabla [d \Psi]^\top = 0,$ however, when $\Psi$ is not an isometry, I don't understand how $\nabla$ to acts on $[d \Psi]^\top.$ Furthermore, since the covariant Hessian is symmetric by definition, how is $$\left[\nabla [d \Psi]^\top\right]\big \vert_{\bf x} \left[d f\right]\vert_{\Psi({\bf x})}$$ symmetric?
I'll consider a slightly more general case than yours: Let $\Psi:M\to N$ be a smooth map, where $(M,g)$ and $(N,\widetilde{g})$ are riemaniian manifolds. Throughout I'll use $x^a$ to denote local coordinates in $M$, $y^\alpha$ to denote local coordinates in $N$, and tildes ($\widetilde\nabla$, etc.) to denote objects in $N$. We can push forward tangent vectors with the differential $(d\Psi)_p :T_pM\to T_{\Psi(p)}N$, and pull back covariant tensors fields from $N$ to $M$, by the formula $(\Psi^*T)(X_1,\dots,X_k)=T(d\Psi X_1,\dots,d\Psi X_k)$. In particular, the pullback of a smooth function $f:N\to\mathbb{R}$ is $\Psi^*f=f\circ\Psi$ and the pullback of a covector $\omega$ is $\Psi^*\omega=\omega\circ d\Psi$.
First, we can describe how the connections differ under $\Psi$. To do this, given a covector field $\omega$ on $N$, we can either take the covariant derivative in $N$ and pull the result back to $M$, or pull $\omega$ back to $M$ and then take the covariant derivative there. We define $A(\omega)$ as the difference of these two operations. $$ A(\omega)=\nabla(\Psi^*\omega)-\Psi^*(\widetilde\nabla\omega) $$ One can show that $A(\omega)$ depends tensorially on $\omega$ (i.e. $A(\omega)|_p$ is a linear function of $\omega|_{f(p)}$). This allows us to give a well-defined coordinate representation of $A$. Given $\omega=\omega_\gamma dy^\gamma$, we have $A(\omega)=A^{\gamma}{}_{ab}\omega_\gamma dx^a\otimes dx^b$. With a bit of computation, one can show $$ A^{\gamma}{}_{ab}=\frac{\partial y^\gamma}{\partial x^a\partial x^b}+\frac{\partial y^\alpha}{\partial x^a}\frac{\partial y^\beta}{\partial x^b}\widetilde\Gamma^{\gamma}{}_{\alpha\beta}-\frac{\partial y^\gamma}{\partial x^c}\Gamma^{c}{}_{ab} $$ We can now do the same thing for the second covariant derivative: given a smooth function $f:N\to\mathbb{R}$, we can compute the Hessian of $f$ in $N$ and pull it back to $M$, or pull $f$ back to $M$ then compute the Hessian. We can use the fact that the differential commutes with pullbacks $d(\Psi^* f)=\Psi^*(df)$ to write the difference between these two in terms of $A$ $$\begin{align} \operatorname{Hess}(\Psi^*f)-\Psi^*\widetilde{\operatorname{Hess}}(f)&=\nabla d(\Psi^*f)-\Psi^*(\widetilde\nabla df) \\ &=\nabla(\Psi^*df)-\Psi^*(\widetilde\nabla df) \\ &=A(df) \end{align}$$ Or, in coordinates, $$ \nabla_a\nabla_b(f\circ\Psi)-\frac{\partial y^\alpha}{\partial x^a}\frac{\partial y^\beta}{\partial x^b}\widetilde{\nabla}_\alpha\widetilde{\nabla}_\beta f=A^\gamma{}_{ab}\frac{\partial f}{\partial y^\gamma} $$ We see that $A$ plays the role of $\nabla d\Psi$ in your computation, and indeed it is a kind of "second derivative" of $\Psi$. However, we typically only define covariant derivatives of tensors on a single manifold, so $\nabla d\Psi$ doesn't make sense, since $d\Psi$ is a map between different tangent spaces. One could attempt to find a generalization of the covariant derivative that allows one to interpret things like $\nabla d\Psi$, bu since all we need here is $A$, that's probably more trouble than it's worth.