The new Qiskit Textbook beta is now available. Try it out now
Measuring the Qubit ac-Stark Shift

# Measuring the Qubit ac-Stark Shift

### Physics Background

Let's consider a qubit with frequency $\omega_q$ strongly coupled to a resonator with frequency $\omega_r$ with $\omega_q<\omega_r$; the qubit-resonator coupling strength is $g$ and the detuning is $\Delta=\omega_q-\omega_r$. In the dispersive limit, the system can be described using the following Hamiltonian:

$H_{JC(disp)}/\hbar=\omega_r (a^\dagger a+\frac{1}{2}) + \frac{1}{2} (\omega_q + \frac{g^2}{\Delta} + \frac{2g^2}{\Delta} a^\dagger a) \sigma_z$

where $a$ and $a^\dagger$ are the raising a lowering operators of the resonator photons, and $\sigma_z$ is the Pauli-Z operator acting on the qubit. In this frame, the qubit frequency

$\tilde{\omega}_q=\omega_q + \frac{g^2}{\Delta} + \frac{2g^2}{\Delta} \bar{n}$

experiences a constant Lamb shift of $g^2/\Delta$ induced by the vacuum fluctuations in the resonator, and an ac-Stark shift of $(2g^2/\Delta) \bar{n}$ where $\bar{n}=\langle a^\dagger a \rangle$ is the average number of photons present in the resonator. For more details checkout this paper. In this tutorial, we investigate the ac-Stark shift of the qubit caused by the photon population in the resonator using Qiskit Pulse.

### 0. Getting started

We'll first set up our basic dependencies so we're ready to go.

# Importing standard Qiskit libraries and configuring account
from qiskit import QuantumCircuit, execute, Aer, IBMQ
from qiskit.compiler import transpile, assemble
provider = IBMQ.get_provider(hub='ibm-q')
backend = provider.get_backend('ibmq_armonk')


We then extract the default backend configuration and settings for the selected chip.

backend_config = backend.configuration()
backend_defaults = backend.defaults()
dt=backend_config.dt  # hardware resolution
backend.configuration().parametric_pulses = [] # will allow us to send a larger waveform for our experiments


Next, we define some helper functions that we will use for fitting and interpreting our data.

from scipy.optimize import leastsq,minimize, curve_fit

# samples need to be multiples of 16 to accommodate the hardware limitations
def get_closest_multiple_of_16(num):
return int(num + 8 ) - (int(num + 8 ) % 16)

# lorentzian function
def lorentzian(f, f0, k, a, offs):
return a*k/(2*np.pi)/((k/2)**2+(f-f0)**2)+offs

#fit_lorentzian takes two arrays that contain the frequencies and experimental output values of each frequency respectively.
#returns the lorentzian parameters that best fits this output of the experiment.
#popt are the fit parameters and pcov is the covariance matrix for the fit
def fit_lorentzian(freqs,values):
p0=[freqs[np.argmax(values)],(freqs[-1]-freqs[0])/2,max(values),0]
bounds=([freqs[0],0,0,-np.inf],[freqs[-1],freqs[-1]-freqs[0],np.inf,np.inf])
popt,pcov=curve_fit(lorentzian, freqs, values, p0=p0, bounds=bounds)
return popt,pcov

# Gaussian function
def gaussian(f, f0, sigma, a, offs):
return a*np.exp(-(f-f0)**2/(2*sigma**2))+offs

#fit_gaussian takes two arrays that contain the frequencies and experimental output values of each frequency respectively.
#returns the gaussian parameters that best fits this output of the experiment.
#popt are the fit parameters and pcov is the covariance matrix for the fit
def fit_gaussian(freqs,values):
p0=[freqs[np.argmax(values)],(freqs[-1]-freqs[0])/2,max(values),0]
bounds=([freqs[0],0,0,-np.inf],[freqs[-1],freqs[-1]-freqs[0],np.inf,np.inf])
popt,pcov=curve_fit(gaussian, freqs, values, p0=p0, bounds=bounds)
return popt,pcov

