# Qubit Spectroscopy and Readout (by nick brønn)

Cahyati Supriyati Sangaji (My Note)

# Jaynes-Cummings Hamiltonian

For qubits, we can describe the dynamics of the qubit with the Pauli Matrices

The Pauli matrices let us define raising and lowering operators

They raise and lower qubit states

The interaction between qubit and resonator is descibed by the Jaynes-Cummings Hamiltonian

In order to understand it, we must block-diagonalize it using the Schrieffer-Wolff transformation.

# The Schrieffer-Wolff Transformation

You have a diagonalized Hamiltonian plus a perturbation

# The Schrieffer-Wolff Transformation

We seek an anti-hermitian, block off-diagonal operator 𝑆 such that

is the generalized commutator

to keep track of the order 𝜆

Then,

up to second order in 𝜆

# The Schrieffer-Wolff Transformation

to each order of 𝜆

Expanding

Going to second order yields

# The Schrieffer-Wolff Transformation

Substituting the equations needed to cancel the block off-diagonal terms and letting 𝜆→1

The recent manuscript A Systematic Method for Schrieffer-Wolff Transformation and Its Generalizations systematically provides the ansatz by calculating

and set the coefficients such that

# S-W Tranformation of the J-C Hamiltonain

For ease of tedious calculations, we will use the Python package sympy for symbolic mathematics.

# import SymPy and define symbolsimport sympy as spsp.init_printing(use_unicode=True)wr = sp.Symbol('\omega_r') # resonator frequencywq = sp.Symbol('\omega_q') # qubit frequencyg = sp.Symbol('g', real=True) # vacuum Rabi couplingDelta = sp.Symbol('Delta', real=True) # wr - wq; defined later# import operator relations and define themfrom sympy.physics.quantum.boson import BosonOpa = BosonOp('a') # resonator photon annihilation operatorfrom sympy.physics.quantum import pauli, Dagger, Commutatorfrom sympy.physics.quantum.operatorordering import normal_ordered_form

# S-W Tranformation of the J-C Hamiltonian

# Pauli matricessx = pauli.SigmaX()sy = pauli.SigmaY()sz = pauli.SigmaZ()# qubit raising and lowering operatorssplus = pauli.SigmaPlus() sminus = pauli.SigmaMinus()# define J-C Hamiltonian in terms of diagonal and # non-block diagonal termsH0 = wr*Dagger(a)*a - (1/2)*wq*sz; H1 = 0H2 = g*(Dagger(a)*sminus + a*splus); HJC = H0 + H1 + H2; HJC # print

# S-W Tranformation of the J-C Hamiltonian

# using the above method for finding the ansatzeta = Commutator(H0, H2); eta

We will need to used the methods doit(), expand, normal_ordered_form, and qsimplify_pauli

pauli.qsimplify_pauli(normal_ordered_form(eta.doit().expand()))

# S-W Tranformation of the J-C Hamiltonian

A = sp.Symbol('A')B = sp.Symbol('B')eta = A * Dagger(a) * sminus - B * a * splus;pauli.qsimplify_pauli(normal_ordered_form(Commutator(H0, eta).doit().expand()))
H2

# S-W Tranformation of the J-C Hamiltonian

S1 = eta.subs(A, g/Delta)S1 = S1.subs(B, g/Delta); S1.factor()
Heff = H0 + H1 + 0.5*pauli.qsimplify_pauli(normal_ordered_form(Commutator(H2, S1).doit().expand())).simplify(); Heff

This is typically written as

# Qiskit Pulse — Spectroscopy

Import Necessary Libraries

from qiskit.tools.jupyter import *from qiskit import IBMQIBMQ.save_account('<Your Token')IBMQ.load_account()provider = IBMQ.get_provider(hub='ibm-q', group='open', project='main')backend = provider.get_backend('ibmq_armonk')

# Verify Backend is Pulse-enabled

backend_config = backend.configuration()assert backend_config.open_pulse, "Backend doesn't support Pulse"dt = backend_config.dtprint(f"Sampling time: {dt*1e9} ns") backend_defaults = backend.defaults()

Out:

Sampling time: 0.2222222222222222 ns

# Qiskit Pulse — Spectroscopy

