# Problems when computing the two-electron reduced density matrix (2-RDM) by qiskit

The definition of a 2-RDM is $$D_{ijkl}=\langle\psi|a_i^{\dagger}a_j^{\dagger}a_k a_l|\psi\rangle,i,,j,k,l\in[0,N-1]$$, where $$N$$ is the number of qubits/ orbitals in the quantum system.I am trying to use the following code to get the 2-RDM for a given state and the eigenvalue for the 2-RDM,

from qiskit import *
import numpy as np

#Operator Imports
from qiskit.opflow import Z, X, I

#Circuit imports
from qiskit_nature.circuit.library import HartreeFock
from qiskit import Aer
from qiskit_nature.mappers.second_quantization import JordanWignerMapper
from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit.algorithms import VQE, NumPyMinimumEigensolver
import matplotlib.pyplot as plt
import matplotlib
from qiskit.tools.visualization import circuit_drawer
from qiskit_nature.operators.second_quantization import FermionicOp
from qiskit.opflow import Z, X, I, StateFn, CircuitStateFn, SummedOp
from qiskit.opflow.converters import CircuitSampler
from qiskit.utils import QuantumInstance
from qiskit.opflow import expectations
import itertools
matplotlib.use('Agg')

mapper = JordanWignerMapper()
converter = QubitConverter(mapper=mapper, two_qubit_reduction=False)
num_particles = (2, 2)
num_spin_orbitals = 6
init_state = HartreeFock(num_spin_orbitals, num_particles, converter)

init_state1 = QuantumCircuit(num_spin_orbitals)

init_state1.h(0)
init_state1.h(1)
print(init_state1)

def get_rdm(init_state):
two_rdm = np.zeros((num_spin_orbitals,) * 4)
for i, j, k, l in itertools.product(range(num_spin_orbitals), repeat=4):
# if i != j and k != l:
s = "+_{i} +_{j} -_{k} -_{l}".format(i=str(i),j=str(j),k=str(k),l=str(l))
fermi_term = FermionicOp(s, register_length=num_spin_orbitals)
qubit_term = converter.convert(fermi_term, num_particles=num_particles)

# Evaluate the Hamiltonian term w.r.t. the given state
# temp = (~init_state @ qubit_term @ init_state).eval()
temp = ~StateFn(qubit_term) @ CircuitStateFn(primitive=init_state, coeff=1.)
temp = int(temp.eval())

two_rdm[i,j,k,l] = temp

return two_rdm

# rdm1 = rdm(init_state)
# print('--------------')
two_rdm = get_rdm(init_state1)

# Transpose of the 2-RDM
two_rdm = np.transpose(two_rdm, (0, 1, 3, 2))
two_rdm = two_rdm.reshape((num_spin_orbitals ** 2, num_spin_orbitals ** 2))
# SVD of the 2-RDM
a, b = np.linalg.eig(two_rdm)
print(a)

u, s, vh = np.linalg.svd(two_rdm, full_matrices=True)
print(s)
print(s.sum())


For the Hartree Fock state, the eigenvalues of the corresponding 2-RDM output six non-zero elements, which all are $$2$$. In my opinion, the eigenvalues of the 2-RDM stand for the meaning that probability of two electron residing in the corresponding two orbitals. Therefore, the eigenvalue should be from zero to one.