# normalize the data points to fall in the range of [0,1]
def normalize(a):
a= a-min(a)
return a/max(a)


### 1. ac-Stark Shifting the qubit

In order to ac-Stark shift the qubit we need to populate the resonator with photons using an on-resonance drive. For a drive amplitude $\epsilon$, and a resonator decay rate of $\kappa$, the number of photons in the resonator $\bar{n}=\langle a^\dagger a \rangle = \frac{\epsilon^2}{\Delta^2 +(\kappa/2)^2}$. As a reminder $\tilde{\omega}_q=\omega_q + \frac{g^2}{\Delta} + \delta \omega_q$ where the shift in frequency due to ac-Stark shift is $\delta \omega_q = \frac{2g^2}{\Delta} \bar{n}$. Since $\Delta=\omega_q-\omega_r<0$ the qubit frequency gets smaller as we increase the of photons in the resonator.

from qiskit import pulse            # This is where we access all of our Pulse features!
from qiskit.pulse import Play, Acquire
import qiskit.pulse.library as pulse_lib
import numpy as np

qubit=0   # qubit used in our experiment

inst_sched_map = backend_defaults.instruction_schedule_map
# Get the default measurement pulse sequence
measure = inst_sched_map.get('measure', qubits=backend_config.meas_map[0])

qubit_drive_sigma = 100e-9           #the width of the qubit spectroscopy drive
resonator_drive_sigma=10e-9          #the width of the resonator drive
drive_duration=10*qubit_drive_sigma  #the resonator drive duration

