量子コンピュータでニュートリノ振動をシミュレーションする

## 0. これは何の記事？
- ゲート型の量子コンピュータを使ったシミュレーションをやってみた
- シミュレーションの対象はニュートリノ振動と呼ばれる物理現象
- 下記の文献を参考に２フレーバーの場合をシミュレーションしてみた  
  Neutrino oscillations in a quantum processor  
  https://journals.aps.org/prresearch/pdf/10.1103/PhysRevResearch.1.033176

  参考文献の内容をかいつまんで説明しながら、私がやった結果を掲載します。
  なお今回作成したコードは[こちらのGitHub](https://github.com/p-circle-kuh/qc_neutrino_oscillation)にあります。

## 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)$は下図のように振動します。

![Two_flavor_result_graph_fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1219280/a84c8d42-33b2-96a1-1ef8-6ecc59bfbece.png)

ここで、縦軸は観測確率$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)式に一致します。量子回路は下図のようになります。

![circuit.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1219280/14b64697-80c2-f43c-e09b-cfb852de6cca.png)

この量子回路を古典コンピュータ上でシミュレーションした結果が下図の緑の点です。

![Two_flavor_result_graph_fig3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1219280/669b0c31-bd6b-dcf1-7519-9486f808ff13.png)
 
シミュレーションにはQiskit Aerのqasm_simulatorを使用しました。

- qasm_simulator  
  https://qiskit.org/documentation/stubs/qiskit.providers.aer.QasmSimulator.html

各点は、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でシミュレーションした結果が下図の橙の点です。

![Two_flavor_result_graph_fig4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1219280/9403ca2c-9479-7522-4257-94377fa0b268.png)

各点は、1試行あたり量子測定を4096回実行して$|0\rangle$を測定した割合（つまり$\nu_e$の観測確率）を算出し、10試行分の平均と誤差を評価したものです。 理論計算やシミュレータと同様に、$\nu_e$の観測確率が振動する様子が確認できます。理論計算からのズレが大きいのは、量子コンピュータのノイズの影響と思われます。

## 5. まとめ
ゲート型の量子コンピュータで２フレーバーのニュートリノ振動をシミュレーションしました。