Qiskit Dynamicsを用いて２量子ビット操作シミュレーション

いよいよ２量子ビット操作をシミュレートしていきます．今回はSimulating backends at the pulse-level with DynamicsBackendを参考にしています．

2量子ビット操作について

IBMの量子コンピュータは，周波数の異なるトランズモン同士を結合させ，相方の周波数に共鳴するようなマイクロ波を片方のトランズモンに当てることでCRゲートを実現させています．このようにマイクロ波のみで２量子ビットゲートを操作する手法は以前から提案されていましたが，実験結果を正確に議論する理論というのは，意外にも最近論文として出されたようです．ここではとりあえず「周波数の異なるトランズモン同士を結合させ，相方の周波数に共鳴するようなマイクロ波を片方のトランズモンに当てることで相方の状態操作ができている」ことを確認してみます．

Qiskit Dynamics

量子ビット系のハミルトニアンを手で与えることにより，パルスレベルの操作をシミュレートするモジュールのようです．

まず，効率的にシミュレートするために，JAXを用いるよう設定します．

import jax
jax.config.update("jax_enable_x64", True)
jax.config.update("jax_platform_name", "cpu")
from qiskit_dynamics.array import Array
Array.set_default_backend("jax")

次にハミルトニアンを設定します．
$$
\begin{align}
H(t)&=2 \pi \nu_0 N_0+\pi \alpha_0 N_0\left(N_0-I\right)+2 \pi \nu_1 N_1+\pi \alpha_1 N_1\left(N_1-I\right)+2 \pi J\left(a_0+a_0^{\dagger}\right)\left(a_1+a_1^{\dagger}\right) \\
&\quad +2 \pi r_0 s_0(t)\left(a_0+a_0^{\dagger}\right)\\
N_0&=a_0^{\dagger}a_0,\quad N_1=a_1^{\dagger}a_1
\end{align}
$$
それぞれの量子ビットのエネルギーは$2\pi \nu_i$で与えられ，$\alpha_i$は非線形性を表しています．

import numpy as np

dim = 3

v0 = 4.86e9
anharm0 = -0.32e9
r0 = 0.22e9

v1 = 4.97e9
anharm1 = -0.32e9
r1 = 0.26e9

J = 0.002e9

a = np.diag(np.sqrt(np.arange(1, dim)), 1)
adag = np.diag(np.sqrt(np.arange(1, dim)), -1)
N = np.diag(np.arange(dim))

ident = np.eye(dim, dtype=complex)
full_ident = np.eye(dim**2, dtype=complex)

N0 = np.kron(ident, N)
N1 = np.kron(N, ident)

a0 = np.kron(ident, a)
a1 = np.kron(a, ident)

a0dag = np.kron(ident, adag)
a1dag = np.kron(adag, ident)


static_ham0 = 2 * np.pi * v0 * N0 + np.pi * anharm0 * N0 * (N0 - full_ident)
static_ham1 = 2 * np.pi * v1 * N1 + np.pi * anharm1 * N1 * (N1 - full_ident)

static_ham_full = static_ham0 + static_ham1 + 2 * np.pi * J * ((a0 + a0dag) @ (a1 + a1dag))

drive_op0 = 2 * np.pi * r0 * (a0 + a0dag)
drive_op1 = 2 * np.pi * r1 * (a1 + a1dag)

ここまでの設定で，以下のようにbackendを設定できます．

from qiskit_dynamics import Solver
from qiskit_dynamics import DynamicsBackend

dt = 1/4.5e9

solver = Solver(
    static_hamiltonian=static_ham_full,
    hamiltonian_operators=[drive_op0, drive_op1, drive_op0, drive_op1],
    rotating_frame=static_ham_full,
    hamiltonian_channels=["d0", "d1", "u0", "u1"],
    channel_carrier_freqs={"d0": v0, "d1": v1, "u0": v1, "u1": v0},
    dt=dt,
)


solver_options = {"method": "jax_odeint", "atol": 1e-6, "rtol": 1e-8, "hmax": dt}

backend = DynamicsBackend(
    solver=solver,
    subsystem_dims=[dim, dim],
    solver_options=solver_options
)

さらに，これまで同様以下の設定を行っておきます．

