１次元量子力学における調和振動子
任意の初期運動量分布＋任意の初期中心座標に対する時間発展

文責：遠藤理平（2011年11月23日）カテゴリ：仮想物理実験室(325)｜計算物理学(165)

ここまで、任意の初期空間分布に対する時間発展と任意の初期運動量分布に対する時間発展に対する表式の導出、時間発展の数値計算、並びにアニメーション化を行いました。シュレディンガー方程式は１階の微分方程式なので積分定数は１つで、初期状態から決めることができます。つまり、初期状態としての空間分布あるいは運動量分布を与えることで、積分定数を決めることができるわけです。逆に言えば、空間分布と運動量分布は同時に与えることはできません。量子力学に限らず電磁気学などで用いる波動方程式の場合、空間分布を与えると運動量分布も自動的に決まるため、初期状態としての空間分布あるいは運動量分布のどちらかを与えればよいことになります。しかしながら、単に初期空間分布を与えただけでは、初期運動量分布の中心 $\,p_0$ は 0 となってしまいます。同様に、単に初期運動量分布を与えただけでは空間分布の中心 $\,x_0$ も $\,0$ になります。

そこで本節では、任意の初期運動量分布かつ任意の初期中心座標に対する表式を導出し、時間発展の数値計算を行い、これまでと同様アニメーション化を行います。数値計算の具体例として、初期運動量分布が $\,p=p_0$ のガウス分布かつ、初期中心座標を $\,x=x_0$ とします。

任意の初期運動量分布かつ任意の初期中心座標に対する表式の導出

状態ベクトルの展開の展開係数

状態ベクトルを座標演算子と個数演算子の固有ベクトル $|x\langle$ 、 $|n\langle$ を用いてｍこれまでと同様

(1)

と展開するところから始めます。これまでに得られている $\langle x|n angle$ 、 $\langle n|\psi(t) angle$ を代入すると

(2)

となります。ただし、展開係数 $\psi_n(0)$ は $\,n$ の任意の関数で、後に示すとおり、状態ベクトル $|\psi(t) angle$ の初期状態から決めることができます。一旦、 $\psi_n(0)$ が決まれば、任意の時間における状態は式(5)で計算することができます。

任意の初期運動量分布に対する展開係数を得るため、式(2)の両辺に左から $\langle p|$ を掛け、 $\,t=0$ とします。

(3)

平面波と規格化エルミート多項式の直交性を考慮すると、 $\psi_n(0)$ は、

(4)

となります。任意の初期運動量分布 $\langle p | \psi(0) angle$ に対して、展開係数を得る手続きは前節の「任意の初期運動量分布に対する時間発展」と同様です。前節で示したとおり、 $\langle p | \psi(0) angle$ を $\,p=p_0$ を中心としたガウス分布として与えた場合、空間分布の中心座標は必ず $\,x=0$ です。 $\langle p | \psi(0) angle$ に複素数成分を与えることによって、空間分布の中心座標を移動させることを考えます。

平行移動演算子

$\langle p | \psi(0) angle$ に複素数成分をどのように与えることで、空間分布の中心座標を移動させることができるかわかりません。本節では、平行移動演算子を用いることで、空間分布の中心座標を意図通りに移動させることを考えます。平行移動演算子は、

(5)

で定義されます。 $\hat{U}(a)$ の座標表記における具体形は式(5)の右辺をテーラー展開することで得られます。

(6)

平行移動演算子 $\hat{U}(a)$ は

(7)

の関係式を満たすので、ユニタリー演算子です。この平行移動演算子を用いて関数を下図のように正の方向に平行移動することを考えます。

つまり、平行移動させたい関数に $\hat{U}(-a)=\hat{U}^\dagger(a)$ を作用させればよいことになります。

$\langle p | \psi(0) angle$ の与え方

先述の通り、 $\langle p | \psi(0) angle$ は $\,p$ に関して任意なので、

(8)

の変換を考えます。 $\hat{U}^\dagger(x_0)$ は座標表示における状態を正の方向に $\,x_0$ 平行移動させる意味を持ちます。式(8)の最後の変形は $\hat{p}|p angle = p|p angle$ を用いています。つまり、 $\langle p | \psi(0) angle$ の複素数成分の与え方は、平面波 $xp(-ipx_0/\hbar)$ であることがわかります。

式(8)を式(4)を代入すると、展開係数 $\psi_n(0)$ は

(9)

