So I understand that the delta method allows us to approximate the expectation and variance of a random variable. Mathematically,
Let $Z = g(X,Y)$ where $X$ and $Y$ are random variables with $E[X] = \mu_X$ and $E[Y] = \mu_Y$ and $Var(X) = \sigma_X$ and $Var(Y) = \sigma_Y$. We do not know anything else about the distribution of $X$ and $Y$.
So approximating the $E[Z]$ and $Var(Z)$
Through Taylor's approximation: $$g(X,Y) = g(\mu_x, \mu_y) + (X-\mu_X)\frac{d}{dX}g(\mu_x, \mu_y) + (Y-\mu_Y)\frac{d}{dY}g(\mu_x, \mu_y) + \frac{1}{2}(X-\mu_X)^2\frac{d^2}{dX^2}g(\mu_x, \mu_y) + (X-\mu_X)(Y-\mu_Y)\frac{d^2}{dXdY}g(\mu_x, \mu_y) + \frac{1}{2}(Y-\mu_Y)^2\frac{d^2}{dY^2}g(\mu_x, \mu_y) + \text{higher order terms}$$ We'll drop the higher order terms because we don't know anything about the higher moments of $X$ or $Y$.
So, $$E[Z] = g(\mu_x, \mu_y) + E[(X-\mu_X)]\frac{d}{dX}g(\mu_x, \mu_y) + E[Y-\mu_Y)]\frac{d}{dY}g(\mu_x, \mu_y) + \frac{1}{2}E[(X-\mu_X)^2]\frac{d^2}{dX^2}g(\mu_x, \mu_y) + E[(X-\mu_X)(Y-\mu_Y)]\frac{d^2}{dXdY}g(\mu_x, \mu_y) + \frac{1}{2}E[(Y-\mu_Y)^2]\frac{d^2}{dY^2}g(\mu_x, \mu_y)$$
Notice that $E[Y - \mu_Y] = E[X - \mu_X]= 0$ by linearity. Therefore,
$$E[Z] = g(\mu_x, \mu_y) + \frac{1}{2}Var(X)\frac{d^2}{dX^2}g(\mu_x, \mu_y) + Cov(X,Y)\frac{d^2}{dXdY}g(\mu_x, \mu_y) + \frac{1}{2}Var(Y)\frac{d^2}{dY^2}g(\mu_x, \mu_y)$$
Okay, now I want to write this in linear algebra notation, but when I do this I lose the clear functions "Var()" and "Cov()." I don't know if I'm just formulating this poorly and if I am please let me know
$$E[Z] = g(\mu_X, \mu_Y) + E \left[ \begin{pmatrix} (X - \mu_X) & (Y - \mu_Y)\\ \end{pmatrix} \frac{1}{2}\begin{pmatrix} \frac{d^2}{dX^2} & \frac{d^2}{dXdY} \\ \frac{d^2}{dXdY} & \frac{d^2}{dY^2} \end{pmatrix} \begin{pmatrix} (X - \mu_X) \\ (Y - \mu_Y) \end{pmatrix} \right] $$
I guess my question is can I actually do something here? Some matrix notation that I don't know. I would really love to use $COV()$ and $Var()$ like I did in the scalar formula. Did I just completely mess up the linear algebra formulation of the multivariate method?
I'd appreciate any help :-)
$\DeclareMathOperator{\tr}{\operatorname{tr}} \DeclareMathOperator{\Var}{\operatorname{Var}}$
Ignoring several issues I have with the exposition of your question (e.g. the equations should be approximations, the Hessian is not written correctly, and the derivatives are expressed with respect to random variables instead of the arguments of the function), I think the substance of your question is how to write the second order moment expressions in terms of variance or covariance matrices.
You could use traces. So let $Z=(X-\mu_x, Y-\mu_Y)'$ and let $H$ be half the hessian matrix. Then since we are working with scalars, and using the property $\tr(AB)=\tr(BA)$, we have $$\small E[Z'HZ]=E[\tr(Z'HZ)]=E[\tr(HZZ')]=\tr(E[HZZ'])=\tr(HE[ZZ'])=\tr(H\Var((X,Y)')).$$
where $\Var((X,Y)')$ denotes the variance matrix of column random vector $(X,Y)'$.