Given $d>1$, $\mathbf{x}\in{\rm I\!R}^d$, and two multivariate functions $\mathbf{u},\mathbf{v} : {\rm I\!R}^d\rightarrow{\rm I\!R}^d$, what is the second derivative of $f : {\rm I\!R}^d\rightarrow{\rm I\!R} $ $$ f(\mathbf{x}) = \mathbf{u}(\mathbf{x})\cdot\mathbf{v}(\mathbf{x}) = \mathbf{u}(\mathbf{x})^T\mathbf{v}(\mathbf{x}) $$ with respect to $\mathbf{x}$?
The first derivative is $$ \frac{\partial f}{\partial\mathbf{x}}(\mathbf{x}) = \mathbf{u}(\mathbf{x})^T\frac{\partial\mathbf{v}}{\partial\mathbf{x}}(\mathbf{x}) + \mathbf{v}(\mathbf{x})^T\frac{\partial\mathbf{u}}{\partial\mathbf{x}}(\mathbf{x}) $$ which is a row vector of size $d$ when using the numerator layout.
But I do not know how to differentiate again since $\frac{\partial\mathbf{u}}{\partial\mathbf{x}}$ and $\frac{\partial\mathbf{v}}{\partial\mathbf{x}}$ are $3$-by-$3$ matrices...
I tried to develop these equations in a simple situation where $d=2$, but I obtain a quite complicated matrix and I am looking for a clear formula, in any dimension $d>1$, where only $\mathbf{u},\mathbf{v}$ and their derivatives appear.
It's easier to work in index notation. Write $f=\sum_{i=1}^du_iv_i$ so$$(\nabla f)_j=\sum_i[(\partial_ju_i)v_i+u_i\partial_jv_i]=[(\nabla u)v+(\nabla v)u]_j,$$i.e. $\nabla f=(\nabla u)v+(\nabla v)u$. Now for the second derivative:$$(\nabla^2f)_{jk}=\partial_k(\nabla f)_j=\sum_i[(\partial_j\partial_ku_i)v_i+\partial_ju_i\partial_kv_i+\partial_ku_i\partial_jv_i+u_i\partial_j\partial_kv_i].$$Note this includes rank-$3$ tensors we might denote $\nabla^2u,\,\nabla^2v$, but writing in terms of those without indices requires some clarification on the notation. For example, if we take $(\nabla^2u)_{jki}:=(\partial_j\partial_ku_i)$, we get$$\nabla^2f=(\nabla^2u)\cdot v+(\nabla u)\cdot(\nabla v)+(\nabla v)\cdot(\nabla u)+(\nabla^2v)\cdot u,$$where each $\cdot$ contracts the factors' last index. This is the general Leibniz rule with non-commuting factors (see here for a treatment of the commuting case, e.g. if $u,\,v$ were scalars).