For ease of typing, let's use $\,x={\pmb\delta}\,$ as the independent variable.
We'll also use the fact that the Kronecker product of two vectors can be expanded in two ways:
$$a\otimes b = (I_a\otimes b)\,a = (a\otimes I_b)\,b$$ where $I_a$ is the identity matrix whose dimensions are compatible with the $a$ vector, while $I_b$ is compatible with the $b$ vector.
The differential and Jacobian of first function can be calculated as
$$\eqalign{
y &= x\otimes x \cr
dy &= dx\otimes x + x\otimes dx = \big(I\otimes x + x\otimes I\big)\,dx \cr
J = \frac{\partial y}{\partial x} &= I\otimes x + x\otimes I \cr
}$$
Using the previous result $\big(dy=J\,dx\big),\,$ your second function can be easily dispatched
$$\eqalign{
\phi &= x^TAy = y^TA^Tx \cr
d\phi
&= x^TA\,dy + y^TA^T\,dx \cr
&= \big(x^TAJ + y^TA^T\big)\,dx \cr
&= \big(Ay + J^TA^Tx\big)^T\,dx \cr
\frac{\partial\phi}{\partial x}
&= Ay + J^TA^Tx \cr
&= A(x\otimes x) + \big(I\otimes x^T + x^T\otimes I\big)A^Tx \cr
}$$
The first term is the same between you, Frank, and me.
My second term appears to be the transpose of Frank's and different from yours.
Short answer: the derivative of $Q\otimes Q\otimes Q\otimes Q$ with respect to $Q$ is a mess, at first sight...
Let's start simple. Let $Q$ be a $K\times K$ matrix with entries $Q_{ij}$ and let $E^{ab}$ be the $K\times K$ matrix with all $0$ entries, except the entry $(a,b)$ which is $1$; in other words, $(E^{ab})_{ij} = \delta_a^i\delta_b^j$.
Then I claim that
$$
\frac{\partial(Q\otimes Q)}{\partial Q_{ij}} = E^{ij}\otimes Q+Q\otimes E^{ij} .
$$
I leave it to you to see why, because trying to write out the involved matrices will probably crash the entire Stack Exchange network...
Jokes aside, this is really immediate to see: just write $Q\times Q$ as in the first formula of the definition and think which elements are affected by $Q_{ij}$. There is the entire $(i,j)$th block, so you get $E^{ij}\otimes Q$, but there is also the $(i,j)$th entry in each block, which gives you $Q\otimes E^{ij}$.
Now, if $A$ and $B$ are matrices which are functions of $Q$, by the same reasoning you get
$$
\frac{\partial(A\otimes B)}{\partial Q_{ij}} = \frac{\partial A}{\partial Q_{ij}}\otimes B + A\otimes \frac{\partial B}{\partial Q_{ij}} .
$$
So you can iterate for instance
$$
\begin{split}
\frac{\partial (Q\otimes Q\otimes Q)}{\partial Q_{ij}}
&= \frac{\partial(Q\otimes Q)}{\partial Q_{ij}}\otimes Q
+ (Q\otimes Q)\otimes \frac{\partial Q}{\partial Q_{ij}} \\
&= (E^{ij}\otimes Q+Q\otimes E^{ij})\otimes Q + (Q\otimes Q)\otimes E^{ij} \\
&= E^{ij}\otimes Q\otimes Q + Q\otimes E^{ij}\otimes Q + Q\otimes Q\otimes E^{ij}.
\end{split}
$$
Now you can prove by induction that
$$
\frac{\partial \bigl(\bigotimes_{n=1}^N Q\bigr)}{\partial Q_{ij}}
= \sum_{n=1}^N \left(\bigotimes_{h=1}^{n-1} Q\right) \otimes E^{ij} \otimes \left(\bigotimes_{h=n+1}^{N} Q\right).
$$
Written more concisely,
$$
\frac{\partial Q^{\otimes N}}{\partial Q_{ij}}
= \sum_{n=1}^N Q^{\otimes (n-1)}\otimes E^{ij} \otimes Q^{\otimes (N-n)} .
$$
Best Answer
For typing convenience, define the following variables $$x=\delta,\qquad Y=I_N\otimes x,\qquad A = Y^TMY,\qquad\lambda=\ln(\det(A))$$ Assume you are given an analytic function of a scalar variable and its derivative $$f=f(s),\qquad p(s)=\frac{df}{ds}$$ These functions can be applied to a square matrix argument to yield the matrix values $$F=f(A),\qquad P=p(A)$$ The gradient (or equivalently the differential) of the trace of the matrix function can be written as $$\frac{\partial {\rm Tr}(F)}{\partial A} = P^T \quad\iff\quad d\,{\rm Tr}(F) = P^T:dA$$ In this particular case, the scalar functions that will be needed are $$f(s)=\ln(s),\qquad p(s)=s^{-1}$$ and the matrices are symmetric.
Therefore $$\eqalign{ d\lambda &= d\,\ln(\det(A)) \\ &= d\,{\rm Tr}(\ln(A)) &\quad\big({\rm Jacobi\,formula}\big) \\ &= A^{-1}:dA \\ &= A^{-1}:(Y^TM\,dY+dY^TMY) \\ &= 2A^{-1}:Y^TM\,dY \\ &= 2MYA^{-1}:dY \\ &= 2MYA^{-1}:(I_N\otimes dx) &\quad\big({\rm Your\,formula}\big) \\ &= 2\;{\rm vec}\big(MYA^{-1}\big):\big({\rm vec}(I_N)\otimes I_p\big)\,dx\\ &= 2\,\big({\rm vec}(I_N)\otimes I_p\big)^T{\rm vec}\big(MYA^{-1}\big):dx\\ \frac{\partial\lambda}{\partial x} &= 2\,\big({\rm vec}(I_N)\otimes I_p\big)^T{\rm vec}\big(MYA^{-1}\big) \\ }$$
In several of the steps above, a colon is used to denote the trace/Frobenius product, i.e. $$\eqalign{ A:B = {\rm Tr}(A^TB) = {\rm Tr}(B^TA) = B:A }$$