# We use a Gaussian shape pulse to drive the qubit for spectroscopy
qubit_drive = pulse_lib.gaussian(duration = get_closest_multiple_of_16(drive_duration//dt),
amp = .1,
sigma = get_closest_multiple_of_16(qubit_drive_sigma//dt),
name = 'qubit tone')

drive_chan = pulse.DriveChannel(qubit)   # qubit drive channel
meas_chan = pulse.MeasureChannel(qubit)  # resonator channel
acq_chan = pulse.AcquireChannel(qubit)   # readout signal acquisition channel

measurement_delay=20e-9*3

resonator_tone_amplitude = np.linspace(0,1,11) #change to amplitude
resonator_tone_pulses = []
for amp in resonator_tone_amplitude:
# we use a square pulse with Gaussian rise and fall time to populate the resonator with photons
temp_resonator_tone=pulse_lib.GaussianSquare(duration = get_closest_multiple_of_16(drive_duration//dt),
amp = amp,
sigma = get_closest_multiple_of_16(resonator_drive_sigma//dt),
width = get_closest_multiple_of_16((drive_duration-4*resonator_drive_sigma)//dt),
name = 'resonator tone')
# pulse sequence for the experiment at different amplitudes
with pulse.build(name=f"resonator tone amplitude = {np.round(amp,2)} V") as temp_pulse:
pulse.play(qubit_drive, drive_chan)
pulse.play(temp_resonator_tone, meas_chan)
pulse.delay(int(measurement_delay//dt), meas_chan)
pulse.call(measure)

resonator_tone_pulses.append(temp_pulse)

resonator_tone_pulses[2].draw()

start=4.960e9  # qubit spectroscopy start frequency
stop=4.980e9   # qubit spectroscopy stop frequency
freqs = np.linspace(start, stop, 41)-500e3
# list of qubit drive frequencies for spectroscopy
schedule_frequencies = [{drive_chan: freq , meas_chan: backend_defaults.meas_freq_est[qubit]} for freq in freqs]


Here, we send our pulse sequence to the hardware.

from qiskit.tools.monitor import job_monitor

num_shots = 4*1024

resonator_tone_experiments=[(resonator_tone_pulses[i],
{'meas_level':  1,
'meas_return': 'avg',
'shots':        num_shots,
'schedule_los': schedule_frequencies}) for i in range(len(resonator_tone_pulses))]

resonator_tone_results=[]
for experiment, args in resonator_tone_experiments:
job = backend.run(experiment, **args)
job_monitor(job)
resonator_tone_results.append(job.result(timeout=120))

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


And then we access the measurement data.

import matplotlib.pyplot as plt

skip_idx=0

resonator_tone_values = []
for result in resonator_tone_results:
result_values=[]
for i in range(len(result.results)):
result_values.append(result.get_memory(i)[qubit])

resonator_tone_values.append(normalize(result_values))

plot_extent=[freqs[0]/1e9,freqs[-1]/1e9,resonator_tone_amplitude[skip_idx],resonator_tone_amplitude[-1]]
plt.imshow(np.abs(resonator_tone_values[skip_idx:]), aspect='auto', origin='lower', cmap='viridis',
extent=plot_extent)

plt.xlabel('Qubit tone frequency [GHz]')
plt.ylabel('Resonator tone amplitude [V]')
plt.title('Qubit ac-Stark shift')
plt.show()


### 2. Qubit frequency shift and linewidth broadening

Using the Jaynes-Cummings model we expect a qubit frequency shift of $\delta \omega_q = \frac{2g^2}{\Delta} \bar{n}$. The qubit frequency experiences fluctuations due the photon shot-noise which leads to qubit linewidth broadening and a dephasing rate of $\Gamma_\phi=\frac{4 \chi^2}{\kappa} \bar{n}$.

show_individual_traces=False
skip_idx=3  # number of points to skip

center=[]
fwhm=[]
for i in range(len(resonator_tone_values)):
popt,pcov=fit_gaussian(freqs,np.abs(np.real(resonator_tone_values[i])))
center.append(popt[0])
fwhm.append(2.355*popt[1])
if show_individual_traces:
plt.plot(freqs/1e3, np.real(resonator_tone_values[i]))
plt.plot(freqs/1e3, gaussian(freqs,*popt), '--')
if show_individual_traces: plt.show()

center_fit=np.polyfit(resonator_tone_amplitude[skip_idx:], (center[skip_idx:]-center[0]),1)
plt.plot(resonator_tone_amplitude[skip_idx:], np.poly1d(center_fit/1e6)(resonator_tone_amplitude[skip_idx:]), '--', lw=2, color='grey')
plt.plot(resonator_tone_amplitude[skip_idx:], (center[skip_idx:]-center[0])/1e6, 'o', color='black')
plt.xlabel(r'Resonator tone amplitude [V]')
plt.ylabel(r'$\delta \omega_q (MHz)$')
plt.show()

fwhm_fit=np.polyfit(resonator_tone_amplitude[skip_idx:], np.array(fwhm[skip_idx:]),1)
plt.plot(resonator_tone_amplitude[skip_idx:], np.poly1d(fwhm_fit/1e6)(resonator_tone_amplitude[skip_idx:]), '--', lw=2, color='orange')
plt.plot(resonator_tone_amplitude[skip_idx:], np.array(fwhm[skip_idx:])/1e6, 'o', color='red')
plt.xlabel(r'Resonator tone amplitude [V]')
plt.ylabel(r'FWHM (MHz)')
plt.show()


In this chapter, we discuss the ac-Stark shift that the qubit experiences due to the presence of photons in the resonator. We use Qiskit Pulse to measure the qubit frequency shift and linewidth broadening.

import qiskit.tools.jupyter
%qiskit_version_table


### Version Information

Qiskit SoftwareVersion
Qiskit0.27.0
Terra0.17.4
Aer0.8.2
Ignis0.6.0
Aqua0.9.2
IBM Q Provider0.14.0
System information
Python3.7.7 (default, May 6 2020, 04:59:01) [Clang 4.0.1 (tags/RELEASE_401/final)]
OSDarwin
CPUs8
Memory (Gb)32.0
Thu Jun 17 10:50:12 2021 BST