I have these expressions, which are the gradient of a number of constraints on an optimization problem
$$\frac{\partial c}{\partial z_r}=B(I\otimes z_r + z_r\otimes I - iI\otimes z_i + iz_i\otimes I) \\ \frac{\partial c}{\partial z_i}=B(iI\otimes z_r - iz_r\otimes I + I\otimes z_i + z_i\otimes I).$$
where
$$ q=z_r+iz_i. $$
Furthermore, $c\in\mathbb{R}^{100}$, $(q, z_r, z_i)\in\ \mathbb{C}^{9}$, $B\in\mathbb{C}^{100\times 81}$ and $I\in\mathbb{R}^{9\times 9}$ is the identity matrix. Consequently, ($\frac{\partial c}{\partial z_r},\frac{\partial c}{\partial z_i})\in\mathbb{R}^{9\times100}$. I want to calculate the Hessian and therefore i require the gradient of the gradient (and thus $\frac{\partial^2 c}{\partial z_r^2}$ and the others).
This is what i have tried:
$$ d(X\otimes Y) = dX\otimes Y +X\otimes dY $$ and so
$$ f=B(I\otimes z_r) \\ df = B(I\otimes dz_r). $$
But i am not sure how to go to $\frac{\partial f}{\partial z_r}$ from here, since $dz_r$ is still part of the kronecker product.