I want to compute this <Zi Zj> - <Zi><Zj> for an entangled n-qubit initial state under the application of a general XY Hamiltonian for a range of time 0 to T. This is to be done to study how the correlations between the qubits decay over time. This quantity should be plotted <Zi Zj> - <Zi><Zj> as a function of time from 0 to T.
Hamiltonian Code
def get_hamiltonian(params, **kwargs):
'''
Returns a general XY Hamiltonian.
'''
## Extract parameters
L, alpha, J = params['L'], params['alpha'], params['J']
### Define the XX and YY tuples with corrected list comprehension
XX_tuples = [("XX", [i, j], -J / np.power(np.abs(i - j), alpha)) for i in range(L - 1) for j in range(i + 1, L)]
YY_tuples = [("YY", [i, j], -J / np.power(np.abs(i - j), alpha)) for i in range(L - 1) for j in range(i + 1, L)]
# Combine the tuples to create the Hamiltonian
hamiltonian = SparsePauliOp.from_sparse_list([*XX_tuples, *YY_tuples], num_qubits=L)
return hamiltonian.simplify()
N-qubit Entangled State
def get_qc_for_n_qubit_GHZ_state(params, **kwargs):
"""This function will create a qiskit.QuantumCircuit (qc) for an n-qubit GHZ state.
Args:
params - only requires 'L' parameter to construct Q.circuit and entangled state.
Returns:
QuantumCircuit: Quantum circuit that generate the n-qubit GHZ state, assuming all qubits start in the 0 state
Entangled n-qubit GHZ state
"""
n = params['L']
initial_state = '0' * n
sv1 = Statevector.from_label(initial_state)
if isinstance(n, int) and n >= 2:
qc = QuantumCircuit(n)
qc.h(0)
for i in range(n-1):
qc.cx(i, i+1)
else:
raise Exception("n is not a valid input")
new_sv1 = sv1.evolve(qc)
return qc, new_sv1
Observables Construction
def observables_construction(params, **kwargs):
'''
Returns the list of observables
'''
L = params['L']
Z_op = SparsePauliOp.from_sparse_list(
[("Z", [i], 1.0) for i in range(0, L)], num_qubits=L
)
correlation_op = SparsePauliOp.from_sparse_list(
[("ZZ", [i, j], 1.0) for i in range(0, L - 1) for j in range(i+1, L)], num_qubits=L
)
return Z_op, correlation_op
Part-I
Using Qiskit Trotterization Algorithm to study <Zi Zj> - <Zi><Zj> as a function of time
params['T'] = 40.0
params['L'] = 5
params['time_step'] = 60
## Assigning Hamiltonian and Observables
H = get_hamiltonian(params)
Z_op, correlation_op = observables_construction(params)
GHZ_qc, initial_state = get_qc_for_n_qubit_GHZ_state(params)
Defining Time Evolution Problem
problem = TimeEvolutionProblem(
H,
initial_state=initial_state,
time=params['T'],
aux_operators=[Z_op, correlation_op],
)
Performing Trotterization
trotter = TrotterQRTE(num_timesteps=params['time_step'], estimator=Estimator())
result = trotter.evolve(problem)
Extracting results of the observables
observables = np.array(np.array(result.observables)[:, :, 0])
observables.shape
Plotting Data
fig, axes = plt.subplots(2, sharex=True)
times = np.linspace(0, params['T'], params['time_step'] + 1) # includes initial state
axes[0].plot(
times, observables[:, 0], label="First order", marker="x", c="darkmagenta", ls="-", lw=0.8
)
axes[1].plot(
times, np.round(observables[:, 1].astype(float), decimals = 3), label="First order", marker="x", c="darkmagenta", ls="-", lw=0.8
)
axes[0].set_ylabel(f"$<Z_i>$")
axes[1].set_ylabel(f"$<Z_iZ_j>$")
axes[1].set_xlabel("Time")
fig.suptitle("Observable evolution")
Packages used in Part-1
from qiskit.quantum_info import SparsePauliOp
from math import sin, cos, pi
import numpy as np
import matplotlib.pyplot as plt
from qiskit.quantum_info import Statevector
from qiskit_algorithms import TimeEvolutionProblem
from qiskit_algorithms import TrotterQRTE
from qiskit.primitives import Estimator
import scipy as sc
from qiskit.circuit import QuantumCircuit
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter
from qiskit.synthesis import SuzukiTrotter
from qiskit.synthesis import SuzukiTrotter
from matplotlib import cm
Queries in Part-1
The expectation values of <Z_i> are coming out to be zero. And I guess the data I obtain for these two sets of observables <Zi Zj> and <Z_i> are the summation over all possible combinations for each time step. So, I am unable to compute the required quantity <Zi Zj> - <Zi><Zj>. Could somebody please tell me how to obtain this quantity?
Part-II: Execution on Real Hardware Part-I is only for testing purposes before running on real hardware.
Initial Set-up
params['T'] = 40.0
params['L'] = 5
params['time_step'] = 60
## Assigning Hamiltonian and Observables
H = get_hamiltonian(params)
Z_op, correlation_op = observables_construction(params)
GHZ_qc, initial_state = get_qc_for_n_qubit_GHZ_state(params)
Optimizing the Problem
service = QiskitRuntimeService(channel="ibm_quantum", token=IBM_API)
backend = service.least_busy(simulator=False, operational=True, min_num_qubits=params['L'])
pm = generate_preset_pass_manager(optimization_level=1, backend=backend)
isa_circuit = pm.run(GHZ_qc)
isa_operators_list = [op.apply_layout(isa_circuit.layout) for op in correlation_op]
Execution on Hardware
options = EstimatorOptions()
options.resilience_level = 1
estimator = Estimator(backend, options=options)
Submit the circuit to Estimator
job = estimator.run([(isa_circuit, isa_operators_list)])
job_id = job.job_id()
print(job_id)
Intrepreting Results
result = job.result()[0]
values = result.data.evs
Queries in Part-2
I am not able to understand the hardware results are with respect to which quantity (for eg, with respect to time). In case they are not with respect to time, what must be done to obtain the quantity <Zi Zj> - <Zi><Zj> as a function of time?
