0. これは何の記事？
- ゲート型の量子コンピュータを使ったシミュレーションをやってみた
- シミュレーションの対象はニュートリノ振動と呼ばれる物理現象
- 下記の文献を参考に２フレーバーの場合をシミュレーションしてみた
Neutrino oscillations in a quantum processor
https://journals.aps.org/prresearch/pdf/10.1103/PhysRevResearch.1.033176
参考文献の内容をかいつまんで説明しながら、私がやった結果を掲載します。
なお今回作成したコードはこちらのGitHubにあります。
1. ニュートリノ振動
ニュートリノは素粒子で、これまでに３種類が知られています。それぞれを電子ニュートリノ$\nu_e$、ミューニュートリノ$\nu_\mu$、タウニュートリノ$\nu_\tau$といい、これらのニュートリノをフレーバー固有状態といいます。また、それぞれを質量の軽い順に並べて$\nu_1$, $\nu_2$, $\nu_3$と表記することもでき、これらのニュートリノを質量固有状態といいます（現時点ではニュートリノの質量はわかっていません）。
ただし、フレーバー固有状態と質量固有状態は１対１には対応しておらず、互いに線形結合で関係づけられます。この関係を、フレーバー固有状態を$\nu_\alpha (\alpha = e, \mu, \tau)$、質量固有状態を$\nu_j (j = 1,2,3)$として数式で表現すると
$$
\nu_\alpha = \Sigma_j U_{\alpha j} \nu_j, \quad \nu_j = \Sigma_\alpha U^\ast_{\alpha j} \nu_\alpha$$
または
$$
\begin{bmatrix}
\nu_e \newline \nu_\mu \newline \nu_\tau
\end{bmatrix} =
\begin{bmatrix}
U_{e1} & U_{e2} & U_{e3} \newline
U_{\mu 1} & U_{\mu 2} & U_{\mu 3} \newline
U_{\tau 1} & U_{\tau 2} & U_{\tau 3}
\end{bmatrix}
\begin{bmatrix}
\nu_1 \newline \nu_2 \newline \nu_3
\end{bmatrix}, \quad
\begin{bmatrix}
\nu_1 \newline \nu_2 \newline \nu_3
\end{bmatrix} =
\begin{bmatrix}
U_{e1}^* & U_{\mu 1}^* & U_{\tau 1}^* \newline
U_{e2}^* & U_{\mu 2}^* & U_{\tau 2}^* \newline
U_{e3}^* & U_{\mu 3}^* & U_{\tau 3}^*
\end{bmatrix}
\begin{bmatrix}
\nu_e \newline \nu_\mu \newline \nu_\tau
\end{bmatrix}$$
と書けます。これをニュートリノ混合といい、変換行列$U$をPMNS（Pontecorvo–Maki–Nakagawa–Sakata）行列といいます。
フレーバー固有状態－つまり質量固有状態の線形結合－にあるニュートリノは、時間とともに結合の割合が変化し、異なるフレーバー固有状態に遷移します。さらに時間がたつとニュートリノは始めのフレーバー固有状態に戻ります。
このようにニュートリノが異なるフレーバー固有状態間で遷移をくりかえす現象をニュートリノ振動といいます。ニュートリノ振動のため、ニュートリノを観測したときに得られるフレーバー固有状態の観測確率も時間とともに変化します。ニュートリノ振動は実験的に確認されており、2015年のノーベル物理学賞にも選ばれました。
今回は量子コンピュータを使ったシミュレーションが主な目的ということもあり、理論面はシンプルなほうが良かったので、以下では$\nu_e$, $\nu_\mu$ および $\nu_1$, $\nu_2$ の部分だけに限った（＝２フレーバーの）ニュートリノ振動について考えていきます。
2. 理論計算
数値計算の前に、フレーバー固有状態の観測確率の理論値を確認しておきます。２フレーバーの場合、PMNS行列$U$は2x2行列となり、その独立な成分は１つになります。それを質量固有状態の混合角$\theta$とすると、ニュートリノ混合は
$$
\begin{bmatrix}
\nu_e \newline \nu_\mu
\end{bmatrix} =
\begin{bmatrix}
\cos\theta & -\sin\theta \newline
\sin\theta & \cos\theta
\end{bmatrix}
\begin{bmatrix}
\nu_1 \newline \nu_2
\end{bmatrix}
$$
と書けます。$\nu_\alpha$として発生したニュートリノが空間を伝播した後、時刻$t$に$\nu_\beta$として観測される確率$P(\nu_\alpha \rightarrow \nu_\beta)$は
$$
\begin{align}
P(\nu_\alpha \rightarrow \nu_\beta) = |\Sigma_j U_{\beta j} Prop(\nu_j) U_{j\alpha}|^2, \quad
Prop(\nu_j) = \exp \left[ -i \frac{m_j^2 c^4}{2E\hbar} t \right]
\tag{1}
\end{align}
$$
となります。ここで、$m_j$は$\nu_j$の質量、$c$は光速、$\hbar$はDirac定数、$E$は伝播していくニュートリノのエネルギーです。例えば、観測確率$P(\nu_e \rightarrow \nu_e)$は
$$
\begin{align}
P(\nu_e \rightarrow \nu_e)
&= \left| \exp \left[ -i \frac{m_1^2c^4}{2E\hbar} t \right] \cos^2\theta + \exp \left[ -i \frac{m_2^2c^4}{2E\hbar} t \right] \sin^2\theta \right|^2 \newline
&= 1 - \sin^2(2 \theta) \sin^2 \left(\frac{\Delta m_{12}^2c^4t}{4E\hbar} \right) \tag{2}
\newline
\Delta m_{12}^2 &= |m_1^2 - m_2^2|
\end{align}
$$
と計算できます。この式から、$P(\nu_e \rightarrow \nu_e)$は周波数$\omega = \Delta m_{12}^2c^4 / 2E\hbar$で振動することがわかります。具体的な数値を代入してみると、$P(\nu_e \rightarrow \nu_e)$は下図のように振動します。
ここで、縦軸は観測確率$P(\nu_e \rightarrow \nu_e)$、横軸はニュートリノのエネルギー$E[{\rm MeV}]$です。各設定は
$$
\begin{align}
t &= \frac{L[{\rm km}]}{c[{\rm km/sec}]} = \frac{180}{3\times10^5}[\sec], \newline
\hbar &= \frac{4.135667696\times10^{-15}}{2 \pi} [{\rm eV}\cdot \sec], \newline
\Delta m_{12}^2c^4 &= 7.42\times 10^{-5}[{\rm eV}^2], \newline
\theta &= 33.44[°]
\end{align}
$$
としました。$L$はニュートリノが伝播する距離です。Dirac定数$\hbar$はCODATAを、$\Delta m_{12}^2c^4$と$\theta$はNuFITを参考にしました。文献ではKamLAND（2003年）のデータを使用していますが、今回はNuFIT（2020年）のデータを使用しました。
- CODATA
https://physics.nist.gov/cgi-bin/cuu/Results?search_for=planck+constant - NuFIT
http://www.nu-fit.org/?q=node/228
(2)式によると、２種類のニュートリノに質量差が無い（$\Delta m_{12}^2 = 0$）場合、$P(\nu_e \rightarrow \nu_e) = 1$となってニュートリノ振動は起こりません。逆に言えば、ニュートリノ振動によりニュートリノが質量をもつことがわかります。
3. シミュレータ
続いて、量子回路計算のフレームワークであるQiskitを利用して、２フレーバーのニュートリノ振動を計算します。実際の量子コンピュータで計算する前に、古典コンピュータ上で量子回路をシミュレーションしてみます。
はじめに、$\nu_e$と$\nu_\mu$を量子ビットの$|0\rangle$と$|1\rangle$に対応させます。
$$
|\nu_e\rangle
= |0\rangle
= \begin{bmatrix} 1 \newline 0 \end{bmatrix}, \quad
|\nu_\mu\rangle
= |1\rangle
= \begin{bmatrix} 0 \newline 1 \end{bmatrix}
$$
次に、(1)式を量子回路で計算するために、PMNS行列とニュートリノの伝播Propを量子ゲートとして用意します。どちらもQiskitのUゲートで表現できます（文献ではU3ゲートとU1ゲート）。
$$
PMNS: U(2\theta,0,0) =
\begin{bmatrix} \cos\theta & -\sin\theta \newline
\sin\theta & \cos\theta
\end{bmatrix}, \quad
Prop: U(0,0,\phi) =
\begin{bmatrix}
1 & 0 \newline
0 & e^{i\phi}
\end{bmatrix}$$
例えば、観測確率$P(\nu_e \rightarrow \nu_e)$は
$$
\begin{align}
P(\nu_e \rightarrow \nu_e)
&= |\langle 0 |U(2\theta,0,0)U(0,0,\phi)U^\dagger (2\theta,0,0) |0\rangle |^2 \newline
&= |\cos^2\theta + e^{i\phi} \sin^2\theta|^2 \newline
&= 1 -\sin^2(2\theta)\sin^2 \left( \frac{\phi}{2} \right)
\end{align}$$
と計算できて、$\phi = \Delta m_{12}^2c^4t / 2E\hbar$とすると(2)式に一致します。量子回路は下図のようになります。
この量子回路を古典コンピュータ上でシミュレーションした結果が下図の緑の点です。
シミュレーションにはQiskit Aerのqasm_simulatorを使用しました。
各点は、1試行あたり量子測定を4096回実行して$|0\rangle$を測定した割合（つまり$\nu_e$の観測確率）を算出し、10試行分の平均と誤差を評価したものです。理論計算と概ね一致することが確認できます。
4. 量子コンピュータ
それでは本題である量子コンピュータを利用した２フレーバーのニュートリノ振動を計算します。今回は、IBMの量子コンピュータ（IBM Quantum Experience）を利用しました。前節で作成したQiskitのコードがそのまま利用できます。
- IBM Quantum Experience (IBMQ)
https://quantum-computing.ibm.com/
IBMQの使い方は、例えば以下が参考になります。
- Quantum Native Dojo
https://dojo.qulacs.org/ja/latest/notebooks/3.2_Qiskit_IBMQ.html
IBMQでシミュレーションした結果が下図の橙の点です。
各点は、1試行あたり量子測定を4096回実行して$|0\rangle$を測定した割合（つまり$\nu_e$の観測確率）を算出し、10試行分の平均と誤差を評価したものです。 理論計算やシミュレータと同様に、$\nu_e$の観測確率が振動する様子が確認できます。理論計算からのズレが大きいのは、量子コンピュータのノイズの影響と思われます。
5. まとめ
ゲート型の量子コンピュータで２フレーバーのニュートリノ振動をシミュレーションしました。