There seems to be an easy way to get Hessian in a batched manner:
All you need to do is ensure that your function does not operate across samples in a batch
Here is a code snippet:
import torch
from torch.autograd.functional import jacobian

p=3 # dimension
n=10 # number of samples in the batch
x = torch.randn(n, p, requires_grad=True)

Q_12 = torch.randn(p, p)
Q = Q_12.T @ Q_12 # random symmetric positive-semidefinite matrix

# quadratic function. 0.5*xT Q x where x is every row of input tensor
f = lambda x: 0.5*torch.tensordot(x, x @ Q) 

def get_sum_of_gradients(x):
    h = f(x)
    h.backward(create_graph=True)
    return x.grad.sum(0)

Hessian_per_sample = jacobian(get_sum_of_gradients, x).swapaxes(0, 1)

# Check if Hessian for every sample is Q
print(torch.all(torch.isclose(Hessian_per_sample, torch.tile(Q, (n, 1, 1)))).item())