I solved this integral using Mathematica:
$H(z,t)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}dxdy(x\cdot y)^aE_{-t}(-ix)E_{-t}(-iy)E_{-z}(ix)E_{-z}(iy) $
Where $E_{n}(x)$ is the Exponential Integral function (ExpIntegralE[n,x] in Mathematica)
And the solution was:
$\frac{2\pi^2 i}{(1-a+z+t)^2sin(\pi a)}(\frac{\Gamma(a-z)^2}{\Gamma(-z)^2}+\frac{\Gamma(a-t)^2}{\Gamma(-t)^2}-\frac{2e^{i\pi a}\Gamma(a-z)\Gamma(a-t)}{\Gamma(-z)\Gamma(-t)})$
As you can see, this integral satisfy this symmetry: $H(z,t)=H(t,z)^*$ When I assumed that a,t,z are real (this assumption was included in Mathematica). But weirdly enough the solution does not satisfy this symmetry.
What happen?
Few notes:
the only condition I got was $a>-1$
I know a little bit more about $a$. It has the value of $a=C+z+t$ where $C$ is just a real constant that I can freely fix.
the solution above was not written like this in Mathematica. I used the identity $\Gamma(-x)\Gamma(1+x)sin(\pi x)=-\pi$ to write it in this form.
Moreover,this integral cannot be solved separately, you get different answers (for more info about it you can ask).
The answer is, in fact, correct. You only need to evaluate the integral $$I(z, t) = \int_0^\infty x^a E_{-z}(i x) E_{-t}(-i x) dx = \\ \int_0^{i \infty} x^a E_{-z}(i x) E_{-t}(-i x) dx = i^{a + 1} \int_0^\infty x^a E_{-z}(-x) E_{-t}(x) dx,$$ which can be done by using the general formula for the integral of a product of two Meijer G-functions: $$I(z, t) = i^{a + 1} \int_0^\infty x^a G_{1, 2}^{2, 0} \left(x \middle| {-t \atop -t - 1, 0} \right) G_{1, 2}^{2, 0} \left(-x \middle| {-z \atop -z - 1, 0} \right) dx = \\ i^{a + 1} G_{3, 3}^{2, 2} \left(-1 \middle| {-a, t + 1 - a, -z \atop 0, -z - 1, t - a} \right),$$ which reduces to a sum of hypergeometric functions with the argument equal to $1$ and then, by Gauss's theorem, to gamma functions.
Then, integrating over each of the four quadrants separately, we obtain $$H(z, t) = I(z, t)^2 + 2 (-1)^a I(z, t) I(t, z) + I(t, z)^2,$$ which simplifies to the solution that you provided. $I(z, t) = I(t, z)^*$, but the relation doesn't hold for $H(z, t)$.
The convergence conditions, determined by the behavior of the exponential integral at zero and at infinity, are $$ -1 < a < 1 \land \max(t, -1) + \max(z, -1) < a - 1.$$