# Applying Trotterization to a Hamiltonian for Time Evolution in Qiskit

I'm currently working on a project where I need to simulate the time evolution of a quantum system using Qiskit. The Hamiltonian of my system is given by:

$$H = -J \sum_{j=1}^{N-1} (\sigma_j^x \sigma_{j+1}^x + \sigma_j^y \sigma_{j+1}^y) + U \sum_{j=1}^{N-1} \sigma_j^z \sigma_{j+1}^z + \sum_{j=1}^{N} h_j \sigma_j^z$$

I want to use the Trotter-Suzuki decomposition to approximate the time evolution operator, which is given by: $$e^{iH\Delta t} = \prod_j A^j \prod_{j \, \text{even}} B^j \prod_{j \, \text{odd}} B^j + O((\Delta t)^2)$$ where $$A^j = e^{-ih_j\sigma^z_j \Delta t}$$ and $$B^j = e^{-i(U\sigma^z_j\sigma^z_{j+1} - J(\sigma^x_j\sigma^x_{j+1} + \sigma^y_j\sigma^y_{j+1}))\Delta t}$$

In a paper that I'm working on they shown this scheme as the following:

I've tried to solve this by using the following code

The part where I've created my Hamiltonian and its corresponding circuit:

def construct_TFIM_circuit(N_qubit=6,initial_spin_config="domain_wall",
Type="XX_chain",J=1,t=1):
qc = QuantumCircuit(N_qubit)
if Type=="XX_chain":
h=np.random.rand(N_qubit)*0
U=0

if initial_spin_config=="neel":
for j in range(N_qubit):
if j % 2 == 0:
qc.x(j)
elif initial_spin_config=="domain_wall":
for j in range(N_qubit):
if j < (N_qubit/2):
qc.x(j)

state_vector=[]
loc_magnetization=[]
Infos={}

for j in range(N_qubit):
# Build the evolution gate
H=  (U * Z^Z ) - (J*X^X)-(J*Y^Y)
evo_H= PauliEvolutionGate(H, time=t)
if j!=N_qubit-1:
qc.append(evo_H, [j,j+1])
for j in range(N_qubit):
qc.append(PauliEvolutionGate(h[j]*Z, time=t), [j])



It seems like it crates the desired circuit since it matches with the given figure

Then I tried to get local magnetization expectations with the variable loc_magnetization for jth spin by the formula $$M_j(t) = \langle \psi(t) | \sigma^z_j | \psi(t) \rangle$$ with the code part(in the same function):

state_vector=Statevector.from_instruction(qc)
Exact=[]
Sampled=[]
for j in range(N_qubit):
Op_circ=QuantumCircuit(N_qubit)
Op_circ.z(j)
op=CircuitOp(Op_circ)

psi = CircuitStateFn(qc)
psi._statevector = state_vector   # Set the statevector manually

# Exact & Sampled expectations
(M_j_exact,M_j_sampled)=find_expectation(op,psi)
# loc_magnetization have (Exact,Sampled) pair values
loc_magnetization.append((M_j_exact,M_j_sampled))

Infos["JT"]=J*t
Infos["State"]=state_vector
Infos["loc_magnetization"]=loc_magnetization
return (qc,Infos)


However my result doesn't match with the desired expectation, can anybody tell where I'm doing wrong?