GHz = 1.0e9 # Gigahertz
MHz = 1.0e6 # Megahertz
us = 1.0e-6 # Microseconds
ns = 1.0e-9 # Nanoseconds

# 以下の量子ビットについて量子ビット周波数を探索する
qubit = 0

# ガウシアンパルスのパラメタ設定(us = microseconds)
drive_sigma_sec = 0.3 * us    # ガウシアンの幅を決定，大きめに取った
drive_duration_sec = drive_sigma_sec * 8  # 切り捨てパラメータ。ガウシアンを有限の長さで切り捨てる
drive_amp = 0.05

def frequencies_series(center_frequency_Hz):

    # 推定周波数を中心に100MHzのスパンでスイープする
    frequency_span_Hz = 100 * MHz
    # 1MHzの単位で
    frequency_step_Hz = 1 * MHz

    # 推定周波数より20 MHz高く、20 MHz低くスイープ
    frequency_min = center_frequency_Hz - frequency_span_Hz / 2
    frequency_max = center_frequency_Hz + frequency_span_Hz / 2
    # 実験のために周波数のnp配列をGHz単位で作成
    frequencies_GHz = np.arange(frequency_min / GHz,
                                frequency_max / GHz,
                                frequency_step_Hz / GHz)
    
    return frequencies_GHz

今回は$0$番目の量子ビットにマイクロ波を照射することを考えます．

# デフォルトの推定周波数を中心値としてスイープを行う
center_frequency_Hz = v1
print(f"量子ビット {qubit} の推定周波数は {center_frequency_Hz / GHz} GHz")

# 実験のために周波数のnp配列をGHz単位で作成
frequencies_GHz = frequencies_series(center_frequency_Hz)

print(f"スイープは {frequencies_GHz[0]} GHz から {frequencies_GHz[-1]} GHz まで、行われる")

import qiskit.pulse as pulse
from qiskit.circuit import Parameter
# ベーススケジュールを作成
# ドライブパルスがドライブチャネルで実行されるように開始する
freq = Parameter('freq')
with pulse.build(backend=backend, default_alignment='sequential', name='Frequency sweep') as sweep_sched:
    drive_duration = pulse.seconds_to_samples(drive_duration_sec)
    drive_sigma = pulse.seconds_to_samples(drive_sigma_sec)
    drive_chan = pulse.DriveChannel(qubit)

    pulse.set_frequency(freq, drive_chan)
    
    pulse.play(pulse.Gaussian(duration=drive_duration,
                              sigma=drive_sigma,
                              amp=drive_amp,
                              name='freq_sweep_excitation_pulse'), drive_chan)

この設定したpulseを，実際に回路に載せます．

from qiskit import pulse
from qiskit.circuit import QuantumCircuit, Gate

sweep_gate = Gate("sweep", 1, [freq])

qc_sweep = QuantumCircuit(2, 2)
qc_sweep.append(sweep_gate, [qubit])
qc_sweep.measure([0, 1], [0, 1])

qc_sweep.add_calibration(gate = sweep_gate, qubits = [qubit], schedule = sweep_sched, params = [freq])

frequencies_Hz = frequencies_GHz*GHz
exp_sweep_circs = [qc_sweep.assign_parameters({freq: f}, inplace=False) for f in frequencies_Hz]

さて，自作したbackendで実行してみましょう．

num_shots_per_frequency = 1024

job = backend.run(exp_sweep_circs,
                  meas_level=1,
                  meas_return='avg',
                  shots=num_shots_per_frequency)

この結果を見てみます．まずは$0$番目の量子ビットについて．

import matplotlib.pyplot as plt

frequency_sweep_results = job.result() 

sweep_values = []
for i in range(len(frequency_sweep_results.results)):
    res = frequency_sweep_results.get_memory(i)
    sweep_values.append(res[0])

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

$\nu_1$周りでスイープしているため，$0$番目の量子ビットには何も起きていないことが分かります．一方，res[1]として$1$番目の量子ビットをプロットすると，

周波数が$\nu_1$となる箇所で異なる測定結果を得ていることが分かります．

最後に

次は理想的なCRゲートとどの程度異なったゲートになっているのか，また，それを修正する手法をシミュレートしてみたいと思います．