となります。つまり、任意の $\langle p | \psi(0) angle$ に対して、式(9)を計算することで展開係数を得ることができ、さらに、式(2)の両辺に左から $\langle x|$ を掛けることで得られる状態ベクトルの座標表記

(10)

から、時間発展を計算することができます。

ガウス型運動量分布の時間発展

運動量表記における初期状態が次のようなガウス分布の場合を考えます。

(1)

また本節では、平行移動量を $x_0=1.0 imes10^{-9}[ m m]$ 、ガウス分布の幅を決めるパラメータ $\sigma = 10^{-24}$ に対して、 $\,p_0 = 0$ 、 $p_0 = 1.0 imes 10^{-24}$ 、 $p_0 = 2.0 imes 10^{-24}$ の３つの場合について、時間発展を計算します。

積分定数 $\psi_n(0)$ の計算結果

下の図は、式(8)を用いて各運動量分布に対する $\psi_n(0)$ を計算した結果です。横軸は n、縦軸は $\psi_n(0)$ の実部と虚部をそれぞれプロットしています。

$\,p_0 = 0$

$p_0 = 1.0 imes 10^{-24}$

$p_0 = 2.0 imes 10^{-24}$

積分定数 $\psi(x,t)quiv\langle x |\psi (t) angle$ の計算結果

下図は式(9)を計算した結果です。図の横軸が x、縦軸が $\psi(x,t)$ です。赤線、緑線と青線は、それぞれ $\psi(x,t)$ の実部、虚部と絶対値を表します。アニメーションの時間間隔は、 $\Delta t = 1.0 imes 10^{-16}[ m s]$ です。

$\,p_0 = 0$

$p_0 = 1.0 imes 10^{-24}$

$p_0 = 2.0 imes 10^{-24}$

計算結果の考察

計算結果は意図通り、 $\,t=0$ の初期空間分布の中心が $\,x=x_0$ となっていることが確認できます。つまり、式(9)の展開係数の与え方が正しいことが確認できました。初期状態の与え方のポイントを整理します。

①初期空間分布の中心座標 $\,x_0$ を平面波として挿入する。
② $\langle p|\psi(0) angle$ に初期運動量分布とする。

空間分布の中心座標の最大値も、前節同様

$x_{\max} = rac{p_0}{2m} = 1.11 imes10^{-9} [ m m]$

となり、計算結果と一致しています。また、空間分布が初期状態に戻る周期は、「任意の初期状態に対する時間発展」で示したとおり、絶対値では T[s]、実部では 2T[s]であることも確認できます。ただし、 $T = 2\pi/\omega = 62.8 imes 10^{-16}[ m s]$ です。

これで本稿の目的である、１次元調和振動子における任意の初期運動量分布かつ空間分布の中心座標に対する時間発展を計算することができました。

C言語プログラムソース

複素数を含んだ関数の多重積分については、「シンプソン法を用いた複素数を含んだ２重積分(C++)」をご覧ください。

/*１次元量子力学の調和振動子における_任意の初期運動量分布＋任意の初期中心座標に対する時間発展
（2011.11.23公開）*/
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <fstream>
#include <sstream>
#include <iostream>
#include <string>
#include <cstdio>
#include <iomanip>
#include <stdio.h>
#include <complex>
#if defined(_OPENMP)
  #include <omp.h>
#endif
#if defined(_MSC_VER)
  #include <direct.h>   // Windowsフォルダ作成用
#elif defined(__GNUC__)
  #include <sys/stat.h> //  UNIX系ディレクトリ作成用
#endif
using namespace std;

const double PI = acos(-1.0);
const double e = 2.7182818284590452354;
const complex<double> I = complex<double>(0.0,1.0);
const double c = 2.99792458E+8;
const double lambda0 = 500.0E-9;
const double mu0 = 4.0*PI*1.0E-7;
const double epsilon0 = 1.0/(4.0*PI*c*c)*1.0E+7;
const double h  = 6.6260896 * 1.0E-34;
double hbar = h/(2.0*PI);
const double me = 9.10938215 * 1.0E-31;
const double eV = 1.60217733 * 1.0E-19;
const int kizami = 3000;

double omega = 1.0 * 1.0E+15;
double A = sqrt(me*omega/hbar);
double dx = 1.0 * 1.0E-11;
double Dx = 1.0 * 1.0E-9;
double dp = 1.0 * 1.0E-26;
double Dp = 1.0 * 1.0E-24;
double p0 = 0.0 * 1.0E-24;
double sigma = 1.0 * 1.0E-24;
const int n_max = 300;
complex<double> Psi_n[n_max+1];

