HOME > natural science Laboratory > コンピュータ・シミュレーション講座 > 計算物理学

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

文責:遠藤 理平 (2011年11月23日) カテゴリ:仮想物理実験室(268)計算物理学(127)

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

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


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

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

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

(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) 
anglep=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) 
anglepに関して任意なので、

(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 = 0p_0 = 1.0	imes 10^{-24}p_0 = 2.0	imes 10^{-24} の3つの場合について、時間発展を計算します。

積分定数 \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] です。

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


C言語プログラムソース

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

/*1次元量子力学の調和振動子における_任意の初期運動量分布+任意の初期中心座標に対する時間発展
(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;
}

ファイルたち

C言語プログラム

1次元調和振動子における任意の初期運動量分布に対する時間発展

gnuplotテンプレート

Harmonic1_1D.plt
Harmonic1_1D-1.plt

Illustrator

1次元調和振動子における任意の初期運動量分布に対する時間発展

Origin

1次元調和振動子のn励起状態の運動量表示.OPJ



▲このページのトップNPO法人 natural science トップ

関連記事

仮想物理実験室







計算物理学







▲このページのトップNPO法人 natural science トップ