import numpy as np# unit conversion factors -> all backend properties returned in SI # (Hz, sec, etc)GHz = 1.0e9 # GigahertzMHz = 1.0e6 # MegahertzkHz = 1.0e3us = 1.0e-6 # Microsecondsns = 1.0e-9 # Nanoseconds# We will find the qubit frequency for the following qubit.qubit = 0# The sweep will be centered around the estimated qubit frequency.# The default frequency is given in Hz                                                                    # warning: this will change in a future releasecenter_frequency_Hz = backend_defaults.qubit_freq_est[qubit]        print(f"Qubit {qubit} has an estimated frequency of {center_frequency_Hz / GHz} GHz.")

Out:

Qubit 0 has an estimated frequency of 4.974445862780936 GHz.

# Qiskit Pulse — Spectroscopy

# scale factor to remove factors of 10 from the datascale_factor = 1e-14# We will sweep 40 MHz around the estimated frequencyfrequency_span_Hz = 40 * MHz# in steps of 1 MHz.frequency_step_Hz = 1 * MHz# We will sweep 20 MHz above and 20 MHz below the estimated # frequencyfrequency_min = center_frequency_Hz - frequency_span_Hz / 2frequency_max = center_frequency_Hz + frequency_span_Hz / 2# Construct an np array of the frequencies for our experimentfrequencies_GHz = np.arange(frequency_min / GHz,                             frequency_max / GHz,                             frequency_step_Hz / GHz)print(f"The sweep will go from {frequency_min / GHz} GHz to {frequency_max / GHz} GHz \in steps of {frequency_step_Hz / MHz} MHz.")

Out:

The sweep will go from 4.954445862780936 GHz to 4.994445862780936 GHz in steps of 1.0 MHz.

# Qiskit Pulse — Spectroscopy

# This is where we access all of our Pulse features!from qiskit import pulseinst_sched_map = backend_defaults.instruction_schedule_mapmeasure = inst_sched_map.get('measure', qubits=[qubit])x_pulse = inst_sched_map.get('x', qubits=[qubit])### Collect the necessary channelsdrive_chan = pulse.DriveChannel(qubit)meas_chan = pulse.MeasureChannel(qubit)acq_chan = pulse.AcquireChannel(qubit)# Create the base schedule# Start with drive pulse acting on the drive channelschedule = pulse.Schedule(name='Frequency sweep')schedule += x_pulse# The left shift << is special syntax meaning to # shift the start time of the schedule by some durationschedule += measure << schedule.duration# Create the frequency settings for the sweep (MUST BE IN HZ)frequencies_Hz = frequencies_GHz*GHzschedule_frequencies = [{drive_chan: freq} for freq in frequencies_Hz]

# Qiskit Pulse — Spectroscopy

schedule.draw(label=True, scaling=0.8)

# Qiskit Pulse — Spectroscopy

from qiskit import assemblefrequency_sweep_program = assemble(schedule,                                   backend=backend,                                    meas_level=1,                                   meas_return='avg',                                   shots=1024,                                   schedule_los=schedule_frequencies)# RUN the job on a real devicejob = backend.run(frequency_sweep_program)print(job.job_id())from qiskit.tools.monitor import job_monitorjob_monitor(job)# OR retreive result from previous run#job = backend.retrieve_job('<Job ID>')

Out:

<Your Jon ID>Job Status: job has successfully run

# Qiskit Pulse — Spectroscopy

frequency_sweep_results = job.result()import matplotlib.pyplot as pltplt.style.use('dark_background')sweep_values = []for i in range(len(frequency_sweep_results.results)):    # Get the results from the ith experiment    res = frequency_sweep_results.get_memory(i)*scale_factor    # Get the results for qubit from this experiment    sweep_values.append(res[qubit])plt.scatter(frequencies_GHz, np.real(sweep_values), color='white') # plot real part of sweep valuesplt.xlim([min(frequencies_GHz), max(frequencies_GHz)])plt.xlabel("Frequency [GHz]")plt.ylabel("Measured signal [a.u.]")plt.show()

# Qiskit Pulse — Spectroscopy

