Qiskit pulseでRabi振動

[Qiskit Pulseによる量子ビットの較正](https://qiskit.org/textbook/ja/ch-quantum-hardware/calibrating-qubits-pulse.html#calib)に従って，Rabi振動を見てみることにします．

# Rabi振動
これまで同様，次のハミルトニアンを考えます．

$$
\begin{align}
    H(t)=-\frac{1}{2} \hbar (\omega_q-\omega_d) \sigma^z-\hbar \Omega(t) \sigma^x
\end{align}
$$

マイクロ波のエンベロープ関数はガウシアンとしましょう[^2]．

$$
\begin{align}
    \hbar\Omega(t)&\propto
\begin{cases}
A  \frac{f(t)-f(-1)}{1-f(-1)}& 0\leq t<T \\\\
0& \text { others }
\end{cases}\\\\
f(t)&= \exp \left(-\frac{1}{2} \frac{(t-\frac{T}{2})^2}{\sigma^2}\right)
\end{align}
$$

さて，共鳴周波数$\omega_d=\omega_q$では，時間発展演算子は

$$
\begin{align}
    U(t,0)=e^{i\int_0^t \Omega(t')dt' \sigma^x}
\end{align}
$$

と表せます．これより，ガウシアンパルスの面積によって，$x$軸回転の回転角が変わる(=Rabi振動)ことが分かります．ガウシアンパルスの面積を変えるために，様々な振幅$A$を用いて量子ビットを操作してみます．

# 実験

コードは[前回](https://qiita.com/tare_/items/ccef2856903eef886cf6)とほぼ同じです．

```python
from qiskit import IBMQ
import numpy as np
from qiskit import pulse                  
from qiskit.circuit import Parameter      
from qiskit.circuit import QuantumCircuit, Gate


#IBMQ.save_account("") # 最初はtokenの入力が必要
IBMQ.load_account()
provider = IBMQ.get_provider(hub='ibm-q')
backend = provider.get_backend('') # 使えるbackendを選んでください

backend_config = backend.configuration()
backend_defaults = backend.defaults()

# 単位変換係数（バックエンドプロパティはHz、秒で返される）
GHz = 1.0e9 # Gigahertz
MHz = 1.0e6 # Megahertz
us = 1.0e-6 # Microseconds
ns = 1.0e-9 # Nanoseconds

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

# value を base_number の倍数に丸める関数
def get_closest_multiple_of(vaule, base_number):
    return int(vaule + base_number/2) - (int(vaule + base_number/2) % base_number)

# num を granularity の倍数に丸める関数
def get_closest_multiple_of_granu(num):
    return get_closest_multiple_of(num, granularity)

granularity = backend.configuration().timing_constraints['granularity']
print(granularity)
# 8

# デフォルトの推定周波数をそのまま用いる
center_frequency_Hz = backend_defaults.qubit_freq_est[qubit]  

# ガウシアンパルスのパラメタ設定(us = microseconds)
drive_sigma_sec = 15 * ns   # ガウシアンの幅を決定
drive_duration_sec = 120 * ns   # パルスの全時間幅を決定
```
ここまではほぼ同じです．次に振幅の値を変えて実行するコードです[^1]．
```python
# 駆動パルスの振幅の種類
num_rabi_points = 50

# 駆動パルスの振幅値：0から0.75の間で等間隔に配置された50個の振幅
drive_amp_min = 0
drive_amp_max = 0.75
drive_amps = np.linspace(drive_amp_min, drive_amp_max, num_rabi_points)

# ベーススケジュールを作成
amp = Parameter('amp')
with pulse.build(backend=backend, default_alignment='sequential', name='Amplitude sweep') as sweep_sched:
    drive_duration = get_closest_multiple_of_granu(pulse.seconds_to_samples(drive_duration_sec))
    drive_sigma = pulse.seconds_to_samples(drive_sigma_sec)
    drive_chan = pulse.DriveChannel(qubit)
    pulse.set_frequency(center_frequency_Hz, drive_chan)
    # Drive pulse samples
    pulse.play(pulse.Gaussian(duration=drive_duration,
                              sigma=drive_sigma,
                              amp=amp,
                              name='amp_sweep_excitation_pulse'), drive_chan)
```
これ以降は，前回同様ゲートを作って，作用させて実行していきます．
```python
sweep_gate = Gate("amp", 1, [amp])
qc_sweep = QuantumCircuit(1, 1)
qc_sweep.append(sweep_gate, [0])
qc_sweep.measure(0, 0)
qc_sweep.add_calibration(sweep_gate, (0,), sweep_sched, [amp])

exp_sweep_circs = [qc_sweep.assign_parameters({amp: drive_amp}, inplace=False) for drive_amp in drive_amps]

# 実験
num_shots_per_amp = 1024

job = backend.run(exp_sweep_circs,
                  meas_level=1,
                  meas_return='avg',
                  shots=num_shots_per_amp)
```
この結果，次のような結果が得られました．
![osc.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/887541/6ba98116-d4ac-80c9-c695-a5fc522238b6.png)

確かに振幅を変更することによって，$x$軸周りの回転角が変わっていることが見て取れます．

さて，$|f(-1)|\ll1$とし，さらに$\hbar \Omega(t)=xAf(t)$として比例係数$x$を導入すると，時間発展演算子は
$$
\begin{align}
    U(t,0)&=e^{i\frac{xA}{\hbar}\int_0^t f(t')dt' \sigma^x}\\\\
    &=\cos\left(\frac{xA}{\hbar}\int_0^t f(t')dt'\right)I+i\sin\left(\frac{xA}{\hbar}\int_0^t f(t')dt'\right)\sigma_x
\end{align}
$$

と表されます．つまり，
$$
\begin{align}
    \frac{xA}{\hbar}\int_0^t f(t')dt'=\pi
\end{align}
$$

で，状態が（全体位相を除いて）元に戻ることが分かります．今回用いたパラメータを代入すると，
$$
\begin{align}
    xA=\frac{\hbar\pi}{37.6\mathrm{ns}}&\simeq 0.28x\\\\
    x&\simeq 1.96\times 10^{-7}\mathrm{eV}\simeq 0.3\times \frac{\hbar}{\mathrm{ns}}
\end{align}
$$
と求まります．これより，``drive_amp=1``と設定すると，ハミルトニアンに含まれる$\Omega(t)$の最大値は$x$となることが分かりました．

# 参考文献 
- [QiskitのOpenPulseでパルス制御【ラビ振動】](https://sjsy.hatenablog.com/entry/20200823/1598177520)
- [Qiskit Pulseによる量子ビットの較正](https://qiskit.org/textbook/ja/ch-quantum-hardware/calibrating-qubits-pulse.html#calib)


[^2]:コードで指定するガウシアンパルスの振幅を$A$と表しています．一方，[トランズモン物理入門](https://qiskit.org/textbook/ja/ch-quantum-hardware/transmon-physics.html)を参照すると，$\Omega(t)$の出自は双極子相互作用となっており，ガウシアンパルスの振幅がそのまま$\Omega(t)$と表せるわけではありません．

[^1]:ここのコードの書き方が合っているのかどうかは未検証です．