# QPE algorithm to the Hubbard Model

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

$$U=\exp(-\frac{i\tau\hat{H}}{\hslash})=\exp{(-i\theta\hat{h})}$$

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):
tmp.rx(-theta*x/4,k)
for k in range(n//2):
tmp.cx(2*k,2*k+1)
tmp.rz(-theta*x/4,2*k+1)
tmp.cx(2*k,2*k+1)
for k in range(n//2):
tmp.cx(k,k+2)
tmp.rx(-2*theta,k+2)
tmp.cx(k,k+2)
tmp.cx(k+2,k)
tmp.rx(-2*theta,k)
tmp.cx(k+2,k)
for k in range(n//2):
tmp.cx(2*k,2*k+1)
tmp.rz(-theta*x/4,2*k+1)
tmp.cx(2*k,2*k+1)
for k in range(n):
tmp.rx(-theta*x/4,k)
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('------------------------------------------------------------------')
print('Phase measured =', phase, 'for x =', x)
print('Phase theoretical =', -theta*(x-mt.sqrt(x**2+16))/2, 'for x =', x)
print('------------------------------------------------------------------')
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.