Taking second derivative of multivariate normal density wrt covariance matrix In attempting to compute the second derivative of the density of a $d$-dimensional $\mathrm{MVN}(\pmb0,\Sigma)$ random variable with respect to $\Sigma$, I am running into an issue. In particular, I am having trouble figuring out the order of multiplication from element-wise notation. I am also having difficulties figuring out how to translate an expression from element-wise notation to matrix notation. I am looking for an answer to the order of multiplication and at least some tips on how to move forward with the final expression below.
For $\Sigma=(\sigma_{ij})$, set $\nabla_\Sigma=(\partial_{\sigma_{ij}})$. For the pdf, write
$$p(x)=(2\pi)^{-d/2}|\Sigma|^{-1/2}\exp\left(-\frac{1}{2}x^\top\Sigma^{-1}x\right).$$
Then we have
$$\nabla_\Sigma p(x)=(2\pi)^{-d/2}\left[|\Sigma|^{-1/2}\nabla_\Sigma\exp\left(-\frac{1}{2}x^\top\Sigma^{-1}x\right)\\
+\exp\left(-\frac{1}{2}x^\top\Sigma^{-1}x\right)\nabla_\Sigma|\Sigma|^{-1/2}\right]$$
$$=-\frac{1}{2}(2\pi)^{-1/2}\exp\left(-\frac{1}{2}x^\top\Sigma^{-1}x\right)\left[|\Sigma|^{-1/2}\Sigma^{-1}xx^\top\Sigma^{-1}-|\Sigma|^{-3/2}|\Sigma|\Sigma^{-1}\right]$$
$$=\frac{p(x)}{2}\left[\Sigma^{-1}-\Sigma^{-1}xx^\top\Sigma^{-1}\right].$$
Hence, $\nabla_\Sigma\nabla_\Sigma p(x)=[\nabla_\Sigma(p(x)\Sigma^{-1})-\nabla_\Sigma(p(x)\Sigma^{-1}xx^\top\Sigma^{-1})]/2$. Beginning with $\nabla_\Sigma(p(x)\Sigma^{-1})$, I go component-wise to get
$$\partial_{\sigma_{ij}}p(x)(\Sigma^{-1})_{kl}=p(x)\partial_{\sigma_{ij}}(\Sigma^{-1})_{kl}+(\Sigma^{-1})_{kl}\partial_{\sigma_{ij}}p(x)$$
$$=-p(x)(\Sigma^{-1})_{ik}(\Sigma^{-1})_{lj}+(\Sigma^{-1})_{kl}\partial_{\sigma_{ij}}p(x).$$
Here I come to my first question. Should $(\Sigma^{-1})_{kl}\partial_{\sigma_{ij}}p(x)$ correspond to $\Sigma^{-1}\otimes\nabla_\Sigma p(x)$ or $\nabla_\Sigma p(x)\otimes\Sigma^{-1}$? How can I tell?
Assuming the former provides
$$\nabla_\Sigma(p(x)\Sigma^{-1})=-p(x)\Sigma^{-1}\otimes\Sigma^{-1}+\Sigma^{-1}\otimes\nabla_\Sigma p(x)$$
$$=\frac{p(x)}{2}\Sigma^{-1}\otimes\left[\Sigma^{-1}-\Sigma^{-1}xx^\top\Sigma^{-1}\right].$$
For $\nabla_\Sigma[p(x)\Sigma^{-1}xx^\top\Sigma^{-1}]$, I also proceed component-wise:
$$\partial_{\sigma_{ij}}\left[p(x)\sum_{k,l}(\Sigma^{-1})_{il}x_lx_k(\Sigma^{-1})_{kj}\right]=p(x)\partial_{\sigma_{ij}}\sum_{k,l}(\Sigma^{-1})_{il}x_lx_k(\Sigma^{-1})_{kj}+\left[\partial_{\sigma_{ij}}p(x)\right]\sum_{k,l}(\Sigma^{-1})_{il}x_lx_k(\Sigma^{-1})_{kj}.$$
Since the right-most component corresponds to either $\nabla_\Sigma p(x)\otimes\Sigma^{-1}xx^\top\Sigma^{-1}$ or $\Sigma^{-1}xx^\top\Sigma^{-1}\otimes\nabla_\Sigma p(x)$, I focus on the left component, dropping the $p(x)$ for brevity
$$\partial_{\sigma_{ij}}\sum_{k,l}(\Sigma^{-1})_{il}x_lx_k(\Sigma^{-1})_{kj}=\sum_{k,l}(\Sigma^{-1})_{il}x_lx_k\partial_{\sigma_{ij}}(\Sigma^{-1})_{kj}+(\Sigma^{-1})_{kj}x_lx_k\partial_{\sigma_{ij}}(\Sigma^{-1})_{il}$$
$$=-\sum_{k,l}(\Sigma^{-1})_{il}x_lx_k(\Sigma^{-1})_{ki}(\Sigma^{-1})_{jj}+(\Sigma^{-1})_{kj}x_lx_k(\Sigma^{-1})_{ii}(\Sigma^{-1})_{jl}$$
$$=-(\Sigma^{-1})_{jj}\sum_{k,l}(\Sigma^{-1})_{il}x_lx_k(\Sigma^{-1})_{ki}-(\Sigma^{-1})_{ii}\sum_{k,l}(\Sigma^{-1})_{kj}x_lx_k(\Sigma^{-1})_{jl}.$$
Since I am not sure how to write this last expression in terms of matrices and Kronecker products, this is where my journey has ended so far. Any ideas for finishing this calculation would be greatly appreciated.
 A: Okay I think I figured it out:
