I have the following function (just used as an example):
$y_t=g(y_{t-1},\epsilon_t, \sigma)$
of which I have the following second-order Taylor expansion around a point such that $y=\bar{y}, \ \epsilon_t=0=\sigma$:
$y_t=g_y(y_{t-1}-\bar{y})+ g_{\epsilon}\epsilon_t+g_{\sigma}\sigma+0.5g_{yy}(y_{t-1}-\bar{y})\otimes(y_{t-1}-\bar{y})+g_{\epsilon y}(y_{t-1}-\bar{y})\otimes \epsilon_t+g_{\sigma \epsilon}(y_{t-1}-\bar{y})\sigma+0.5g_{\epsilon \epsilon}(\epsilon_t \otimes \epsilon_t)+g_{\sigma \epsilon}\sigma \epsilon_t+0.5g_{\sigma \sigma}\sigma^2$
where $y_t$ and $\epsilon_t$ are vectors of conformable dimension and the pedices in $g$ denote partial derivatives. $\sigma$ is instead a scalar.
My struggle is in understanding why the notation uses Kronecker product instead of simple square of the type $(y_{t-1}-\bar{y})'(y_{t-1}-\bar{y})$ for the part related to the second-order derivatives.
For a vector-valued function $f(\mathbf z): \mathbb{R}^n \rightarrow \mathbb{R}$, the Taylor series up to the second order is given by
$$ f(\mathbf z) \approx f(\mathbf z_0) + \frac{\partial f}{\partial \mathbf z}^{\rm T} (\mathbf z - \mathbf z_0) + \frac 12 (\mathbf z - \mathbf z_0)^{\rm T} \frac{\partial^2 f}{\partial \mathbf z \mathbf z^{\rm T}} (\mathbf z - \mathbf z_0),$$
which we can alternatively write as
$$ f(\mathbf z) \approx f(\mathbf z_0) + \frac{\partial f}{\partial \mathbf z}^{\rm T} (\mathbf z - \mathbf z_0) + \frac 12 {\rm vec}\left( \frac{\partial^2 f}{\partial \mathbf z \mathbf z^{\rm T}}\right)\Big( (\mathbf z - \mathbf z_0) \otimes (\mathbf z - \mathbf z_0)\Big),$$
where $\frac{\partial^2 f}{\partial \mathbf z \mathbf z^{\rm T}}$ represents the matrix of second derivatives of $f$ with respect to its arguments and ${\rm vec}\left(\cdot\right)$ the vectorization operator.
Why do we need this form and not simply something like $\mathbf b^{\rm T}\left(\mathbf z-\mathbf z_0\right)^2$, where the $(\cdot)^2$ is element-wise? Simply, because this way we would be missing all the cross terms.
Consider the example $n=2$ and $\mathbf z = [x, y]^{\rm T}$. Then, $\mathbf b^{\rm T}\left(\mathbf z-\mathbf z_0\right)^2 = b_1 (x-x_0)^2 + b_2 (y-y_0)^2$, however
$$ \frac 12 (\mathbf z - \mathbf z_0)^{\rm T} \frac{\partial^2 f}{\partial \mathbf z \mathbf z^{\rm T}} (\mathbf z - \mathbf z_0) = \frac 12 f_{xx}(\mathbf z_0)\cdot (x-x_0)^2 + \frac 12 f_{yy}(\mathbf z_0)\cdot (y-y_0)^2 + \color{red}{f_{xy}(\mathbf z_0)\cdot(x-x_0)(y-y_0)}. $$
The red part would be missing if you did not use the Kronecker product relation. You can skip it only if you know that the mixed derivative $f_{xy}$ is zero at $\mathbf z_0$, which it will in general not be.