Inverse Fourier Transform
Suppose in time domain a Gaussian kernel is expressed as:
f(x,y) = exp (-r.^2/r0.^2);

where r.^2 = x.^2 + y.^2;, r0 is the radius of the exp(-1) point.
Then the 2D FFT should be:
F(u,v) = pi/r0.^2 * exp(-w.^2/w0.^2);

where w.^2 = u.^2 + v.^2; , w0 = 1/ (pi * r0);
As a result, I would change your function to:
function H = gaussian_kernel(size, sigma)

    half_size = floor(size/2);
    [u,v]=meshgrid(-half_size:half_size,-half_size:half_size);
    H = exp(-(u.^2+v.^2) ./ (2 * sigma^2));
    r0 = 2 * sigma^2 * pi^2;
    H = H * r0 * pi;
end

And,
H = gaussian_kernel(16, 2);
subplot(2,1,1),imagesc(H)  % frequency domain
subplot(2,1,2),imagesc(real(fftshift((ifft2(fftshift(H))))))  % time domain

result:

Suppose it is (-N/2+1 : N/2) /N * fs in the frequency axis (N is the sampling point number, and fs is the sampling rate), then it is supposed to be (0:N-1)/(N * fs) in spatial axis.