はじめに

2020年7月に開催されたQiskit Global Summer Schoolに参加してきました。その中のセクションでパルス制御のセッションがあり、そこでQiskitのOpenPulseを知ったので、今回はOpenPulseで遊んでみたいと思います。

はじめに
参考文献
2準位系の量子状態の操作
共鳴周波数を測定してみる
まとめ

$\def\bra#1{\mathinner{\left\langle{#1}\right|}} \def\ket#1{\mathinner{\left|{#1}\right\rangle}} \def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}$

参考文献

このQiskit公式ノートを参考にします。
https://qiskit.org/textbook/ch-quantum-hardware/calibrating-qubits-openpulse.htmlqiskit.org

2準位系の量子状態の操作

f:id:SJSY:20200803214300p:plain

量子ビットの2準位 $\ket{0}, \ket{1}$ の量子状態を操作するには、2準位間のエネルギー差 $\Delta E$ に共鳴する周波数の電磁波を照射します。一般的に共鳴周波数はGHzのオーダのため、マイクロ波を照射することになります。
この共鳴するマイクロ波のパルス列を用いることで、ゲート方式の量子コンピュータで用いられる $X$ ゲートや $H$ ゲートなどのゲート操作は実現されています。詳しくは核磁気共鳴（NMR）あたりを勉強するといいかもしれません。

では早速、QiskitのOpenPulseを用いてIBMの実機で量子状態の操作をしてみましょう。

共鳴周波数を測定してみる

IBMの実機を使うので、IBM Quantum Experienceのアカウント登録が必要です。初回実行時にAPI Tokenを登録します。API TokenはIBM Quantum ExperienceのMy Accountから入手できます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

from qiskit.tools.jupyter import *
from qiskit import IBMQ
from qiskit import pulse
from qiskit.pulse import Play
from qiskit.pulse import pulse_lib
from qiskit import assemble
from qiskit.tools.monitor import job_monitor

#単位の設定
GHz = 1.0e9 # Gigahertz
MHz = 1.0e6 # Megahertz
us = 1.0e-6 # Microseconds
ns = 1.0e-9 # Nanoseconds
scale_factor = 1e-14

IBMQ.save_account(<API_TOKEN>) # 初回
#IBMQ.load_account() # 2回目以降

OpenPulseを実行できるバックエンド（実機）は決まっているので、適切なバックエンドを選択してあげます。2020年8月時点では、ibmq_armonkのみ。ibmq_armonkは1量子ビットの量子コンピュータですね。

IBMQ.load_account()
provider = IBMQ.get_provider(hub='ibm-q', group='open', project='main')
backend = provider.get_backend('ibmq_armonk')

# 設定したバックエンドがOpenPulseが使えるか確認
backend_config = backend.configuration()
assert backend_config.open_pulse, "Backend doesn't support Pulse"

バックエンドには、パルスのサンプリング時間間隔 $dt$ が設定されており、パルス制御において重要なパラメタになってくるので確認しておきます。

dt = backend_config.dt
print(f"Sampling time: {dt*1e9} ns")

Sampling time: 0.2222222222222222 ns

backend.defaultsには量子ビットの共鳴周波数の推定値が格納されているので確認しておきます。

backend_defaults = backend.defaults()
print(backend_defaults)

<PulseDefaults(<InstructionScheduleMap(1Q instructions:
  q0: {'measure', 'u1', 'x', 'u3', 'u2', 'id'}
Multi qubit instructions:
)>Qubit Frequencies [GHz]
[4.974450544012941]
Measurement Frequencies [GHz]
[6.993427855] )>

4.974GHzあたりに共鳴周波数があることがわかりますね。実機の共鳴周波数は量子ビットの物性や量子ビットの周りの環境によって微妙に変わってくるので、測定して確かめることが大事です。

ここでは、4.974GHzの前後20MHz周波数を掃引して共鳴周波数を確認してみましょう。共鳴する周波数のときには2準位間で遷移が起きて大きく量子状態が変化するので、測定値の変化する周波数を見つけます。

掃引する周波数の設定をしておきます。

# 操作する量子ビットの設定。1量子ビットなので0。
qubit = 0

# 掃引する周波数の中心周波数
center_frequency_Hz = backend_defaults.qubit_freq_est[qubit]
print(f"Qubit {qubit} has an estimated frequency of {center_frequency_Hz / GHz} GHz.")

# 掃引する周波数のレンジとステップを設定。ステップ数多いと実行時にエラー出る。
frequency_span_Hz = 40 * MHz
frequency_step_Hz = 1 * MHz

frequency_min = center_frequency_Hz - frequency_span_Hz / 2
frequency_max = center_frequency_Hz + frequency_span_Hz / 2

frequencies_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.")

Qubit 0 has an estimated frequency of 4.974450544012941 GHz.
The sweep will go from 4.954450544012941 GHz to 4.994450544012941 GHz in steps of 1.0 MHz.

続いて量子ビットの2準位を遷移(drive)させるマイクロ波パルスの設定をします。

