Let $f:\mathbb{R}^{m\times n}\to\mathbb{R}$ be $C^3$, then according to Taylor’s formula in Banach spaces and also chapter 4 of Zeidler's functional analysis, we can write the Taylor's formula for $f$ as $$f(X)=f(A)+Df(A)H+\frac{1}{2}D^2f(A)H^2+R_2,$$ where $Df(A)$ and $D^2f(A)$ are the first and second order Fréchet derivatives and are linear and bilinear maps, respectively, from $\mathbb{R}^{m\times n}\to\mathbb{R}$. We also have the Taylor's theorem for multivariate functions (though the inputs are scalar).
Now let $F:\mathbb{R}^{m\times n}\times\mathbb{R}^{p\times q} \to\mathbb{R}$ be a smooth function. Can we approximate the function at $X\in\mathbb{R}^{m\times n}$ and $Y\in\mathbb{R}^{p\times q}$ by $$F(X+\Delta X,Y+\Delta Y)\simeq F(X,Y)+D_XF(X,Y)\Delta X^2+D_YF(X,Y)\Delta Y^2+H(X,Y,\Delta X,\Delta Y)$$ where, $$H(X,Y,\Delta X,\Delta Y)=\frac{1}{2}( D^2_{X^2} F(X,Y)\Delta X^2+D^2_{Y^2}F(X,Y)\Delta Y^2)+D^2_{XY}F(X,Y)\Delta Y\Delta X+D^2_{YX}F(X,Y)\Delta X\Delta Y,$$
and $D^2_{XY}F(X,Y)$ is a bilinear form that is the Fréchet derivative of $D_XF(X,Y)$ with respect to $Y$.
More importantly, can we investigate the critical points of $F$ similar to "Second partial derivative test" using the sign of $H$?
I think one way to prove this is to vectorize both matrices and concatenate them into a long vector and use Taylor's theorem for vector valued functions. But before getting into its details I just wanted to know if there is a source dealing with Taylor approximation of this type of functions or is there a more straightforward way to prove this.
Zeidler, Eberhard, Applied functional analysis. Main principles and their applications, Applied Mathematical Sciences. 109. New York, NY: Springer-Verlag. xvi, 404 p. (1995). ZBL0834.46003.