You can do this using D[]. Given a function with (say) arguments x and y you can write:

D[f[x,y], {{x,y},2}]

You can use Grad[] twice, of course, but D[] will get you there directly.

Some background info here:

https://mathworld.wolfram.com/Hessian.html

On that page, they present an implementation of a Hessian operator:

HessianH[f_, x_] := D[f, {x,2}]

You'd call it like this:

HessianH[f[x,y], {x,y}]

If you don't want to have to repeat the argument list, you could do something like:

H[f_, x_]:= HessianH[Apply[f, x], x]

So now you can say:

H[f, {x,y}]

As of v12, you can access a pre-defined version of this function as follows:

ResourceFunction["HessianMatrix"][f[x,y], {x,y}]

which emits a matrix in the normal way. More info here:

https://resources.wolframcloud.com/FunctionRepository/resources/HessianMatrix