#サンプル数は16の倍数で設定します。
def get_closest_multiple_of_16(num):
    return int(num + 8 ) - (int(num + 8 ) % 16)

#ガウシアンパルスのパラメタ設定
drive_sigma_us = 0.075
drive_samples_us = drive_sigma_us*8

drive_sigma = get_closest_multiple_of_16(drive_sigma_us * us /dt)
drive_samples = get_closest_multiple_of_16(drive_samples_us * us /dt) 

drive_amp = 0.3

#パルス生成
drive_pulse = pulse_lib.gaussian(duration=drive_samples,
                                 sigma=drive_sigma,
                                 amp=drive_amp,
                                 name='freq_sweep_excitation_pulse')

続いて読み出し（測定）パルスの設定です。実機では測定マップを確認しておく必要があります。ハードウエアの制約によって、1つの量子ビットに対して読み出しが実行されると他の量子ビットに対しても読み出しが行われてしまいます。測定する量子ビットがどのグループにあるかをチェックしてみます。

meas_map_idx = None
for i, measure_group in enumerate(backend_config.meas_map):
    if qubit in measure_group:
        meas_map_idx = i
        break
assert meas_map_idx is not None, f"Couldn't find qubit {qubit} in the meas_map!"

読み出しパルスはバックエンドで校正されているので、デフォルトのままが正確ぽいです。

inst_sched_map = backend_defaults.instruction_schedule_map
measure = inst_sched_map.get('measure', qubits=backend_config.meas_map[meas_map_idx])

最後にパルスを照射するチャンネルを設定しておきます。

drive_chan = pulse.DriveChannel(qubit)
meas_chan = pulse.MeasureChannel(qubit)
acq_chan = pulse.AcquireChannel(qubit)

ここまでで、照射する各パルスの形状の設定ができました。続いて、パルスを照射するタイミングと順番、つまりパルス列の設定をしてきます。

schedule = pulse.Schedule(name='Frequency sweep')
schedule += Play(drive_pulse, drive_chan)
schedule += measure << schedule.duration
frequencies_Hz = frequencies_GHz*GHz
schedule_frequencies = [{drive_chan: freq} for freq in frequencies_Hz]

設定したパルス列を確認しておきましょう。

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

f:id:SJSY:20200804211942p:plain

さあ、いよいよ実機にジョブを追加して実行してみましょう。 meas_level=1でショットごとに1つの複素数値を返し、meas_return='avg'で平均値を求めます。

num_shots_per_frequency = 8192
frequency_sweep_program = assemble(schedule,
                                   backend=backend, 
                                   meas_level=1,
                                   meas_return='avg',
                                   shots=num_shots_per_frequency,
                                   schedule_los=schedule_frequencies)

job = backend.run(frequency_sweep_program)

#print(job.job_id())
#ジョブの実行状況を監視
job_monitor(job)

Job Status: job has successfully run

無事、実機で実行できたみたいです。
結果を確認してみましょう。

frequency_sweep_results = job.result(timeout=120) 
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='black') # plot real part of sweep values
plt.xlim([min(frequencies_GHz), max(frequencies_GHz)])
plt.xlabel("Frequency [GHz]")
plt.ylabel("Measured signal [a.u.]")
plt.show()

f:id:SJSY:20200804215213p:plain

4.975 MHzあたりにピークがあり、このピークが共鳴周波数に対応しています。このピークの形状は、ドライブパルスのガウシアンのパラメタを変えると変わります。

たとえば、drive_sigma_usを0.075から0.05に変えて、実行してみます。

f:id:SJSY:20200804215226p:plain

ガウスアンが少しシャープになりました。実行してみます。

f:id:SJSY:20200804215759p:plain

ピークの形が変わりました。ちなみに、ガウスアンのパワーが強ければ良いという訳ではありません。いろいろパラメタを変えて確認してみましょう。理由は次のラビ振動を知ると理解できます。

測定した結果をフィッティングして、共鳴周波数に対応するピーク周波数を求めましょう。共振のフィッティングに一般的に用いられているローレンツ関数でフィッティングをかけます。フィッティングの初期値はピーク波形をみてそれっぽい値を入れておきます。(でないと、エラーでフィッティングできません)

def 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_fit

fit_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, 0.2, 1] #ピーク波形をみてそれっぽい値を入れておく                                )

plt.scatter(frequencies_GHz, np.real(sweep_values), color='black')
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()

f:id:SJSY:20200804220530p:plain

A, rough_qubit_frequency, B, C = fit_params
rough_qubit_frequency = rough_qubit_frequency*GHz 
print(f"We've updated our qubit frequency estimate from "
      f"{round(backend_defaults.qubit_freq_est[qubit] / GHz, 5)} GHz to {round(rough_qubit_frequency/GHz, 5)} GHz.")

We've updated our qubit frequency estimate from 4.97445 GHz to 4.97453 GHz.

共鳴周波数を測定した結果、この量子ビットの共鳴周波数は4.97453 GHzだと分かりました。