$$\nabla_\Sigma p(z)=(2\pi)^{-d/2}\left[|\Sigma|^{-1/2}\nabla_\Sigma\exp\left(-\frac{1}{2}z^\top\Sigma^{-1}z\right)\\
+\exp\left(-\frac{1}{2}z^\top\Sigma^{-1}z\right)\nabla_\Sigma|\Sigma|^{-1/2}\right]$$
$$=\frac{1}{2}(2\pi)^{-1/2}\exp\left(-\frac{1}{2}z^\top\Sigma^{-1}z\right)\left(|\Sigma|^{-1/2}\Sigma^{-1}zz^\top\Sigma^{-1}-|\Sigma|^{-3/2}|\Sigma|\Sigma^{-1}\right)$$
$$=\frac{p(z)}{2}\left(\Sigma^{-1}zz^\top\Sigma^{-1}-\Sigma^{-1}\right).$$
Since $z$ is a column vector, $zz^\top$ is a $d\times d$ matrix which we denote $Z$. Hence, $\nabla_\Sigma\nabla_\Sigma p(z)=[\nabla_\Sigma(p(z)\Sigma^{-1}Z\Sigma^{-1})-\nabla_\Sigma(p(z)\Sigma^{-1})]/2$. Beginning with $\nabla_\Sigma(p(z)\Sigma^{-1})$, I go component-wise to get
$$\partial_{\sigma_{ij}}p(z)(\Sigma^{-1})_{kl}=p(z)\partial_{\sigma_{ij}}(\Sigma^{-1})_{kl}+(\Sigma^{-1})_{kl}\partial_{\sigma_{ij}}p(z)$$
$$=-p(z)(\Sigma^{-1})_{ik}(\Sigma^{-1})_{lj}+(\Sigma^{-1})_{kl}\partial_{\sigma_{ij}}p(z).$$
Choosing the convention to have letters earlier in the alphabet go last (so that $kl$ goes before $ij$ and $lj$ before $ik$) provides
$$\nabla_\Sigma(p(z)\Sigma^{-1})=-p(z)\Sigma^{-1}\otimes\Sigma^{-1}+\Sigma^{-1}\otimes\nabla_\Sigma p(z)$$
$$=-p(z)\Sigma^{-\otimes2}+\frac{p(z)}{2}\Sigma^{-1}\otimes(\Sigma^{-1}Z\Sigma^{-1}-\Sigma^{-1})$$
$$=\frac{p(z)}{2}\left[\Sigma^{-1}\otimes(\Sigma^{-1}Z\Sigma^{-1})-3\Sigma^{-\otimes2}\right],$$
Where $\Sigma\otimes\Sigma=\Sigma^{\otimes2}$. For $\nabla_\Sigma[p(z)\Sigma^{-1}Z\Sigma^{-1}]$, also proceed component-wise:
$$\partial_{\sigma_{ij}}\left[p(z)\sum_{m,n}(\Sigma^{-1})_{kn}z_nz_m(\Sigma^{-1})_{ml}\right]=p(z)\partial_{\sigma_{ij}}\sum_{m,n}(\Sigma^{-1})_{kn}z_nz_m(\Sigma^{-1})_{ml}+\left[\partial_{\sigma_{ij}}p(z)\right]\sum_{m,n}(\Sigma^{-1})_{kn}z_nz_m(\Sigma^{-1})_{ml}.$$
Using the established convention, the right-most component corresponds to $\Sigma^{-1}Z\Sigma^{-1}\otimes\nabla_\Sigma p(z)$. Focusing on the left component and dropping the $p(z)$ for brevity:
$$\partial_{\sigma_{ij}}\sum_{m,n}(\Sigma^{-1})_{kn}z_nz_m(\Sigma^{-1})_{ml}=\sum_{m,n}(\Sigma^{-1})_{kn}z_nz_m\partial_{\sigma_{ij}}(\Sigma^{-1})_{ml}+(\Sigma^{-1})_{ml}z_nz_m\partial_{\sigma_{ij}}(\Sigma^{-1})_{kn}$$
$$=-\sum_{m,n}(\Sigma^{-1})_{kn}z_nz_m(\Sigma^{-1})_{im}(\Sigma^{-1})_{lj}+(\Sigma^{-1})_{ml}z_nz_m(\Sigma^{-1})_{ik}(\Sigma^{-1})_{nj}$$
$$=-(\Sigma^{-1})_{lj}\sum_{m,n}(\Sigma^{-1})_{kn}z_nz_m(\Sigma^{-1})_{im}-(\Sigma^{-1})_{ik}\sum_{m,n}(\Sigma^{-1})_{nj}z_nz_m(\Sigma^{-1})_{ml}.$$
$$=-\Sigma^{-1}\otimes[\Sigma^{-1}Z\Sigma^{-1}]-[\Sigma^{-1}Z\Sigma^{-1}]\otimes\Sigma^{-1},$$
where we made use of $\Sigma=\Sigma^\top$ and followed the established convention to determine the order multiplication for each component. Putting it all together provides
$$\nabla_\Sigma\nabla_\Sigma p(z)=\frac{p(z)}{4}\left[(\Sigma^{-1}Z\Sigma^{-1})^{\otimes2}+3\Sigma^{-\otimes2}-3\Sigma^{-1}\otimes(\Sigma^{-1}Z\Sigma^{-1})-3(\Sigma^{-1}Z\Sigma^{-1})\otimes\Sigma^{-1}\right].$$