double x0 = 1.0 * 1.0E-9;

const int ts = 0, te = 0;                   //<------------ステップ数を指定
double dt  = 1.0 * 1.0E-16;                   //<------------時間間隔を指定
int Lx = 1000;
int Lp = 1000;
string folder = "harmonic5_0.0";//作成するフォルダ名

//関数プロトタイプ
double NormalizedHermite( int n, double x );//規格化エルミート多項式
complex<double> Simpson(double a, double b, int n, double p);
complex<double> Simpson2(complex<double> **f, const int m, const int n, double h1, double h2);


//被積分関数
complex<double> Integrand(double x, double p, int n){ 
  return  sqrt(A)*NormalizedHermite(n, A*x)* exp(I*p*(x-x0)/hbar) * 1.0/sqrt(sqrt(PI)*sigma) * exp(-1.0/2.0*pow((p-p0)/sigma,2)) ;
}

char str[200];
int main(){
  #if defined(_MSC_VER)
    _mkdir(folder.c_str());   // Windowsフォルダ作成
  #elif defined(__GNUC__)
    mkdir(folder.c_str(), S_IRWXU | S_IRWXG | S_IRWXO); // UNIX系のディレクトリ作成
  #endif
  #if defined(_OPENMP)
    omp_set_num_threads(8);
    cout << "OpenMPを利用します（最大スレッド数："<< omp_get_max_threads() << "）" <<  endl;
  #endif


  string f1 = folder + "/Psi_n.data";
  ofstream ofs( f1.c_str() );
#pragma omp parallel
  {
    #pragma omp for //schedule(dynamic, 100)
    for(int n=0; n<=n_max; n++){
      double xs = -Lx/2*dx , xe = Lx/2*dx;   //xの積分範囲
      double ps = -Lp/2*dp , pe = Lp/2*dp ;  //pの積分範囲
      complex<double> **f= new complex<double>*[kizami+1];  // complex<double>型 m+1 個分の領域を動的確保
      for (int i= 0; i<=kizami; ++i) f[i] = new complex<double>[kizami+1];  // complex<double>型 n+1 個分の領域を動的確保
      double h1 = (xe-xs) / double(kizami);
      double h2 = (pe-ps) / double(kizami);
      for(int i = 0; i <= kizami; i++){
        for (int j = 0; j <= kizami; j++)  {
          double x = xs + (xe-xs)/double(kizami) * double(i);
          double p = ps + (pe-ps)/double(kizami) * double(j);
          f[i][j] = Integrand(x,p,n)/sqrt(2.0*PI*hbar);
        }
      }
      Psi_n[n] = Simpson2(f,kizami,kizami,h1,h2);
      //cout << n << " " << Psi_n[n].real() << " " << Psi_n[n].imag() << " "  << abs(Psi_n[n])<< endl;
      for (int i= 0; i<=kizami; ++i) delete[] f[i]; // 動的に確保した領域をそれぞれ解放
      delete[] f;
    }
  }
  for(int n=0; n<=n_max; n++){
    ofs << n << " " << Psi_n[n].real() << " " << Psi_n[n].imag() << " "  << abs(Psi_n[n])<< endl;
  }  
  #pragma omp parallel
  {
    #pragma omp for //schedule(dynamic, 100)
    for(int tn=ts; tn<=te; tn++){
      double t = tn * dt;
      char str[200];
      string str1;
      ofstream fout_s;
      sprintf(str, "/%d.data", tn); str1 = folder + str;
      fout_s.open(str1.c_str());
      for(int i=-Lx/2; i<=Lx/2; i++){
        double x = double(i) * dx; 
        complex<double> Psi = complex<double>(0.0,0.0);
        for(int n=0; n<=n_max; n++){
          Psi += sqrt(A)*NormalizedHermite(n, A*x) * Psi_n[n] * exp(-I*omega*(n+1.0/2.0)*t);
        }
        Psi = Psi/sqrt(A);
        fout_s << x/Dx << " " << Psi.real() << " " << Psi.imag() << " " << abs(Psi) << endl;
      }
    }
  }
}
double NormalizedHermite( int n, double x )
{
  double y0, y1, y2;
  y0 = sqrt(1.0/sqrt(PI)) * exp(-pow(x,2)/2);
  if( n==0 ) return y0;
  y1 = sqrt(2.0/sqrt(PI)) * exp(-pow(x,2)/2) * x;
  if( n==1 ) return y1;
  for(int m = 2; m<=n; m++ ){
    y2 = sqrt(2.0/double(m)) * x * y1 - sqrt(double(m-1)/double(m)) * y0;
    y0 = y1;
    y1 = y2;
  }
  return y2;
}
complex<double> Simpson(double a, double b, int n, double p)
{
  complex<double> ss1 =0.0;
  for(int i=1; i<=kizami/2-1; i++){
    double x = a + (b-a)/double(kizami) * double(2*i);
    ss1 += Integrand(x, p, n);
  }
  complex<double> ss2 =0.0;
  for(int i=1; i<=kizami/2; i++){
    double x = a + (b-a)/double(kizami) * double(2*i-1);
    ss2 += Integrand(x, p, n);
  }
  return (b-a)/(3.0*double(kizami)) * ( Integrand(a, p, n) + Integrand(b, p, n) + 2.0 * ss1 + 4.0 * ss2 );
}
complex<double> Simpson2(complex<double> **f, const int m, const int n, double h1, double h2)
{
  int i, j;
  complex<double> v;
  complex<double> *temp;
  temp = new complex<double>[m+1];  //double型 m+1 個分の領域を動的確保

  for(i = 0; i <= m; i++){
    v = - f[i][0] + f[i][n];
    for(j = 0 ; j < n - 1; j += 2)
      v += (2.0 * f[i][j] + 4.0 * f[i][j + 1]);
    temp[i] = v;
  }
  v = - temp[0] + temp[m];
  for(i = 0; i < m - 1; i += 2)  v += (2.0 * temp[i] + 4.0 * temp[i + 1]);
  delete [] temp;
  return v * h1 * h2 / 9.0;
}

