A question about qiskit's WeightedPauliOperator.evolve method for IPEA (IQPE) implementation

Question

I'm doing an implementation/tutorial of the iterative phase estimation algorithm (IPEA) and have a problem. My Hamiltonian is a diagonal matrix $H = \begin{pmatrix} E_1 & 0 \\ 0 & E_2\end{pmatrix}$, I want to estimate both eigenvalues. When I do it instead of $E_1$ the estimated value is 0, and instead of $E_2$ the estimated value is $E_2 - E_1$. Is it possible that when I call H.evolve(t=1) method it doesn't create a circuit for $\mathrm{e}^{-iHt}$, but instead, creates a circuit for $\mathrm{e}^{-i\hat{H}t}$ , where $\hat{H} = \begin{pmatrix}0 & 0 \\ 0 & E_2 - E_1\end{pmatrix}$?

In other words, if I have $H = \begin{pmatrix} E_1 & 0 \\ 0 & E_2\end{pmatrix}$, what should do the H.evolve(t=1) method?

a) Create a circuit for $\mathrm{e}^{-iHt}$ unitary operator.

b) Create a circuit for $\mathrm{e}^{-i\hat{H}t}$ unitary operator.

This is the link to my tutorial and the short version of it is presented here:

In the tutorial I am writing

E_1, E_2 = (0, random())

and it finds $E_2$, because when $E_1 = 0$, $E_2 = E_2 - E_1$.

Here in the following code, this line is changed:

E_1, E_2 = (random(), random())

and the algorithm estimates $E_2 - E_1$ value instead of $E_2$.

import numpy as np
from random import random
from qiskit import *
from qiskit.aqua.operators import WeightedPauliOperator, MatrixOperator
from qiskit.aqua.operators.op_converter import to_matrix_operator
from qiskit.aqua.utils.controlled_circuit import get_controlled_circuit

backend = BasicAer.get_backend('qasm_simulator')
q = QuantumRegister(1)
a = QuantumRegister(1)
c = ClassicalRegister(1)


def simple_hamiltonian_operator(E_1, E_2):
    """
    Creates 0.5*(E_1 + E_2)*I + 0.5*(E_1 - E_2)*Z pauli sum
    that will be our "simple" Hamiltonian operator. The corresponding
    matrix for Hamiltonian is [[E_1, 0], [0, E_2]].

    """
    pauli_dict = {
        'paulis': [{"coeff": {"imag": 0.0, "real": 0.5 * (E_1 - E_2)}, "label": "Z"},
                   {"coeff": {"imag": 0.0, "real": 0.5 * (E_1 + E_2)}, "label": "I"}
                   ]
    }
    return WeightedPauliOperator.from_dict(pauli_dict)


E_1, E_2 = (random(), random())
print("We want to estimate E_2 = {}".format(E_2))

H = simple_hamiltonian_operator(E_1, E_2)
print("The Hamiltonian in matrix form:")
print(to_matrix_operator(H).dense_matrix)

t = 1
H_circuit = H.evolve(evo_time=t, quantum_registers=q)
H_circuit.data.__delitem__(-1) # deleting a barrier at the end of circuit

# control version of the circuit
control_H = get_controlled_circuit(H_circuit, a[0])

num_bits_estimate = 10
phase = 0
for k_precision in reversed(range(num_bits_estimate)):
    # Create a Quantum Circuit acting on the q register
    k_circ = QuantumCircuit(q, a, c)

    # (1) |1> eigenstate initialization
    k_circ.x(q[0])

    # (2) Initial Hadamard gate applied on ancillary qubit.
    k_circ.h(a[0])

    # (3) The control Hamiltonian applied on the qubits where control qubit is the ancillary qubit.
    for order in range(2 ** k_precision):
        k_circ += control_H

    # (4) The phase gate and final Hadamard gate on ancillary qubit.
    phase_shift = 2 * np.pi * phase * 2 ** k_precision
    k_circ.u1(-phase_shift, a[0])
    k_circ.h(a[0])

    # (5) Measurement of ancillary qubit (findig the bit)
    k_circ.measure(a[0], c[0])

    # (6) executing on Quantum Computer and finding a bit from the phase
    job = execute(k_circ, backend, shots=1000)
    result = job.result()
    counts = result.get_counts()
    value = int(max(counts, key=counts.get))

    # (7) phase after this iteration
    phase += value / 2 ** (k_precision + 1)

# obtaining E_2 from the estimated phase
if E_2 > E_1:
    eigenvalue = 2 * np.pi * (1 - phase) / t
else:
    eigenvalue = -2 * np.pi * phase / t

print("Eigenvalue of Hamiltonian that we want to estimate: E_2 = " + str(E_2))
print("Meanwhile, E_2 - E_1 = {}".format(E_2 - E_1))
print("Estimated eigenvalue of Hamiltonian: " + str(eigenvalue))