from scipy.optimize import curve_fitdef fit_function(x_values, y_values, function, init_params):    fitparams, conv = curve_fit(function, x_values, y_values, init_params)    y_fit = function(x_values, *fitparams)        return fitparams, y_fitfit_params, y_fit = fit_function(frequencies_GHz,                                 np.real(sweep_values),                                  lambda x, A, q_freq, B, C: (A / np.pi) * (B / ((x - q_freq)**2 + B**2)) + C,                                 [5, 4.975, 1, 3] # initial parameters for curve_fit                                )

# Qiskit Pulse — Spectroscopy

plt.scatter(frequencies_GHz, np.real(sweep_values), color='white')plt.plot(frequencies_GHz, y_fit, color='red')plt.xlim([min(frequencies_GHz), max(frequencies_GHz)])plt.xlabel("Frequency [GHz]")plt.ylabel("Measured Signal [a.u.]")plt.show()

# Qiskit Pulse — Dispersive Shift

# Create the schedules for 0 and 1schedule_0 = pulse.Schedule(name='0')schedule_0 += measureschedule_1 = pulse.Schedule(name='1')schedule_1 += x_pulseschedule_1 += measure << schedule_1.duration

# Qiskit Pulse — Dispersive Shift

schedule_0.draw()

# Qiskit Pulse — Dispersive Shift

schedule_1.draw()

# Qiskit Pulse — Dispersive Shift

frequency_span_Hz = 320 * kHzfrequency_step_Hz = 8 * kHzcenter_frequency_Hz = backend_defaults.meas_freq_est[qubit]print(f"Qubit {qubit} has an estimated readout frequency of {center_frequency_Hz / GHz} GHz.")frequency_min = center_frequency_Hz - frequency_span_Hz / 2frequency_max = center_frequency_Hz + frequency_span_Hz / 2frequencies_GHz = np.arange(frequency_min / GHz,                             frequency_max / GHz,                             frequency_step_Hz / GHz)print(f"The sweep will go from {frequency_min / GHz} GHz to {frequency_max / GHz} GHz\      in steps of {frequency_step_Hz / MHz} MHz.")

Out:

Qubit 0 has an estimated readout frequency of 6.993427855 GHz.The sweep will go from 6.993267855 GHz to 6.993587855 GHz      in steps of 0.008 MHz.

# Qiskit Pulse — Dispersive Shift

num_shots_per_frequency = 2048frequencies_Hz = frequencies_GHz*GHzschedule_los = [{meas_chan: freq} for freq in frequencies_Hz]cavity_sweep_0 = assemble(schedule_0,                               backend=backend,                                meas_level=1,                               meas_return='avg',                               shots=num_shots_per_frequency,                               schedule_los=schedule_los)cavity_sweep_1 = assemble(schedule_1,                               backend=backend,                                meas_level=1,                               meas_return='avg',                               shots=num_shots_per_frequency,                               schedule_los=schedule_los)# RUN the job on a real devicejob_0 = backend.run(cavity_sweep_0)job_monitor(job_0)job_0.error_message()job_1 = backend.run(cavity_sweep_1)job_monitor(job_1)job_1.error_message()# OR retreive result from previous run#job_0 = backend.retrieve_job('Job ID1')#job_1 = backend.retrieve_job('Job ID2')cavity_sweep_0_results = job_0.result()cavity_sweep_1_results = job_1.result()

Out:

Job Status: job has successfully runJob Status: job has successfully run

# Qiskit Pulse — Dispersive Shift

scale_factor = 1e-14sweep_values_0 = []for i in range(len(cavity_sweep_0_results.results)):    res_0 = cavity_sweep_0_results.get_memory(i)*scale_factor    sweep_values_0.append(res_0[qubit])sweep_values_1 = []for i in range(len(cavity_sweep_1_results.results)):    res_1 = cavity_sweep_1_results.get_memory(i)*scale_factor    sweep_values_1.append(res_1[qubit])    plotx = frequencies_Hz/kHzploty_0 = np.abs(sweep_values_0)ploty_1 = np.abs(sweep_values_1)

# Qiskit Pulse — Dispersive Shift

# plot real part of sweep valuesplt.plot(plotx, ploty_0, color='blue', marker='.')# plot real part of sweep values plt.plot(plotx, ploty_1, color='red', marker='.') plt.legend([r'$\vert0\rangle$', r'$\vert1\rangle$'])plt.grid()plt.xlabel("Frequency [kHz]")plt.ylabel("Measured signal [a.u.]")plt.yscale('log')plt.show()

References:

Qiskit (Global Summer School), Introduction to Quantum Computing and Quantum Hardware — Lab 7.