ファイルたち

	【受講生募集】『natural science 科学・技術講座』新講座のご案内 2024.01.09 【大草芳江｜TOPICS】
	サイエンス・デイオブザイヤー2023贈賞式（文部科学大臣賞表彰等）を開催しました 2023.11.07 【大草芳江｜TOPICS】
	■「natural science 科学・技術講座」夏休み短期講座のご案内 2023.07.14 【遠藤理平｜TOPICS】
	「国際イノベーションコンテスト2022世界大会」でnatural science『科学・技術講座』チームが世界3等入賞 2023.07.12 【大草芳江｜TOPICS｜パブリシティ】
	サイエンス・デイ限定『学都「仙台・宮城」サイエンスマップ～光編～』第7版プレゼント申込開始（先着500名）！ 2023.07.08 【遠藤理平｜TOPICS】
	学都「仙台・宮城」サイエンス・デイ2023　7月16日（日）開催！ 2023.07.08 【大草芳江｜TOPICS】
	人材募集のご案内 2022.08.19 【大草芳江｜TOPICS】
	学都「仙台・宮城」サイエンス・デイ2022　7月17日（日）開催！ 2022.07.01 【大草芳江｜TOPICS】
	「第13回国際イノベーションコンテスト2022」国内予選大会　natural science チームが第３位入賞、４年連続通算８回目の世界大会出場へ 2022.05.23 【大草芳江｜TOPICS】
	「キャンパスベンチャーグランプリ2021」東北大会で最優秀賞、全国大会で日刊工業新聞社賞を受賞 2022.03.09 【大草芳江｜TOPICS】
	【プレスリリース】「第12回国際イノベーションコンテスト2021」世界大会　natural science チームが世界３等入賞 2022.01.17 【大草芳江｜TOPICS】
	【第８話】有限の高さの障壁へ照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.29 【遠藤理平｜仮想物理実験室】
	【第７話】無限に高い障壁に向けた電子パルスの照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.24 【遠藤理平｜仮想物理実験室】
	【第６話】無限に高い障壁へ照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.23 【遠藤理平｜仮想物理実験室】
	【第５話】電子パルスの運動アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.21 【遠藤理平｜仮想物理実験室】
	【第４話】電子パルスの作り方【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.20 【遠藤理平｜仮想物理実験室】
	【第３話】自由粒子の運動アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.19 【遠藤理平｜仮想物理実験室】
	【第２話】自由粒子のスナップショット【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.18 【遠藤理平｜仮想物理実験室】
	【第１話】プログラムの動作確認【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.17 【遠藤理平｜仮想物理実験室】
	「natural science 科学・技術講座」講師３人で設立した株式会社weCAN。アイディアの商品化第一弾「美姿勢メガネ」を朝日新聞さんに取材いただきました。 2021.08.23 【遠藤理平｜パブリシティ】

１次元量子力学における調和振動子
任意の初期運動量分布＋任意の初期中心座標に対する時間発展