高校で習う二重スリット実験をシミュレーションしてみよう

## はじめに
今回は、QMsolveというPythonライブラリを使って、高校の物理で習う二重スリット実験をシミュレーションする方法を紹介します。QMsolveは、一粒子や二粒子のシュレーディンガー方程式を解くことができるライブラリで、視覚的なグラフやアニメーションを生成することができます。QMsolveの詳細については、[公式サイト](https://github.com/quantum-visualizations/qmsolve)や[PyPI](https://pypi.org/project/qmsolve/)を参照してください。
## QMsolveのインストールと使い方

QMsolveはpipで簡単にインストールできます。コマンドプロンプトで以下のコマンドを実行してください。

```
pip install qmsolve
```

3Dプロットをする場合は、Mayaviというライブラリも必要です。Mayaviもpipでインストールできますので、以下のコマンドを実行してください。

```
pip install mayavi
```

QMsolveを使うには、まずHamiltonianというクラスを使って、シミュレーションしたい系のハミルトニアンを定義します。ハミルトニアンとは、系のエネルギーを表す物理量です。QMsolveでは、ハミルトニアンの定義に必要なパラメータとして、以下のものを指定します。

- particles: シミュレーションしたい粒子の種類や数。SingleParticleやTwoParticlesなどが用意されています。
- potential: シミュレーションしたいポテンシャル（力場）の関数。粒子の位置や運動量などを引数に取り、ポテンシャルエネルギーを返す関数を定義します。
- spatial_ndim: 空間次元数。1, 2, 3のいずれかを指定します。
- N: 空間グリッドの数。空間を等間隔に分割した点の数です。
- extent: 空間グリッドの範囲。空間グリッドがカバーする空間の大きさです。

例えば、調和振動子というポテンシャルにおける一粒子系のハミルトニアンは、以下のように定義できます。

```python
import numpy as np
from qmsolve import Hamiltonian, SingleParticle, init_visualization, Å, eV

#ポテンシャルの定義
def harmonic_oscillator(particle):
    k = 100. * eV / Å**2
    return 0.5 * k * particle.x**2

#ハミルトニアンの定義
H = Hamiltonian(particles = SingleParticle(), 
                potential = harmonic_oscillator, 
                spatial_ndim = 1, N = 512, extent = 20*Å)
```

次に、Hamiltonian.solveというメソッドを使って、ハミルトニアンを対角化し、系のエネルギー準位や波動関数（固有状態）を求めます。このメソッドには、求めたい固有状態の数をmax_statesという引数で指定します。例えば、以下のようにすると、最低エネルギーから30個の固有状態を求めることができます。

```python
eigenstates = H.solve(max_states = 30)
```

このメソッドは、Eigenstatesというクラスのオブジェクトを返します。このオブジェクトには、energiesという属性があります。energiesは、求めた固有状態のエネルギーを配列で保持しています。数を返します。例えば、以下のようにすると、求めた固有状態のエネルギーを表示できます。

```python
print(eigenstates.energies) # The printed energies are in eV.
```

最後に、visualizationというクラスを使って、求めた固有状態をプロットします。visualization. slider_plotというメソッドは、スライドバーで指定したインデックスの固有状態をプロットします。

```python
visualization = init_visualization(eigenstates)
visualization.slider_plot() #interactive slider
```
![slider_plot.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3468376/b1cc07f5-a9fd-5ca8-a05a-39d7c8a53c13.gif)

以上がQMsolveの基本的な使い方です。次に、高校物理で習う基本的な量子現象を例題として解いてみましょう。

## 例題: 二重スリット実験のシミュレーション

二重スリット実験では、光や電子などの粒子が二つの細い穴から通過した後に干渉縞を作る現象が観測されます。この現象は、粒子が波動性を持つことを示しており、量子力学の基本的な実験として知られています。
前述のように、QMsolveでは、Hamiltonianクラスを使って任意のポテンシャル関数を定義することができます。二重スリット問題では、ポテンシャル関数は以下のようになります。


```python
import numpy as np
from qmsolve import Hamiltonian, SingleParticle, TimeSimulation, init_visualization, femtoseconds, m_e, Å

#二重スリットのポテンシャルを定義
def double_slit(particle):
    b = 2.0* Å # slits separation
    a = 0.5* Å # slits width
    d = 0.5* Å # slits depth

    return np.where( ((particle.x < - b/2 - a) | (particle.x > b/2 + a) | ((particle.x > -b/2)  
                     & (particle.x < b/2))) & ((particle.y < d/2) & (particle.y > -d/2) ),  1e5,  0)

#ハミルトニアンの設定
H = Hamiltonian(particles = SingleParticle(m = m_e), 
                potential = double_slit, 
                spatial_ndim = 2, N = 256, extent = 30 * Å)
```
次に、波動関数の初期値を設定します。

```python
def initial_wavefunction(particle):
    #This wavefunction correspond to a gaussian wavepacket with a mean Y momentum equal to p_y0
    σ = 1.0 * Å
    v0 = 80 * Å / femtoseconds
    p_y0 = m_e * v0
    return np.exp( -1/(4* σ**2) * ((particle.x-0)**2+(particle.y+8* Å)**2)) / np.sqrt(2*np.pi* σ**2)  *np.exp(p_y0*particle.y*1j)
```

シミュレーションの条件を設定し、計算を実行します。
```python
total_time = 0.2 * femtoseconds
sim = TimeSimulation(hamiltonian = H, method = "split-step")
store_steps = 800)
sim.run(initial_wavefunction, total_time = total_time, dt = total_time/2000., store_steps = 200)
```

最後に、シミュレーション結果のアニメーションを表示します。

```python
visualization = init_visualization(sim)

visualization.animate(xlim=[-15* Å,15* Å], ylim=[-15* Å,15* Å], potential_saturation = 0.5, wavefunction_saturation = 0.2, animation_duration = 2, fps = 30)
```

実行すると、以下のようなアニメーションが得られます。
![double_slit.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3468376/b17dd882-fbd0-1468-0c5d-b154b67d4bb2.gif)

比較のために一重スリットの実験のシミュレーションを試みます。二重スリットのポテンシャルを提議する関数のb（二つのスリットの間隔）を0に設定してシミュレーションを再度実行します。
```python
def double_slit(particle):
    b = 0.0* Å # slits separation
    a = 0.5* Å # slits width
    d = 0.5* Å # slits depth
    return np.where( ((particle.x < - b/2 - a) | (particle.x > b/2 + a) | ((particle.x > -b/2)  
                     & (particle.x < b/2))) & ((particle.y < d/2) & (particle.y > -d/2) ),  1e5,  0)
```
実行すると、以下のようなアニメーションが得られます。
![single_slit.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3468376/6ce8c999-322d-81a3-2051-5783b8fd740a.gif)


## まとめ

この記事では、QMsolveというPythonライブラリを使って、ヤングの二重スリット実験をシミュレーションしてみました。QMsolveは、一粒子や二粒子のシュレーディンガー方程式を解くことができるライブラリで、視覚的なグラフやアニメーションを生成することができます。QMsolveを使うと、基本的な量子現象をシミュレーションすることができます。例えば、調和振動子や水素原子のエネルギー準位や波動関数、量子トンネル効果などです。

QMsolveは、量子力学に興味を持つ方にとって、理解の助けになるだけでなく、楽しみにもなると思います。ぜひ、自分でパラメータを変えたり、新しいポテンシャルを定義したりして、シミュレーションの結果を観察してみてください。QMsolveで量子力学の世界を楽しく探検してください。