I'm trying to perform the Eigenvalue Estimation algorithm to the Hubbard model with two sites, one spin-up fermion and one spin-down fermion, with Qiskit. Given the Hamiltonian:

$\hat{H}=u\sum_{i=1}^N\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}-t\sum_{i=1}^N\sum_{s=\uparrow,\downarrow}(a_{is}^\dagger a_{i+1\:s}+a_{i+1\:s}^\dagger a_{is})$

I obtained the analytical solution $\epsilon(x)=\frac{1}{2}(x-\sqrt{x^2+16})$ for the ground state energy $E(x)=t\epsilon(x)$ from diagonalization, where $x=u/t$, $u$ is the single-site repulsive energy and $t$ is the hopping energy term. I also found the following eigenstate:

$|\psi(x)\rangle=-\frac{1}{2}\frac{x-\sqrt{x^2+16}}{\sqrt{x^2+16-x\sqrt{x^2+16}}}|\psi_c\rangle+\frac{2}{\sqrt{x^2+16- x\sqrt{x^2+16}}}|\psi_u\rangle$

where $|\psi_c\rangle=|1\uparrow1\downarrow\rangle+|2\uparrow2\downarrow\rangle$ and $|\psi_u\rangle=|1\uparrow2\downarrow\rangle+|1\downarrow2\uparrow\rangle$. I mapped the qubits as $|q_0\rangle=|1\uparrow\rangle$, $|q_1\rangle=|1\downarrow\rangle$, $|q_2\rangle=|2\uparrow\rangle$ and $|q_3\rangle=|2\downarrow\rangle$, plus the measurement qubits and applied Jordan-Wigner transformation to the Hamiltonian. I performed Trotter-Suzuki decomposition at the second order for the time evolution operator


where $\theta$ is now the adimensional time step and $\hat{h}=\hat{H}/t$. The function for creating the time-evolution operator I implemented in Qiskit is:

def evolutionGate(n,theta,x):
tmp = QuantumCircuit(n)
for k in range(n):
for k in range(n//2):
for k in range(n//2):
for k in range(n//2):
for k in range(n):
return tmp.to_gate()

Then I defined the function for creating the circuit of the QPE algorithm:

def QPE(x, theta):
n_qub, n_meas, shots = 4, 6, 2048
qc = QuantumCircuit(n_qub + n_meas, n_meas)
alpha = -(x - mt.sqrt(x**2 + 16))/(2*mt.sqrt(x**2 + 16 - x*mt.sqrt(x**2 + 16)))
beta = 2/mt.sqrt(x**2 + 16 - x*mt.sqrt(x**2 + 16))
state = [0, 0, 0, beta, 0, 0, alpha, 0, 0, alpha, 0, 0, beta, 0, 0, 0]
qc.initialize(state, range(n_qub))
for i in range(n_meas):
    qc.h(n_qub + i)
    ev = QuantumCircuit(n_qub)
    for j in range(2**i):
        ev.append(evolutionGate(n_qub, theta, x), range(n_qub))
    qc.append(ev.to_gate().control(1), [n_qub + i, 0, 1, 2, 3])
qft = QFT(num_qubits = n_meas, inverse = True).to_gate()
qc.append(qft, qargs = qc.qubits[n_qub:n_qub + n_meas])
qc.measure(qc.qubits[n_qub:n_qub + n_meas], range(n_meas))
counts = execute(qc, backend = Aer.get_backend('qasm_simulator'), shots = shots).result().get_counts()
mean = sum(int(key, 2)*value/shots for key, value in counts.items())
phase = 2*np.pi*mean/(2**n_meas)
return phase

In the end, I wanted to perform the QPE by changing the range of the value $x$. If $\phi\in[0,1)$ is the number I get from the algorithm, I understood that the phase measured is $2\pi\phi$ (plus a non-observable phase of the form $2k\pi$, $k$ integer) and then I confronted it with the theoretical phase $-\theta\epsilon(x)$. I have done this with the following code:

x, x_max, theta = 0, 1, 0.01
while x <= x_max:
phase = QPE(x,theta)
phase_th = -theta*(x - np.sqrt(x**2 + 16))/2
print('Phase measured =', phase, 'for x =', x)
print('Phase theoretical =', -theta*(x-mt.sqrt(x**2+16))/2, 'for x =', x)
x = x + 0.1

However, the results I get from the algorithm are very different from the ones theoretical. I would like to ask your help because I can't see the mistake and I'm really struggling with this. Also, I'm new in the platform so I'm sorry if I did something wrong with the posting.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.