What is the convention of indices for the one and two-body integrals in qiskit

I was able to calculate the one and two-body integrals for the H2 with the following code

from qiskit.chemistry.drivers import PySCFDriver, UnitsType

atom = 'H .0 .0 .0; H .0 .0 0.74'
distance_unit = UnitsType.ANGSTROM
basis = 'sto3g'

driver = PySCFDriver(atom, unit=distance_unit, basis=basis)

molecule = driver.run()

h1 = molecule.one_body_integrals
h2 = molecule.two_body_integrals
print( h1 )
print( h2 )


and the results is

[[-1.25330979  0.          0.          0.        ]
[ 0.         -0.47506885  0.          0.        ]
[ 0.          0.         -1.25330979  0.        ]
[ 0.          0.          0.         -0.47506885]]
[[[[-0.33737796  0.          0.          0.        ]
[ 0.         -0.09060523  0.          0.        ]
[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]]

[[ 0.         -0.09060523  0.          0.        ]
[-0.3318557   0.          0.          0.        ]
[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]]

[[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]
[-0.33737796  0.          0.          0.        ]
[ 0.         -0.09060523  0.          0.        ]]

[[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]
[ 0.         -0.09060523  0.          0.        ]
[-0.3318557   0.          0.          0.        ]]]

[[[ 0.         -0.3318557   0.          0.        ]
[-0.09060523  0.          0.          0.        ]
[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]]

[[-0.09060523  0.          0.          0.        ]
[ 0.         -0.34882575  0.          0.        ]
[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]]

[[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]
[ 0.         -0.3318557   0.          0.        ]
[-0.09060523  0.          0.          0.        ]]

[[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]
[-0.09060523  0.          0.          0.        ]
[ 0.         -0.34882575  0.          0.        ]]]

[[[ 0.          0.         -0.33737796  0.        ]
[ 0.          0.          0.         -0.09060523]
[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]]

[[ 0.          0.          0.         -0.09060523]
[ 0.          0.         -0.3318557   0.        ]
[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]]

[[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]
[ 0.          0.         -0.33737796  0.        ]
[ 0.          0.          0.         -0.09060523]]

[[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]
[ 0.          0.          0.         -0.09060523]
[ 0.          0.         -0.3318557   0.        ]]]

[[[ 0.          0.          0.         -0.3318557 ]
[ 0.          0.         -0.09060523  0.        ]
[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]]

[[ 0.          0.         -0.09060523  0.        ]
[ 0.          0.          0.         -0.34882575]
[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]]

[[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]
[ 0.          0.          0.         -0.3318557 ]
[ 0.          0.         -0.09060523  0.        ]]

[[ 0.          0.          0.          0.        ]
[ 0.          0.          0.          0.        ]
[ 0.          0.         -0.09060523  0.        ]
[ 0.          0.          0.         -0.34882575]]]]


This immediately raises the question of how qiskit labels the two-body integrals. If we look at the definition of second quantized Hamiltonian

$$H = \sum_{i,j}h_{ij}a^\dagger_i a_j + \frac{1}{2}\sum_{i,j,k,l}h_{ijkl}a^\dagger_ia^\dagger_ja_ka_l$$

then we realize that $$h_{0000}=0$$, whereas qiskit gives -0.33737796. The answer for the integrals is given in the following table. Any helps are really appreciated.

Something to keep in mind is that integrals often hold two forms of symmetries - not only the symmetries enumerated in that table between the indices, but also the spins of the orbitals. I think this is explained well in Microsoft's documents. For example, if the input file notes $$a^\dagger_0 a^\dagger_1 a_1 a_0$$, we need to consider the potential permutations and the potential spins. So, we'd consider $$a^\dagger_{0, \sigma} a^\dagger_{1, \rho} a_{1, \rho} a_{0, \sigma}$$, and all of the permutations of the numbers, and all combinations of $$\sigma, \rho \in \{ \uparrow, \downarrow \}$$, where $$\{ \uparrow, \downarrow \}$$ are the spins.
So, while the $$a^\dagger_0 a^\dagger_0 a_0 a_0$$ may be invalid, $$a^\dagger_{0, \uparrow}, a^\dagger_{0, \downarrow} a_{0, \downarrow} a_{0, \uparrow}$$ is totally valid after indexing.