FFT of shifted image in Numpy
In the code the OP is comparing the real and the imaginary components of the FFT in each to determine similarity. If a phase shift occurs then the real and imaginary components will also change. What the OP should be comparing is the magnitudes of the complex values in each sample, which should be invariant under such image shifts.
Confirming this; the OP had done the following to compare the FFT results:
plt.bar(['1-1','2-1','3-1'],[np.mean(np.square(H1.imag-H1.imag)),np.mean(np.square(H1.imag-H2.imag)),np.mean(np.square(H1.imag-H3.imag))])

resulting in:

If one does the following instead:
plt.bar(['1-1','2-1','3-1'],[np.mean(np.square(np.abs(H1)-np.abs(H1))),np.mean(np.square(np.abs(H1)-np.abs(H2))),np.mean(np.square(np.abs(H1)-np.abs(H3)))])

We get a result that we expect:

The equivalent magnitude is given by the translation property of the 2D DFT where:
$$x[k-m, l-n] \mapsto X[u,v]e^{-j2\pi (mu/M + nv/N)}$$
This is similar to the time shift property of the 1D DFT (which shows fundamentally what is occurring and can be extended to a 2D DFT):
Given 1D DFT as
$$x[n] \mapsto \sum_{n=0}^{N-1}x[n]W_N^{nk} = X[k]$$
Where $W_N^n=e^{-j2\pi n/N}$
From which:
$$x[n-m] \mapsto \sum_{n=0}^{N-1}x[n-m]W_N^{nk}$$
Substitute $u= n-m$
$$\sum_{u = -m}^{N-1-m}x[u]W_N^{(u+m)k}  = \sum_{u = -m}^{N-1-m}x[u]W_N^{uk}W_N^{mk}$$
$$=W_N^{mk}\sum_{u = -m}^{N-1-m}x[u]W_N^{uk}$$
for any k:
$$\sum_{u = -m}^{N-1-m}x[u]W_N^{uk} = \sum_{u=0}^{N-1}x[u]W_N^{uk}$$
Therefore
$$x[n-m] \mapsto W_N^{mk}X[k]$$
And $W_N^{mk}$ as $e^{-j2\pi mk/N}$ only modifies phase, not magnitude
Also the DFT is circularly periodic in the time domain, so any such rotation must done in circular fashion;  x[mod(n-m,M)] = x [n]