I am working with a function $F : \mathbb{R}^3 \to \mathbb{R}^3$ and need to compute the vector of quadratic forms $Q$ given by $$Q(\textbf{u}) = \frac{1}{2} \textbf{u}^T\frac{d^2 F}{d\textbf{u}^2}(x,y,z)\textbf{u},$$ where $(x,y,z)$ is a point in $\mathbb{R}^3$ that should be consider to be fixed in the computation of $Q$.
I vaguely remember learning about higher order total derivatives in my graduate level introductory analysis course, but cannot seem to find my notes. I guess what would be best is if someone could point me to a nice reference.
Suppose $x \in \mathbb{R}^n$ and $A \in \mathbb{R}^{n \times n}$. We can then write $$x^{\rm T}Ax =\sum_{i, \,j=1}^{n}A_{ij}x_ix_j$$
Now suppose $B \in \mathbb{R}^{n \times n \times n}$ (this means that $B$ is a 3D-matrix). Then we can compute either a number
$$\sum_{i, \,j, \, k=1}^{n}B_{ijk}x_ix_jx_k$$
or something else. For example, we can get a vector $y$ with coordinates
$$ y_i=\sum_{j, \,k=1}^{n}B_{ijk}x_jx_k, \quad {\rm or} \quad y_j=\sum_{i, \,k = 1}^{n}B_{ijk}x_ix_k, \quad {\rm or} \quad y_k=\sum_{i, \,j =1}^{n}B_{ijk}x_ix_j $$
or maybe in some other way.
Your "total second derivative" is in fact a 3D-matrix (each 2D-layer is a derivative of the Jacobian by one of the variables) and so you should do something similar to find the coordinates of the resulting vector.