MiyanTarumi’s blog

MiyanTarumi’s blog

感想とか色々。数物系は書くか微妙。。。

変分モンテカルロ(VMC)をやってみる~ヘリウム原子~

Thijssenの計算物理学の12.2.2の部分の練習問題をやってみる。

www.maruzen-publishing.co.jp

ヘリウム原子

モデル

最後にヘリウム原子についても同じくVMCをやってみよう。
この系のHamiltonianはBorn-Oppenheimer近似の下で\textbf{r}_1を電子1の座標、\textbf{r}_2を電子2の座標とすると、

H=\displaystyle-\frac{1}{2}\nabla_1^2-\frac{1}{2}\nabla_2^2-\frac{2}{r_1}-\frac{2}{r_2}+\frac{1}{|\textbf{r}_1-\textbf{r}_2|}

と書ける。これまでの時と同様に試行関数を与えたいのだが、この場合は良い試行関数を作ること自体が問題になってしまうので、ここでは天下り的であるが試行関数として次のようなものを使用することにしよう。

\displaystyle \psi(\textbf{r}_1,\textbf{r}_2)=e^{-2r_1}e^{-2r_2}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}

ここでr_{12}=|\textbf{r}_1-\textbf{r}_2|とし、これまでと同様に\alphaは変分パラメータである。
次にこの試行関数を用いた時に局所エネルギーE_Lがどのようになるのかを計算してみる*1。 微分演算子の入っていない部分はそのまま係数としてかかるだけなので微分演算子の部分だけを計算すればいいだろう。そこでまず\nabla_1^2について計算することにしよう。合成関数のラプラシアンに対する作用が、

\nabla^2\left(fg\right)=\nabla\cdot\nabla\left(fg\right))=\nabla\cdot\left(\left(\nabla f\right)g+f\nabla g\right)\\ =\left(\nabla\cdot\nabla f\right)g+\nabla f\cdot\nabla g+\nabla f\cdot\nabla g+f\nabla\cdot\nabla g\\ =\nabla^2 f+2\nabla f\cdot\nabla g+f\nabla^2 g

であることに注意すると、\textbf{r}_1に関係する部分の計算は、

\nabla_1^2\left( e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right) =\left(\nabla_1^2 e^{-2r_{1}}\right)e^{\frac{r_{12}}{2(1+\alpha r_{12})}}+\nabla_1 e^{-2r_1}\cdot\nabla_1e^{\frac{r_{12}}{2(1+\alpha r_{12})}}+e^{-2r_1}\nabla_1^2e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\

を計算すればよいことになる。まず初項について計算しよう。この項の計算は簡単で1体の時のラプラシアンの計算と全く同じであるから、

\nabla_1^2 e^{-2r_{1}} =\displaystyle\left[\frac{\partial^2}{\partial r_1^2}+\frac{2}{r_1}\frac{\partial}{\partial r_1}\right]e^{-2r_1}\\ =\displaystyle\left[4-\frac{4}{r_1}\right]e^{-2r_1}

より、

\left(\nabla_1^2 e^{-2r_{1}}\right)e^{\frac{r_{12}}{2(1+\alpha r_{12})}} =\displaystyle\left[4-\frac{4}{r_1}\right]e^{-2r_1}e^{-2r_2}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}

となる。
続いて第二項について計算する。ここでこれ以降の計算の表記を簡単にするためにr_{12}が入っている部分の指数の肩を

U_{12}=\displaystyle\frac{r_{12}}{2(1+\alpha r_{12})}

と書くことにし、その微分を

DU_{12}=\displaystyle\frac{\partial}{\partial r_{12}}U_{12}=-\frac{1}{2}\frac{\alpha r_{12}}{(1+\alpha r_{12})^2}+\frac{1}{2}\frac{1}{1+\alpha r_{12}}

と書くことにする。すると第二項の計算は、

\nabla_1 e^{-2r_1}\cdot\nabla_1e^{\frac{r_{12}}{2(1+\alpha r_{12})}} =\displaystyle\sum_{i}\left(\frac{\partial}{\partial x_{1i}}e^{-2r_1}\right)\left(\frac{\partial }{\partial x_{1i}}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right)\\ =\displaystyle\sum_{i}\left(\frac{\partial r_{1}}{\partial x_{1i}}\frac{\partial}{\partial r_1}e^{-2r_1}\right)\left(\frac{\partial r_{12}}{\partial x_{1i}}\frac{\partial }{\partial r_{12}}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right)\\ =\displaystyle\sum_{i}\left[(-2)\cdot DU_{12}\right]\frac{\partial r_1}{\partial x_{1i}}\frac{\partial r_{12}}{\partial x_{1i}}e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\ =\displaystyle \sum_{i}\left[(-2)\cdot DU_{12}\right]\frac{x_{1i}}{r_1}\frac{x_{1i}-x_{2i}}{r_{12}}e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\ =\displaystyle\frac{-2DU_{12}}{r_{12}}\frac{r_1^2-\textbf{r}_1\cdot\textbf{r}_2}{r_1}e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}

となる。ここでx_{1i}は電子1の座標の成分を、x_{2i}は電子2の座標の成分をそれぞれ表すとした。
最後に第三項の計算を行おう。この項ではr_{12}の入った項の微分を二回行う必要があるのでU_{12}の二回微分を表すものとして

DDU_{12}=\displaystyle\frac{\partial^2}{\partial r_{12}^2}U_{12}=\frac{\partial}{\partial r_{12}}\left(-\frac{1}{2}\frac{\alpha r_{12}}{(1+\alpha r_{12})^2}+\frac{1}{2}\frac{1}{1+\alpha r_{12}}\right)\\ =\displaystyle\frac{\alpha^2 r_{12}}{(1+\alpha r_{12})^3}-\frac{\alpha }{(1+\alpha r_{12})^2}

とDDU_{12}を定義しよう。すると第三項は、

\nabla_1^2e^{\frac{r_{12}}{2(1+\alpha r_{12})}} =\displaystyle \sum_{i}\frac{\partial^2}{\partial x_{1i}^2}e^{\frac{r_{12}}{2(1+\alpha r_{12})}} =\displaystyle\sum_{i}\frac{\partial }{\partial x_{1i}}\left(DU_{12}\frac{x_{1i}-x_{2i}}{r_{12}}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right)\\ =\displaystyle\sum_{i}\left[\frac{DU_{12}}{r_{12}}+\left(\frac{DDU_{12}}{r_{12}}-\frac{DU_{12}}{r_{12}^2}+\frac{(DU_{12})^2}{r_{12}}\right)(x_{1i}-x_{2i})\frac{\partial r_{12}}{\partial x_{1i}}\right]e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\ =\displaystyle\sum_{i}\left[\frac{DU_{12}}{r_{12}}+\left(\frac{DDU_{12}}{r_{12}}-\frac{DU_{12}}{r_{12}^2}+\frac{(DU_{12})^2}{r_{12}}\right)\frac{(x_{1i}-x_{2i})^2}{r_{12}}\right]e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\ =\displaystyle\left[\frac{3DU_{12}}{r_{12}}+\left(\frac{DDU_{12}}{r_{12}}-\frac{DU_{12}}{r_{12}^2}+\frac{(DU_{12})^2}{r_{12}}\right)r_{12}\right]e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\ =\displaystyle\left[\frac{2DU_{12}}{r_{12}}+DDU_{12}+(DU_{12})^2\right]e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\

となる。以上から

\nabla_1^2\left( e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right)\\ =\displaystyle\left[4-\frac{4}{r_1}-\frac{4DU_{12}}{r_{12}}\frac{r_1^2-\textbf{r}_1\cdot\textbf{r}_2}{r_1}+\frac{2DU_{12}}{r_{12}}+DDU_{12}+(DU_{12})^2\right]e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}

となることがわかった。同様にして\nabla_2^2について計算すると、

\nabla_2^2\left( e^{-2r_2}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right)\\ =\displaystyle\left[4-\frac{4}{r_2}+\frac{4DU_{12}}{r_{12}}\frac{\textbf{r}_1\cdot\textbf{r}_2-r_2^2}{r_2}+\frac{2DU_{12}}{r_{12}}+DDU_{12}+(DU_{12})^2\right]e^{-2r_2}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}

を得ることができる*2。
以上から

H\psi(\textbf{r}_1,\textbf{r}_2)\\ =\displaystyle\left[-4-\frac{1}{2}\left(-\frac{4DU_{12}}{r_{12}}\frac{r_1^2-\textbf{r}_1\cdot\textbf{r}_2}{r_1}+\frac{4DU_{12}}{r_{12}}\frac{\textbf{r}_1\cdot\textbf{r}_2-r_2^2}{r_2}\right.\\ \left.+\frac{4DU_{12}}{r_{12}}+2DDU_{12}+2(DU_{12})^2\right)+\frac{1}{r_{12}}\right]\psi(\textbf{r}_1,\textbf{r}_2)

と計算でき、よって局所エネルギーE_Lは、

E_L=\displaystyle-4-\frac{1}{2}\left(-\frac{4DU_{12}}{r_{12}}\frac{r_1^2-\textbf{r}_1\cdot\textbf{r}_2}{r_1}+\frac{4DU_{12}}{r_{12}}\frac{\textbf{r}_1\cdot\textbf{r}_2-r_2^2}{r_2}\right.\\ \left.\displaystyle+\frac{4DU_{12}}{r_{12}}+2DDU_{12}+2(DU_{12})^2\right)+\displaystyle\frac{1}{r_{12}}

と求めることができた。

VMCの方法とソースコード

よってプログラミングで行う操作は次の通りである。

調和振動子のVMC

**REPAT**
始めの変分パラメータ\alphaを決める。 ランダムな場所にN個のwalkerを初期配置として配置する。
**REPEAT**
**FOR 全walker**

1. k番目のサンプリング点をX^{(k)}とする。walkerの移動後の座標としてk+1番目の候補をX^{k}を中心とする正規分布に従ってX'=X^{(k)}+\Delta Xとして提案する。
2. メトロポリステスト:
[0,1]の間の一様乱数rを生成し、r< |\psi(X')/\psi(X)|^2であれば提案を受理してX^{(k+1)}=X'とする。そうでなければ候補を棄却してX^{(k+1)}=X^{(k)}とする。
3. サンプリング点での局所エネルギーE_Lを計算する。

**END FOR**
得られた局所エネルギーの平均を取ってエネルギーの期待値を求める。
Newton法に従って変分パラメータ\alphaを更新する。
**END**

ここでXは座標であり、今の場合X=(x_1,y_1,z_1, x_2, y_2, z_2)である。
このプログラムを書いてみよう。

▶C言語のプログラム（クリックすると展開されます）

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include "MT.h"

int shokihaichi(double data1[], double data2[], double data3[]);/*プロトタイプ宣言*/
int metropolis_walk(double data1[], double data2[], double data3[], double data4[], double data5[], double data6[]);/*プロトタイプ宣言*/
int shokika(double data1[], double data2[], double data3[]);/*プロトタイプ宣言*/

double sigma=0.15;

int main(void)
{
    int sample=50000;
    int i, j;
    int acc[500];

    double alpha_0=0.05;
    double alpha_ini;
    double alpha_fin;
    double x1[500], y1[500], z1[500];
    double x2[500], y2[500], z2[500];
    double x1_backup[500], y1_backup[500], z1_backup[500];
    double x2_backup[500], y2_backup[500], z2_backup[500];
    double r1_ini[500], r2_ini[500], r12_ini[500];
    double r1_fin[500], r2_fin[500], r12_fin[500];
    double r1[500], r2[500], r12[500];
    double psi_ini[500], psi_fin[500];
    double Denom[500];/*エネルギー計算に出てくる分母の関数*/
    double U12[500], DU12[500], DDU12[500];
    double R1DR12[500], R2DR12[500];
    double Sum_E_w[500]={0};
    double E_w[500]={0};
    double E=0;
    double Sum_Edlnpsi_w[500]={0};
    double Edlnpsi_w[500]={0};
    double Edlnpsi=0;
    double Sum_dlnpsi_w[500]={0};
    double dlnpsi_w[500]={0};
    double dlnpsi=0;
    double dE;
    double h=1e-3;/*数値微分の幅*/

    init_genrand((unsigned)time(NULL));

    alpha_fin=alpha_0;

    do
    {
        alpha_ini=alpha_fin;/*前のループでの結果を始まりのalphaに代入*/

        /************************************************/
        /*dEを求める.*************************************/
        /*alpha_iniでの<E>, <Edlnpsi>, <dlnpsi>を計算する.*/
        /************************************************/

        shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*Sumの配列を初期化*/

        shokihaichi(x1, y1, z1);/*ランダムにwalkrを配置*/
        shokihaichi(x2, y2, z2);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x1, y1, z1, x1_backup, y1_backup, z1_backup);/*提案分布に従ってwalkerを変化させる*/
            metropolis_walk(x2, y2, z2, x2_backup, y2_backup, z2_backup);

            /*Metropolis Testを行う*/

            for ( j = 0; j < 500; j++)
            {
                r1_ini[j]=sqrt(x1_backup[j]*x1_backup[j]+y1_backup[j]*y1_backup[j]+z1_backup[j]*z1_backup[j]);
                r2_ini[j]=sqrt(x2_backup[j]*x2_backup[j]+y2_backup[j]*y2_backup[j]+z2_backup[j]*z2_backup[j]);
                r12_ini[j]=sqrt((x1_backup[j]-x2_backup[j])*(x1_backup[j]-x2_backup[j])+(y1_backup[j]-y2_backup[j])*(y1_backup[j]-y2_backup[j])+(z1_backup[j]-z2_backup[j])*(z1_backup[j]-z2_backup[j]));

                psi_ini[j]=exp(-4*(r1_ini[j]+r2_ini[j])+r12_ini[j]/(1+alpha_ini*r12_ini[j]));/*始めの確率分布*/

                r1_fin[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);
                r2_fin[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                r12_fin[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));

                psi_fin[j]=exp(-4*(r1_fin[j]+r2_fin[j])+r12_fin[j]/(1+alpha_ini*r12_fin[j]));/*候補の確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x1[j]=x1_backup[j];
                    y1[j]=y1_backup[j];
                    z1[j]=z1_backup[j];
                    x2[j]=x2_backup[j];
                    y2[j]=y2_backup[j];
                    z2[j]=z2_backup[j];
                }

                /*Metropolis Testの終了*/
                /*walkerについての<E>, <Edlnpsi>, <dlnpsi>を計算する.*/

                if (i>=4000)
                {
                    r1[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);/*更新し終えた位置でパラメタを計算*/
                    r2[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                    r12[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));
                    R1DR12[j]=x1[j]*(x1[j]-x2[j])+y1[j]*(y1[j]-y2[j])+z1[j]*(z1[j]-z2[j]);
                    R2DR12[j]=x2[j]*(x1[j]-x2[j])+y2[j]*(y1[j]-y2[j])+z2[j]*(z1[j]-z2[j]);
                    Denom[j]=1/(1+alpha_ini*r12[j]);
                    U12[j]=0.5*r12[j]*Denom[j];
                    DU12[j]=-0.5*alpha_ini*r12[j]*Denom[j]*Denom[j]+0.5*Denom[j];
                    DDU12[j]=alpha_ini*alpha_ini*r12[j]*Denom[j]*Denom[j]*Denom[j]-alpha_ini*Denom[j]*Denom[j];

                    Sum_E_w[j]=Sum_E_w[j]-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j];
                    E_w[j]=Sum_E_w[j]/(i-3999);

                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]))*(-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j]);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);

                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);

                }
                
            }
            
        }
        
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 500; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }

        E=E/500;
        Edlnpsi=Edlnpsi/500;
        dlnpsi=dlnpsi/500;/*dE/dalphaの計算をするための部品を計算した.*/

        dE=2*(Edlnpsi-E*dlnpsi);/*dE/dalphaの計算ができた.*/

        /*****************************/
        /*ddE/ddalphaの数値微分を行う.*/
        /****************************/

        double alpha_plus=alpha_ini+h;/*微分のためにalphaを変化させる*/
        double alpha_minus=alpha_ini-h;

        /***********************************/
        /***********************************/
        /*alpha_plusでのdE/dalphaを計算する.*/
        /***********************************/
        /***********************************/

        shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*Sumの配列を初期化*/

        shokihaichi(x1, y1, z1);/*ランダムにwalkrを配置*/
        shokihaichi(x2, y2, z2);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x1, y1, z1, x1_backup, y1_backup, z1_backup);/*提案分布に従ってwalkerを変化させる*/
            metropolis_walk(x2, y2, z2, x2_backup, y2_backup, z2_backup);

            /*Metropolis Testを行う*/

            for ( j = 0; j < 500; j++)
            {
                r1_ini[j]=sqrt(x1_backup[j]*x1_backup[j]+y1_backup[j]*y1_backup[j]+z1_backup[j]*z1_backup[j]);
                r2_ini[j]=sqrt(x2_backup[j]*x2_backup[j]+y2_backup[j]*y2_backup[j]+z2_backup[j]*z2_backup[j]);
                r12_ini[j]=sqrt((x1_backup[j]-x2_backup[j])*(x1_backup[j]-x2_backup[j])+(y1_backup[j]-y2_backup[j])*(y1_backup[j]-y2_backup[j])+(z1_backup[j]-z2_backup[j])*(z1_backup[j]-z2_backup[j]));

                psi_ini[j]=exp(-4*(r1_ini[j]+r2_ini[j])+r12_ini[j]/(1+alpha_plus*r12_ini[j]));/*始めの確率分布*/

                r1_fin[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);
                r2_fin[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                r12_fin[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));

                psi_fin[j]=exp(-4*(r1_fin[j]+r2_fin[j])+r12_fin[j]/(1+alpha_plus*r12_fin[j]));/*候補の確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x1[j]=x1_backup[j];
                    y1[j]=y1_backup[j];
                    z1[j]=z1_backup[j];
                    x2[j]=x2_backup[j];
                    y2[j]=y2_backup[j];
                    z2[j]=z2_backup[j];
                }

                /*Metropolis Testの終了*/
                /*walkerについての<E>, <Edlnpsi>, <dlnpsi>を計算する.*/

                if (i>=4000)
                {
                    r1[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);/*更新し終えた位置でパラメタを計算*/
                    r2[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                    r12[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));
                    R1DR12[j]=x1[j]*(x1[j]-x2[j])+y1[j]*(y1[j]-y2[j])+z1[j]*(z1[j]-z2[j]);
                    R2DR12[j]=x2[j]*(x1[j]-x2[j])+y2[j]*(y1[j]-y2[j])+z2[j]*(z1[j]-z2[j]);
                    Denom[j]=1/(1+alpha_plus*r12[j]);
                    U12[j]=0.5*r12[j]*Denom[j];
                    DU12[j]=-0.5*alpha_plus*r12[j]*Denom[j]*Denom[j]+0.5*Denom[j];
                    DDU12[j]=alpha_plus*alpha_plus*r12[j]*Denom[j]*Denom[j]*Denom[j]-alpha_plus*Denom[j]*Denom[j];

                    Sum_E_w[j]=Sum_E_w[j]-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j];
                    E_w[j]=Sum_E_w[j]/(i-3999);

                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]))*(-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j]);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);

                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);

                }
                
            }
            
        }
        
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 500; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }

        E=E/500;
        Edlnpsi=Edlnpsi/500;
        dlnpsi=dlnpsi/500;/*dE/dalpha_plusの計算をするための部品を計算した.*/

        double dE_plus=2*(Edlnpsi-E*dlnpsi);/*dE/dalpha_plusの計算ができた.*/

        /***********************************/
        /***********************************/
        /*alpha_minusでのdE/dalphaを計算する.*/
        /***********************************/
        /***********************************/

        shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*Sumの配列を初期化*/

        shokihaichi(x1, y1, z1);/*ランダムにwalkrを配置*/
        shokihaichi(x2, y2, z2);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x1, y1, z1, x1_backup, y1_backup, z1_backup);/*提案分布に従ってwalkerを変化させる*/
            metropolis_walk(x2, y2, z2, x2_backup, y2_backup, z2_backup);

            /*Metropolis Testを行う*/

            for ( j = 0; j < 500; j++)
            {
                r1_ini[j]=sqrt(x1_backup[j]*x1_backup[j]+y1_backup[j]*y1_backup[j]+z1_backup[j]*z1_backup[j]);
                r2_ini[j]=sqrt(x2_backup[j]*x2_backup[j]+y2_backup[j]*y2_backup[j]+z2_backup[j]*z2_backup[j]);
                r12_ini[j]=sqrt((x1_backup[j]-x2_backup[j])*(x1_backup[j]-x2_backup[j])+(y1_backup[j]-y2_backup[j])*(y1_backup[j]-y2_backup[j])+(z1_backup[j]-z2_backup[j])*(z1_backup[j]-z2_backup[j]));

                psi_ini[j]=exp(-4*(r1_ini[j]+r2_ini[j])+r12_ini[j]/(1+alpha_minus*r12_ini[j]));/*始めの確率分布*/

                r1_fin[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);
                r2_fin[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                r12_fin[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));

                psi_fin[j]=exp(-4*(r1_fin[j]+r2_fin[j])+r12_fin[j]/(1+alpha_minus*r12_fin[j]));/*候補の確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x1[j]=x1_backup[j];
                    y1[j]=y1_backup[j];
                    z1[j]=z1_backup[j];
                    x2[j]=x2_backup[j];
                    y2[j]=y2_backup[j];
                    z2[j]=z2_backup[j];
                }

                /*Metropolis Testの終了*/
                /*walkerについての<E>, <Edlnpsi>, <dlnpsi>を計算する.*/

                if (i>=4000)
                {
                    r1[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);/*更新し終えた位置でパラメタを計算*/
                    r2[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                    r12[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));
                    R1DR12[j]=x1[j]*(x1[j]-x2[j])+y1[j]*(y1[j]-y2[j])+z1[j]*(z1[j]-z2[j]);
                    R2DR12[j]=x2[j]*(x1[j]-x2[j])+y2[j]*(y1[j]-y2[j])+z2[j]*(z1[j]-z2[j]);
                    Denom[j]=1/(1+alpha_minus*r12[j]);
                    U12[j]=0.5*r12[j]*Denom[j];
                    DU12[j]=-0.5*alpha_minus*r12[j]*Denom[j]*Denom[j]+0.5*Denom[j];
                    DDU12[j]=alpha_minus*alpha_minus*r12[j]*Denom[j]*Denom[j]*Denom[j]-alpha_minus*Denom[j]*Denom[j];

                    Sum_E_w[j]=Sum_E_w[j]-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j];
                    E_w[j]=Sum_E_w[j]/(i-3999);

                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]))*(-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j]);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);

                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);

                }
                
            }
            
        }
        
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 500; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }

        E=E/500;
        Edlnpsi=Edlnpsi/500;
        dlnpsi=dlnpsi/500;/*dE/dalpha_minusの計算をするための部品を計算した.*/

        double dE_minus=2*(Edlnpsi-E*dlnpsi);/*dE/dalpha_minusの計算ができた.*/

        /***************************************************/
        /***************************************************/
        /*alphaをNewton法で求めるための部品をすべて計算し終えた*/
        /***************************************************/
        /***************************************************/

        /*alphaの更新を行う*/

        double ddE=(dE_plus-dE_minus)/(2*h);/*dE/dalphaの数値微分ddEが計算できた.*/
        
        alpha_fin=alpha_ini-dE/ddE;/*alphaを更新する.*/
        printf("alpha_fin=%.10f, alpha_ini=%.10f dE=%.10f\n", alpha_fin, alpha_ini, dE);
    } while (fabs(alpha_fin-alpha_ini)>1e-4);

    double alpha_val=alpha_fin;/*最適化した変分パラメータ*/
    printf("変分パラメータはalpha=%.10f", alpha_val);

    /*****************************/
    /*****************************/
    /*alphaでのエネルギーを計算する*/
    /*****************************/
    /*****************************/

    shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*Sumの配列を初期化*/

    shokihaichi(x1, y1, z1);/*ランダムにwalkrを配置*/
    shokihaichi(x2, y2, z2);

    for ( i = 0; i < sample; i++)
    {
        metropolis_walk(x1, y1, z1, x1_backup, y1_backup, z1_backup);/*提案分布に従ってwalkerを変化させる*/
        metropolis_walk(x2, y2, z2, x2_backup, y2_backup, z2_backup);

        /*Metropolis Testを行う*/

        for ( j = 0; j < 500; j++)
        {
            r1_ini[j]=sqrt(x1_backup[j]*x1_backup[j]+y1_backup[j]*y1_backup[j]+z1_backup[j]*z1_backup[j]);
            r2_ini[j]=sqrt(x2_backup[j]*x2_backup[j]+y2_backup[j]*y2_backup[j]+z2_backup[j]*z2_backup[j]);
            r12_ini[j]=sqrt((x1_backup[j]-x2_backup[j])*(x1_backup[j]-x2_backup[j])+(y1_backup[j]-y2_backup[j])*(y1_backup[j]-y2_backup[j])+(z1_backup[j]-z2_backup[j])*(z1_backup[j]-z2_backup[j]));

            psi_ini[j]=exp(-4*(r1_ini[j]+r2_ini[j])+r12_ini[j]/(1+alpha_val*r12_ini[j]));/*始めの確率分布*/

            r1_fin[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);
            r2_fin[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
            r12_fin[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));

            psi_fin[j]=exp(-4*(r1_fin[j]+r2_fin[j])+r12_fin[j]/(1+alpha_val*r12_fin[j]));/*候補の確率分布*/

            double metropolis=genrand_real1();

            if (psi_fin[j]/psi_ini[j]>metropolis)
            {
                acc[j]=acc[j]+1;
            }
            else
            {
                x1[j]=x1_backup[j];
                y1[j]=y1_backup[j];
                z1[j]=z1_backup[j];
                x2[j]=x2_backup[j];
                y2[j]=y2_backup[j];
                z2[j]=z2_backup[j];
            }

            /*Metropolis Testの終了*/
            /*walkerについての<E>, <Edlnpsi>, <dlnpsi>を計算する.*/

            if (i>=4000)
            {
                r1[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);/*更新し終えた位置でパラメタを計算*/
                r2[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                r12[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));
                R1DR12[j]=x1[j]*(x1[j]-x2[j])+y1[j]*(y1[j]-y2[j])+z1[j]*(z1[j]-z2[j]);
                R2DR12[j]=x2[j]*(x1[j]-x2[j])+y2[j]*(y1[j]-y2[j])+z2[j]*(z1[j]-z2[j]);
                Denom[j]=1/(1+alpha_val*r12[j]);
                U12[j]=0.5*r12[j]*Denom[j];
                DU12[j]=-0.5*alpha_val*r12[j]*Denom[j]*Denom[j]+0.5*Denom[j];
                DDU12[j]=alpha_val*alpha_val*r12[j]*Denom[j]*Denom[j]*Denom[j]-alpha_val*Denom[j]*Denom[j];

                Sum_E_w[j]=Sum_E_w[j]-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j];
                E_w[j]=Sum_E_w[j]/(i-3999);
            }
            
        }
        
    }
    
    E=0;/*変数の初期化*/
    /*walker全ての平均を取ることにする*/
    for ( j = 0; j < 500; j++)
    {
        E=E+E_w[j];
    }

    E=E/500;
    
    printf("基底状態のエネルギーは%.10f", E);
    return 0;
}

int shokihaichi(double data1[], double data2[], double data3[])
{
    int i;
    for ( i = 0; i < 500; i++)
    {
        data1[i]=(genrand_real1()-0.5)*2*2;
        data2[i]=(genrand_real1()-0.5)*2*2;
        data3[i]=(genrand_real1()-0.5)*2*2;
    }
    return 0;
}

int metropolis_walk(double data1[], double data2[], double data3[], double data4[], double data5[], double data6[])/*1~3がx,y,z, 4~6はそれらのバックアップ*/
{
    int i;
    double p[500];
    double q[500];
    double pi=asin(1e0)*2e0;
    for ( i = 0; i < 500; i++)
    {
        /*x*/
        data4[i]=data1[i];/*バックアップにxの値を保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data1[i]=data1[i]+sigma*sqrt(-2*log(p[i]))*cos(2*pi*q[i]);

        /*y*/
        data5[i]=data2[i];/*バックアップにyの値を保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data2[i]=data2[i]+sigma*sqrt(-2*log(p[i]))*cos(2*pi*q[i]);

        /*z*/
        data6[i]=data3[i];/*バックアップにxの値を保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data3[i]=data3[i]+sigma*sqrt(-2*log(p[i]))*cos(2*pi*q[i]);
    }
    return 0;
}

int shokika(double data1[], double data2[], double data3[])
{
    int i;
    for ( i = 0; i < 500; i++)
    {
        data1[i]=0;
        data2[i]=0;
        data3[i]=0;
    }
    return 0;
}

このプログラムを実行すると、\alpha=0.1433となり、その時の基底状態のエネルギーとして-2.8783が得られる。これはHartree-Fockの値-2.8617やDFTの値-2.83、厳密値-2.9037と比肩しうる数値になっている。

*1:Thijssenの教科書に載っている式は和訳も原著(2nd Edition)もこの式が間違っているので教科書の式でプログラムを組むと間違った結果が出てしまうので注意。付録のソースコードでは正しい式になってるのにどうして。。。

*2:要はr_{12}のx_{2i}での微分一回につき負号がでるので第二項のみ係数が変わる。

変分モンテカルロ(VMC)をやってみる~水素原子~

Thijssenの計算物理学の12.2.2の部分の練習問題をやってみる。

www.maruzen-publishing.co.jp

水素原子

モデル

調和振動子に続いて水素原子についても同じくVMCをやってみよう。水素原子の問題は1次元の問題に帰着できるが、ここでは3次元の問題として扱うことにする。
この系のHamiltonianはBorn-Oppenheimer近似の下で、

H=\displaystyle-\frac{1}{2}\nabla^2-\frac{1}{r}

と書ける。水素原子の基底状態の厳密解はe^{-r}であるから、ここでは試行関数として\psi(r)=e^{-\alpha r}を選ぶことにしよう。ここで\alphaは変分パラメータである。このように試行関数を決めた時に局所エネルギーはどうなるかというと、

H\psi(r) =\displaystyle\left[-\frac{1}{2}\nabla^2-\frac{1}{r}\right]e^{-\alpha r}\\ =\displaystyle\left[-\frac{1}{2}\left(\frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r}\right)-\frac{1}{r}\right]e^{-\alpha r}\\ =\displaystyle\left[-\frac{1}{r}-\frac{1}{2}\alpha\left(\alpha-\frac{2}{r}\right)\right]e^{-\alpha r}

より、

\displaystyle E_L=-\frac{1}{r}-\frac{1}{2}\alpha\left(\alpha-\frac{2}{r}\right)

となる。

VMCの方法とソースコード

よってプログラミングで行う操作は次の通りである。

調和振動子のVMC

**REPAT**
始めの変分パラメータ\alphaを決める。 ランダムな場所にN個のwalkerを初期配置として配置する。
**REPEAT**
**FOR 全walker**

1. k番目のサンプリング点をX^{(k)}とする。walkerの移動後の座標としてk+1番目の候補をX^{k}を中心とする正規分布に従ってX'=X^{(k)}+\Delta Xとして提案する。
2. メトロポリステスト:
[0,1]の間の一様乱数rを生成し、r< |\psi(X')/\psi(X)|^2=\exp(-2\alpha(r'-r))であれば提案を受理してX^{(k+1)}=X'とする。そうでなければ候補を棄却してX^{(k+1)}=X^{(k)}とする。
3. サンプリング点での局所エネルギーE_Lを計算する。

**END FOR**
得られた局所エネルギーの平均を取ってエネルギーの期待値を求める。
Newton法に従って変分パラメータ\alphaを更新する。
**END**

ここでXは座標であり、今の場合X=(x,y,z)である。
このプログラムを書いてみよう。

▶C言語のプログラム（クリックすると展開されます）

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include "MT.h"

int shokihaichi(double data1[], double data2[], double data3[]);/*プロトタイプ宣言*/
int metropolis_walk(double data1[], double data2[], double data3[], double data4[], double data5[], double data6[]);/*プロトタイプ宣言*/
int shokika(double data1[], double data2[], double data3[]);/*プロトタイプ宣言*/

double sigma=0.4;/*metropoliswalkでの正規乱数の分散*/

int main(void)
{
    int sample=30000;
    int i, j;
    int acc[400]={0};

    double alpha_0=0.8;/*変分を始める初期値*/
    double alpha_ini;/*Newton法の始まり値*/
    double alpha_fin;/*Newton法の更新値*/
    double x[400], y[400], z[400];/*walkerの位置配列*/
    double x_backup[400], y_backup[400], z_backup[400];
    double r_ini[400];
    double r_fin[400];
    double r[400];
    double psi_ini[400], psi_fin[400];
    double Sum_E_w[400]={0};
    double E_w[400]={0};
    double E=0;
    double Sum_Edlnpsi_w[400]={0};
    double Edlnpsi_w[400]={0};
    double Edlnpsi=0;
    double Sum_dlnpsi_w[400]={0};
    double dlnpsi_w[400]={0};
    double dlnpsi=0;
    double dE;
    double h=1e-3;/*数値微分の幅*/

    init_genrand((unsigned)time(NULL));/*メルセンヌツイスターの種として現在時刻を挿入*/

    alpha_fin=alpha_0;

    do
    {
        alpha_ini=alpha_fin;/*前のループの結果を始まりの値に代入*/

        /************************************************/
        /*alpha_iniでのエネルギーの微分dE/dalphaを計算する.*/
        /************************************************/

        shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*sumの配列を初期化する.*/

        shokihaichi(x, y, z);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, y, z, x_backup, y_backup, z_backup);/*提案分布に従ってwalkerを移動させる.*/

            /*Metropolis Testを行う*/
            
            for ( j = 0; j < 400; j++)
            {
                r_ini[j]=sqrt(x_backup[j]*x_backup[j]+y_backup[j]*y_backup[j]+z_backup[j]*z_backup[j]);
                psi_ini[j]=exp(-2*alpha_ini*r_ini[j]);/*元の確率分布*/

                r_fin[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                psi_fin[j]=exp(-2*alpha_ini*r_fin[j]);/*更新値での確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                    y[j]=y_backup[j];
                    z[j]=z_backup[j];
                }
                
                /*Metropolis Testの終了*/
                /*変数を更新し終えた*/

                /*walkerについて<E>や<Edlnpsi>,<dlnpsi>を計算する.*/
                if (i>=4000)
                {   
                    r[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                    Sum_E_w[j]=Sum_E_w[j]-1/r[j]-0.5*alpha_ini*(alpha_ini-2/r[j]);
                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-r[j])*(-1/r[j]-0.5*alpha_ini*(alpha_ini-2/r[j]));
                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-r[j];

                    E_w[j]=Sum_E_w[j]/(i-3999);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);
                }
            }
        }
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 400; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }
        
        E=E/400;
        Edlnpsi=Edlnpsi/400;
        dlnpsi=dlnpsi/400;/*dE/dalphaの計算をするための部品を計算した.*/

        dE=2*(Edlnpsi-E*dlnpsi);/*dE/dalphaの計算ができた.*/

        /*****************************/
        /*ddE/ddalphaの数値微分を行う.*/
        /****************************/

        double alpha_plus=alpha_ini+h;/*微分のためにalphaを変化させる*/
        double alpha_minus=alpha_ini-h;

        /***********************************/
        /***********************************/
        /*alpha_plusでのdE/dalphaを計算する.*/
        /***********************************/
        /***********************************/

        shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*sumの配列を初期化する.*/

        shokihaichi(x, y, z);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, y, z, x_backup, y_backup, z_backup);/*提案分布に従ってwalkerを移動させる.*/

            /*Metropolis Testを行う*/
            
            for ( j = 0; j < 400; j++)
            {
                r_ini[j]=sqrt(x_backup[j]*x_backup[j]+y_backup[j]*y_backup[j]+z_backup[j]*z_backup[j]);
                psi_ini[j]=exp(-2*alpha_plus*r_ini[j]);/*元の確率分布*/

                r_fin[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                psi_fin[j]=exp(-2*alpha_plus*r_fin[j]);/*更新値での確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                    y[j]=y_backup[j];
                    z[j]=z_backup[j];
                }
                
                /*Metropolis Testの終了*/
                /*変数を更新し終えた*/

                /*walkerについて<E>や<Edlnpsi>,<dlnpsi>を計算する.*/
                if (i>=4000)
                {   
                    r[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                    Sum_E_w[j]=Sum_E_w[j]-1/r[j]-0.5*alpha_plus*(alpha_plus-2/r[j]);
                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-r[j])*(-1/r[j]-0.5*alpha_plus*(alpha_plus-2/r[j]));
                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-r[j];

                    E_w[j]=Sum_E_w[j]/(i-3999);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);
                }
            }
        }
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 400; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }
        
        E=E/400;
        Edlnpsi=Edlnpsi/400;
        dlnpsi=dlnpsi/400;/*dE/dalpha_plusの計算をするための部品を計算した.*/

        double dE_plus=2*(Edlnpsi-E*dlnpsi);/*dE/dalpha_plusの計算ができた.*/

        /***********************************/
        /***********************************/
        /*alpha_minusでのdE/dalphaを計算する.*/
        /***********************************/
        /***********************************/

        shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*sumの配列を初期化する.*/

        shokihaichi(x, y, z);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, y, z, x_backup, y_backup, z_backup);/*提案分布に従ってwalkerを移動させる.*/

            /*Metropolis Testを行う*/
            
            for ( j = 0; j < 400; j++)
            {
                r_ini[j]=sqrt(x_backup[j]*x_backup[j]+y_backup[j]*y_backup[j]+z_backup[j]*z_backup[j]);
                psi_ini[j]=exp(-2*alpha_minus*r_ini[j]);/*元の確率分布*/

                r_fin[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                psi_fin[j]=exp(-2*alpha_minus*r_fin[j]);/*更新値での確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                    y[j]=y_backup[j];
                    z[j]=z_backup[j];
                }
                
                /*Metropolis Testの終了*/
                /*変数を更新し終えた*/

                /*walkerについて<E>や<Edlnpsi>,<dlnpsi>を計算する.*/
                if (i>=4000)
                {   
                    r[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                    Sum_E_w[j]=Sum_E_w[j]-1/r[j]-0.5*alpha_minus*(alpha_minus-2/r[j]);
                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-r[j])*(-1/r[j]-0.5*alpha_minus*(alpha_minus-2/r[j]));
                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-r[j];

                    E_w[j]=Sum_E_w[j]/(i-3999);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);
                }
            }
        }
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 400; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }
        
        E=E/400;
        Edlnpsi=Edlnpsi/400;
        dlnpsi=dlnpsi/400;/*dE/dalpha_minusの計算をするための部品を計算した.*/

        double dE_minus=2*(Edlnpsi-E*dlnpsi);/*dE/dalpha_minusの計算ができた.*/

        /***************************************************/
        /***************************************************/
        /*alphaをNewton法で求めるための部品をすべて計算し終えた*/
        /***************************************************/
        /***************************************************/

        /*alphaの更新を行う*/

        double ddE=(dE_plus-dE_minus)/(2*h);/*dE/dalphaの数値微分ddEが計算できた.*/
        
        alpha_fin=alpha_ini-dE/ddE;/*alphaを更新する.*/

        printf("alpha_fin=%.10f, alpha_ini=%.10f\n", alpha_fin, alpha_ini);

    } while (fabs(alpha_fin-alpha_ini)>1e-6);

    double alpha_va=alpha_fin;/*変分で最適化した変分パラメータ*/

    printf("エネルギーが最小になる変分パラメータは%.10f\n", alpha_va);

    /******************************************************/
    /*alpha_vaでの変分で求めた基底状態のエネルギーを計算しよう*/
    /******************************************************/
    
    shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*sumの配列を初期化する.*/

    shokihaichi(x, y, z);

    for ( i = 0; i < sample; i++)
    {
        metropolis_walk(x, y, z, x_backup, y_backup, z_backup);/*提案分布に従ってwalkerを移動させる.*/

        /*Metropolis Testを行う*/
        
        for ( j = 0; j < 400; j++)
        {
            r_ini[j]=sqrt(x_backup[j]*x_backup[j]+y_backup[j]*y_backup[j]+z_backup[j]*z_backup[j]);
            psi_ini[j]=exp(-2*alpha_va*r_ini[j]);/*元の確率分布*/

            r_fin[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
            psi_fin[j]=exp(-2*alpha_va*r_fin[j]);/*更新値での確率分布*/

            double metropolis=genrand_real1();

            if (psi_fin[j]/psi_ini[j]>metropolis)
            {
                acc[j]=acc[j]+1;
            }
            else
            {
                x[j]=x_backup[j];
                y[j]=y_backup[j];
                z[j]=z_backup[j];
            }
            
            /*Metropolis Testの終了*/
            /*変数を更新し終えた*/

            /*walkerについて<E>や<Edlnpsi>,<dlnpsi>を計算する.*/
            if (i>=4000)
            {   
                r[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                Sum_E_w[j]=Sum_E_w[j]-1/r[j]-0.5*alpha_va*(alpha_va-2/r[j]);

                E_w[j]=Sum_E_w[j]/(i-3999);
            }
        }
    }
    E=0;/*変数の初期化*/
    /*walker全ての平均を取ることにする*/
    for ( j = 0; j < 400; j++)
    {
        E=E+E_w[j];
    }
    
    E=E/400;

    printf("基底状態のエネルギーは%.10f\n", E);


    return 0;
}

int shokihaichi(double data1[], double data2[], double data3[])
{
    int i;
    for ( i = 0; i < 400; i++)
    {
        data1[i]=(genrand_real1()-0.5)*2*2;
        data2[i]=(genrand_real1()-0.5)*2*2;
        data3[i]=(genrand_real1()-0.5)*2*2;
    }
    return 0;
}

int metropolis_walk(double data1[], double data2[], double data3[], double data4[], double data5[], double data6[])/*1~3がx,y,z, 4~6はそれらのバックアップ*/
{
    int i;
    double p[400];
    double q[400];
    double pi=asin(1e0)*2e0;
    for ( i = 0; i < 400; i++)
    {
        /*x*/
        data4[i]=data1[i];/*バックアップにxの値を保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data1[i]=data1[i]+sigma*sqrt(-2*log(p[i]))*cos(2*pi*q[i]);

        /*y*/
        data5[i]=data2[i];/*バックアップにyの値を保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data2[i]=data2[i]+sigma*sqrt(-2*log(p[i]))*cos(2*pi*q[i]);

        /*z*/
        data6[i]=data3[i];/*バックアップにxの値を保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data3[i]=data3[i]+sigma*sqrt(-2*log(p[i]))*cos(2*pi*q[i]);
    }
    return 0;
}

int shokika(double data1[], double data2[], double data3[])
{
    int i;
    for ( i = 0; i < 400; i++)
    {
        data1[i]=0;
        data2[i]=0;
        data3[i]=0;
    }
    return 0;
}

このプログラムを実行すると、\alpha=1の時にエネルギーが最小になり、その時のエネルギーは\langle E\rangle=-1/2になることがわかる。この結果は厳密解と一致している。

変分モンテカルロ(VMC)をやってみる~調和振動子~

Thijssenの計算物理学の12.2.2の部分の練習問題をやってみる。

www.maruzen-publishing.co.jp

調和振動子

モデル

まずは一番簡単な例として1次元調和振動子についてVMCをやってみることにしよう。
この系のHamiltonianは、

H=\displaystyle -\frac{1}{2}\frac{d^2}{dx^2}+\frac{1}{2}x^2

である。よく知られているようにこのHamiltonianの基底状態の解析解は\exp(-x^2/2)であるから、ここでは試行関数として\psi(x)=\exp(-\alpha x^2)を用いることにしよう。ただし\alphaは変分パラメータである。
このように試行波動関数を決めたとしてHamiltonianに対する局所エネルギーE_Lを計算してみよう。これは単純に計算すればよく、

H\psi(x)=\left[\displaystyle -\frac{1}{2}\frac{d^2}{dx^2}+\frac{1}{2}x^2\right]e^{-\alpha x^2}\\ =\left[\displaystyle-\frac{1}{2}\frac{d}{dx}\left(-2\alpha x\right)+\frac{1}{2}\frac{d}{dx}\left(-2\alpha x\right)\right]e^{-\alpha x^2}\\ =\left[\displaystyle \alpha+x^2\left(\frac{1}{2}-2\alpha^2\right)\right]e^{-\alpha x^2}

であるから、局所エネルギーE_Lはその定義から

E_L=\displaystyle \alpha+x^2\left(\frac{1}{2}-2\alpha^2\right)

と求まる。

VMCの方法とソースコード

よってプログラミングで行う操作は次の通りである。

調和振動子のVMC

**REPAT**
始めの変分パラメータ\alphaを決める。 ランダムな場所にN個のwalkerを初期配置として配置する。
**REPEAT**
**FOR 全walker**

1. k番目のサンプリング点をx^{(k)}とする。walkerの移動後の座標としてk+1番目の候補をx^{k}を中心とする正規分布に従ってx'=x^{(k)}+\Delta xとして提案する。
2. メトロポリステスト:
[0,1]の間の一様乱数rを生成し、r< |\psi(x')/\psi(x)|^2=\exp(-2\alpha(x'^2-x^2))であれば提案を受理してx^{(k+1)}=x'とする。そうでなければ候補を棄却してx^{(k+1)}=x^{(k)}とする。
3. サンプリング点での局所エネルギーE_Lを計算する。

**END FOR**
得られた局所エネルギーの平均を取ってエネルギーの期待値を求める。
Newton法に従って変分パラメータ\alphaを更新する。
**END**

このプログラムを書いてみよう。

▶C言語のプログラム（クリックすると展開されます）

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include "MT.h"

int shoki(double data1[]);/*プロトタイプ宣言*/
int metropolis_walk(double data1[], double data2[], double stepsize);/*プロトタイプ宣言*/
int shokika(double data1[], double data2[], double data3[]);/*プロトタイプ宣言*/

int main(void)
{
    int i, j;
    int sample=30000;

    double alpha_ini;/*Newton法のためのalphaの元値*/
    double alpha_fin;/*alphaの更新値*/
    double alpha_plus;
    double alpha_minus;
    double alpha_0=0.4;/*alphaの初期値*/
    double stepsize=0.7;
    double x[400];
    double x_backup[400];
    double psi_ini[400];
    double psi_fin[400];
    double ene_walkersum[400]={0};
    double ene_walker[400]={0};
    double E=0;
    double EDpsi_sum[400]={0};
    double EDpsi[400]={0};
    double Edlnpsidalpha=0;
    double Dpsi_sum[400]={0};
    double Dpsi[400]={0};
    double dlnpsidalph=0;
    double acc[400]={0};
    double dEdalpha;
    double h=1e-3;/*数値微分の幅*/

    init_genrand((unsigned)time(NULL));

    alpha_fin=alpha_0;
    alpha_ini=alpha_fin+1;/*while文に引っかかるようにずらしておく*/

    while (fabs(alpha_fin-alpha_ini)>1e-6)
    {
        /*alphaを変化させながら最小になるようなalphaを探しに行く*/
        alpha_ini=alpha_fin;/*前ループの結果を代入*/

        shokika(ene_walkersum, EDpsi_sum, Dpsi_sum);/*sumを収納する変数の初期化を行う.*/
        
        /*メトロポリス法を行ってenergyの微分を求める*/
        shoki(x);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, x_backup, stepsize);/*walkさせる*/

            /*MetropolisTestを行う*/
            for ( j = 0; j < 400; j++)
            {
                psi_ini[j]=exp(-2*alpha_ini*x_backup[j]*x_backup[j]);
                psi_fin[j]=exp(-2*alpha_ini*x[j]*x[j]);

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                }
                /*MetropolisTestの終了*/

                /*各walkerでの<E>, <E dln(psi)/d(alpha)>, <dln(psi)/d(alpha)>の計算 */
                if (i>=4000)
                {
                    ene_walkersum[j]=ene_walkersum[j]+(alpha_ini+x[j]*x[j]*(0.5-2*alpha_ini*alpha_ini));
                    ene_walker[j]=ene_walkersum[j]/(i-3999);

                    EDpsi_sum[j]=EDpsi_sum[j]+(-x[j]*x[j])*(alpha_ini+x[j]*x[j]*(0.5-2*alpha_ini*alpha_ini));
                    EDpsi[j]=EDpsi_sum[j]/(i-3999);

                    Dpsi_sum[j]=Dpsi_sum[j]-x[j]*x[j];
                    Dpsi[j]=Dpsi_sum[j]/(i-3999);
                }
                 
            }
        
        }
        /*walker全ての平均を取る*/
        /*初期化を行う*/
        E=0;
        Edlnpsidalpha=0;
        dlnpsidalph=0;
        for ( j = 0; j < 400; j++)
        {
            E=E+ene_walker[j];
            Edlnpsidalpha=Edlnpsidalpha+EDpsi[j];
            dlnpsidalph=dlnpsidalph+Dpsi[j];
        }
        E=E/400;/*400個の平均だから割っておく*/
        Edlnpsidalpha=Edlnpsidalpha/400;
        dlnpsidalph=dlnpsidalph/400;
        
        dEdalpha=2*(Edlnpsidalpha-E*dlnpsidalph);/*dE/dalphaが計算できた*/


        alpha_plus=alpha_ini+h;
        alpha_minus=alpha_ini-h;

        /***********************************************************************/
        /*dE/dalphaの数値微分を求めるためにalpha+とalpha-でのdE/dalphaの値を求める.*/
        /***********************************************************************/


        /****************************************/
        /****************************************/
        /*alpha+でのdE/dalphaを計算することにする.*/
        /****************************************/
        /****************************************/



        shoki(x);
        shokika(ene_walkersum, EDpsi_sum, Dpsi_sum);/*再利用するためにsumを収納する変数の初期化を行う.*/

        double acc_plus[400]={0};
        double dEdalpha_plus;

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, x_backup, stepsize);/*walkさせる*/

            /*MetropolisTestを行う*/
            for ( j = 0; j < 400; j++)
            {
                psi_ini[j]=exp(-2*alpha_plus*x_backup[j]*x_backup[j]);
                psi_fin[j]=exp(-2*alpha_plus*x[j]*x[j]);

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc_plus[j]=acc_plus[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                }
                /*MetropolisTestの終了*/

                /*各walkerでの<E>, <E dln(psi)/d(alpha)>, <dln(psi)/d(alpha)>の計算 */
                if (i>=4000)
                {
                    ene_walkersum[j]=ene_walkersum[j]+(alpha_plus+x[j]*x[j]*(0.5-2*alpha_plus*alpha_plus));
                    ene_walker[j]=ene_walkersum[j]/(i-3999);

                    EDpsi_sum[j]=EDpsi_sum[j]+(-x[j]*x[j])*(alpha_plus+x[j]*x[j]*(0.5-2*alpha_plus*alpha_plus));
                    EDpsi[j]=EDpsi_sum[j]/(i-3999);

                    Dpsi_sum[j]=Dpsi_sum[j]-x[j]*x[j];
                    Dpsi[j]=Dpsi_sum[j]/(i-3999);
                }
                 
            }
        
        }
        /*walker全ての平均を取る*/
        /*初期化を行う*/
        E=0;
        Edlnpsidalpha=0;
        dlnpsidalph=0;
        for ( j = 0; j < 400; j++)
        {
            E=E+ene_walker[j];
            Edlnpsidalpha=Edlnpsidalpha+EDpsi[j];
            dlnpsidalph=dlnpsidalph+Dpsi[j];
        }
        E=E/400;/*400個の平均だから割っておく*/
        Edlnpsidalpha=Edlnpsidalpha/400;
        dlnpsidalph=dlnpsidalph/400;
        
        dEdalpha_plus=2*(Edlnpsidalpha-E*dlnpsidalph);/*dE/dalpha+が計算できた*/



        /****************************************/
        /****************************************/
        /*alpha-でのdE/dalphaを計算することにする.*/
        /****************************************/
        /****************************************/


        shoki(x);
        shokika(ene_walkersum, EDpsi_sum, Dpsi_sum);/*再利用するためにsumを収納する変数の初期化を行う.*/

        double acc_minus[400]={0};
        double dEdalpha_minus;

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, x_backup, stepsize);/*walkさせる*/

            /*MetropolisTestを行う*/
            for ( j = 0; j < 400; j++)
            {
                psi_ini[j]=exp(-2*alpha_minus*x_backup[j]*x_backup[j]);
                psi_fin[j]=exp(-2*alpha_minus*x[j]*x[j]);

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc_plus[j]=acc_plus[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                }
                /*MetropolisTestの終了*/

                /*各walkerでの<E>, <E dln(psi)/d(alpha)>, <dln(psi)/d(alpha)>の計算 */
                if (i>=4000)
                {
                    ene_walkersum[j]=ene_walkersum[j]+(alpha_minus+x[j]*x[j]*(0.5-2*alpha_minus*alpha_minus));
                    ene_walker[j]=ene_walkersum[j]/(i-3999);

                    EDpsi_sum[j]=EDpsi_sum[j]+(-x[j]*x[j])*(alpha_minus+x[j]*x[j]*(0.5-2*alpha_minus*alpha_minus));
                    EDpsi[j]=EDpsi_sum[j]/(i-3999);

                    Dpsi_sum[j]=Dpsi_sum[j]-x[j]*x[j];
                    Dpsi[j]=Dpsi_sum[j]/(i-3999);
                }
                 
            }
        
        }
        /*walker全ての平均を取る*/
        /*初期化を行う*/
        E=0;
        Edlnpsidalpha=0;
        dlnpsidalph=0;
        for ( j = 0; j < 400; j++)
        {
            E=E+ene_walker[j];
            Edlnpsidalpha=Edlnpsidalpha+EDpsi[j];
            dlnpsidalph=dlnpsidalph+Dpsi[j];
        }
        E=E/400;/*400個の平均だから割っておく*/
        Edlnpsidalpha=Edlnpsidalpha/400;
        dlnpsidalph=dlnpsidalph/400;
        
        dEdalpha_minus=2*(Edlnpsidalpha-E*dlnpsidalph);/*dE/dalpha-が計算できた*/

        double ddE=(dEdalpha_plus-dEdalpha_minus)/(2*h);/*dE/dalphaのalpha_iniでの数値微分*/

        alpha_fin=alpha_ini-dEdalpha/ddE;/*alphaの更新*/

        printf("%.10f %.10f %.10f\n", alpha_fin, alpha_ini, dEdalpha);

    }

    double alpha_va=alpha_fin;/*変分で最小化された変分パラメータ*/
    double E_val;/*変分で求めたエネルギー*/

    /*alpha_valでのエネルギーを計算*/

    shoki(x);
    shokika(ene_walkersum, EDpsi_sum, Dpsi_sum);/*再利用するためにsumを収納する変数の初期化を行う.*/

    double acc_va[400]={0};


    for ( i = 0; i < sample; i++)
    {
        metropolis_walk(x, x_backup, stepsize);/*walkさせる*/

        /*MetropolisTestを行う*/
        for ( j = 0; j < 400; j++)
        {
            psi_ini[j]=exp(-2*alpha_va*x_backup[j]*x_backup[j]);
            psi_fin[j]=exp(-2*alpha_va*x[j]*x[j]);

            double metropolis=genrand_real1();

            if (psi_fin[j]/psi_ini[j]>metropolis)
            {
                acc_va[j]=acc_va[j]+1;
            }
            else
            {
                x[j]=x_backup[j];
            }
            /*MetropolisTestの終了*/

            /*各walkerでの<E>, <E dln(psi)/d(alpha)>, <dln(psi)/d(alpha)>の計算 */
            if (i>=4000)
            {
                ene_walkersum[j]=ene_walkersum[j]+(alpha_va+x[j]*x[j]*(0.5-2*alpha_va*alpha_va));
                ene_walker[j]=ene_walkersum[j]/(i-3999);
            }
             
        }
    
    }
    /*walker全ての平均を取る*/
    /*初期化を行う*/
    E=0;
    for ( j = 0; j < 400; j++)
    {
        E=E+ene_walker[j];
    }
    E=E/400;/*400個の平均だから割っておく*/
    E_val=E;

    printf("変分で求めたエネルギーは%.10f\n", E_val);
    return 0;
}

int shoki(double data1[])
{
    int i;
    for ( i = 0; i < 400; i++)
    {
        data1[i]=(genrand_real1()-0.5)*2*2;/*[-2,2]の間にランダムに初期*/
    }
    return 0;
}

int metropolis_walk(double data1[], double data2[], double stepsize)/*data1がx, data2がx_backup*/
{
    int i;
    double p[400];
    double q[400];
    double pi=asin(1e0)*2e0;
    for ( i = 0; i < 400; i++)
    {
        data2[i]=data1[i];/*data1の情報をdata2に保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data1[i]=data1[i]+stepsize*sqrt(-2*log(p[i]))*cos(2*pi*q[i]);
    }
    return 0;
}

int shokika(double data1[], double data2[], double data3[])
{
    int i;
    for ( i = 0; i < 400; i++)
    {
        data1[i]=0;
        data2[i]=0;
        data3[i]=0;
    }
    return 0;
}

このプログラムを実行すると、\alpha=0.5が得られ、その時にエネルギーは、\langle E\rangle=1/2になる。これは厳密解と一致する。

変分モンテカルロ(VMC)をやってみる~準備~

Thijssenの計算物理学の12.2.2の部分の練習問題をやってみる。

www.maruzen-publishing.co.jp

変分モンテカルロ法(VMC)

量子力学の多くの問題においてSchrödinger方程式を厳密に解くことは難しく、何らかの形で近似的に解く必要がある。そこで変分による方法で基底状態について近似的に解くことを考えてみよう。

VMCのための前準備

量子力学での変分原理

HamiltonianをHとし、そのエネルギー固有値をE_nとしたときにそれに対応した固有関数を\phi_n(X)（Xは座標）と書くことにしよう*1。
ここで唐突だが、ある波動関数\Psi(X,a)（aはパラメータ）を考えてみよう。この波動関数をエネルギー固有関数でモード展開すると次のように展開できる：

\Psi(X,a)=\displaystyle\sum_{n}\Psi_n\phi_n(X)

ここで\Psi_nは\phi_n(X)で展開するときの展開係数であり、それは

\Psi_n=\int dX\, \phi_n^{*}(X)\Psi(X,a)

で与えられる。
量子力学においてエネルギーの期待値\langle E\rangleは、

\langle E\rangle=\langle \Psi|H|\Psi\rangle=\int dX\,\Psi^{*}(X,a)H\Psi(X,a)

で計算できるのであった。ここにモード展開された表式での波動関数を代入してあげると、

\displaystyle \langle E\rangle =\int dX\,\sum_{n,m}\Psi_{m}^{*}\Psi_{n}\phi_m^{*}(X)H\phi_{n}(X)\\ =\displaystyle\sum_{n,m}\Psi_{m}^{*}\Psi_{n}\int dX\,\phi_m^{*}(X)E_{n}\phi_{n}(X)\\ =\displaystyle\sum_{n,m}E_{n}\Psi_{m}^{*}\Psi_{n}\delta_{n,m}\\ =\displaystyle\sum_{n}E_n|\Psi_n|^2\\ \geq\displaystyle\sum_{n}E_{0}|\Psi_n|^2\\ =E_0

となる。ここで\Psi(X,a)が規格化されているとして規格化条件\sum_{n}|\Psi_n|^2=1を用いた。
ここから波動関数の期待値\langle E\rangleは基底状態のエネルギーE_0を下回ることはないことがわかる。なので\Psi(X,a)のパラメータaを期待値\langle E\rangleが最小になるように調整すれば波動関数は基底状態の波動関数に近づくことになり、基底状態のエネルギーの近似値を得ることができる。

量子力学における変分原理

試行波動関数\Psi(X,a)（aは変分パラメータ）について、

\delta\left(\int dX\,\Psi^*(X,a)H\Psi(X,a)\right)=0

となるように変分パラメータを決めると基底状態の良い近似になる。

さて、ここまでは試行波動関数\Psi(X,a)が規格化されているものとして話を進めてきたが、一般に試行波動関数として複雑なものを採用した場合は規格化を行うことは面倒であるし、そもそも積分を解析的に行うことが可能かどうかすら怪しいであろう。そこでこの先では試行波動関数があらかじめには規格化されていないとし、試行波動関数|\Psi(a)\rangleを(||\Psi(a)||\rangle)^{-1/2}|\Psi(a)\rangleと書くことにしよう。すると、エネルギー期待値は、

\displaystyle \langle E\rangle=\frac{\langle\Psi(a)|H|\Psi(a)\rangle}{\langle\Psi(a)|\Psi(a)\rangle}

と書き直される。よって変分原理は、

量子力学における変分原理

試行波動関数\Psi(X,a)（aは変分パラメータ）について、

\delta\left(\displaystyle\frac{\int dX\,\Psi^*(X,a)H\Psi(X,a)}{\int dX\,\Psi^*(X,a)\Psi(X,a)}\right)=0

となるように変分パラメータを決めると基底状態の良い近似になる。

のように書き直せばよい。
原理的にはこの変分原理に従ってエネルギーが最小になるように変分パラメータaを調整すればよいのだが、先に述べたように試行関数やHamiltonianが複雑になれば解析的に実行することは困難であり、積分を区分求積的に数値計算で行うにしても1次元系ならともかく3次元などになれば問題が複雑になればなるほど積分の次元が飛躍的に大きくなり次元の呪いに嵌ってしまう。例えば3次元の2粒子系（Helium原子の電子など）では積分に必要な次元は6次元になる。
そこで高次元でも数値計算を行うためにマルコフ連鎖モンテカルロ法を用いて変分における積分計算を行うことにする（だから変分"モンテカルロ"という。）。

マルコフ連鎖モンテカルロ法(MCMC)

前節での要請に従ってMCMCを導入する。

モンテカルロ法

先にも述べたように高次元の数値積分は難しい。このことを実感するために領域として長さ1のn次元超立方体を積分することを考えてみよう。仮に[0,1]の間を1000分割して区分求積的に積分を行うことを考えると、積分領域は10^{3n}個の領域に分けられることになり、考えなければならない領域の数が指数関数的に増加してしまい高次元では計算不可能になってしまうことがわかるだろう。
この計算コストの指数的増加による次元の呪いから逃れるために積分領域からランダムに値を選んできてそれらの平均をとることで積分を求めるということを考えてみる。同様に[0,1]^nの領域での積分を考えよう。ランダムに点を選んで持ってきて平均を取るというアイデアに従えば積分は、

I=\displaystyle\frac{1}{N}\sum_{i=1}^{N}f(x_i)\hspace{10pt}(x_i\mbox{はランダムに選ぶ})

と表される。この式自体はx_iをランダムに選ぶことを除けば区分求積とほぼ同じ形なのでNが大きくなれば厳密な積分に近づくことは容易にわかる。そこでこの近づき方がどのように振る舞うかを確認してみよう。そのためにはこの式の分散の振る舞いを見ればよい。Iの分散は次の式で表せる。

\sigma^2=\displaystyle\left\langle\left(\frac{1}{N}\sum_{i=1}^{N}f(x_i)\right)^2\right\rangle-\left(\left\langle\frac{1}{N}\sum_{i=1}^{N}f(x_i)\right\rangle\right)^2\\

この式において初項のi\neq jの時は後ろの項と相殺して消え、i=jのときに初項からf^2の期待値がN個でることを考えて計算すると、

\sigma^2=\displaystyle\frac{1}{N}\left(\left\langle f^2\right\rangle-\left\langle f\right\rangle^2\right)

となり、\sigma\sim1/\sqrt{N}のように次元によらない形で振る舞うことがわかる。
以上からモンテカルロ法を用いて積分計算をすれば次元の呪いを回避しながら数値計算できそうだと予想できるだろう。

マルコフ連鎖モンテカルロ法

さて、上の例ではランダムに点を選ぶときに点が選ばれる確率は一様にサンプリングされるとした。しかし、高次元の積分を行うときには一様にサンプリングすることは計算資源の無駄遣いになってしまう。というのも多くの被積分関数は領域によって積分への寄与が全く異なり、次元が上がれば積分への寄与が大きい適切なサンプリング点がサンプリングされる確率が下がってしまうからである。例えばGauss積分e^{-\sum_{i=1}^{n}x_i^2/2}を考え、サンプリングは[-10000,10000]^nから行うとする。よく知られているようにGauss関数は0に近いところでのみ大きな値を持ち原点から離れるにしたがって急速に0に近づく。そこで有意な値を持つ範囲を|x|>2であるとすると、n=1のときの有意なサンプリング領域は2\times10^{-2}\%である。しかし、n=10では有意な領域は全積分領域のうち4\times10^{-6}\%しかない。こうしたことから高次元で安易に一様サンプリングをしてしまうと大きな統計誤差を含んだ結果が出てしまう。つまり一見モンテカルロ法を用いて次元の呪いを回避できたように見えてサンプリングの問題として残ってしまっている。
このようなサンプリングの問題を踏まえると、サンプリング数を高次元でも膨大なものにしないために一様にサンプリングを行うのではなく、積分領域のうち有意な領域を重点的にサンプリングすればいいだろう。この有意な領域への重点的なサンプリングを行うために望んでいる確率分布を構成する方法の一つとしてマルコフ連鎖モンテカルロ法(MCMC)がある。

構成したい確率分布をP(\{x\})=P(x_1,\dots,x_n)としよう。MCMCではx_1, \dots, x_nの値を各点が現れる滞在時間をP(\{x\})に比例するように変化させるようにし、これによって十分長い時間で統計平均の良い近似を得ることができる。そこでこのようなMCMCの方法として今回用いるメトロポリス法によるMCMCを導入することにしよう*2。

メトロポリス法

MCMCの一つの方法としてメトロポリス法の手順を確認する。
先と同様に構成したい確率分布をP(x)とする。メトロポリス法の手続きは次のようなものである。

メトロポリス法

1. k番目のサンプリング点をx^{(k)}とする。実数\Delta xをランダムに選び、k+1番目の候補をx'=x^{(k)}+\Delta xとして提案する。ただし、MCMCの条件を満たすために\Delta x:と-\Delta xが同じ確率で現れるように候補の選び方を提案する。
2. メトロポリステスト:
[0,1]の間の一様乱数rを生成し、r< P(x')/P(x)であれば提案を受理してx^{(k+1)}=x'とする。そうでなければ候補を棄却してx^{(k+1)}=x^{(k)}とする(この条件は\mathrm{min}(1, P(x')/P(x))で提案を受理するということと同じである。)。

1において述べたように\Delta xと\Delta -\Delta xが出現する確率が同じであればよいのでその提案確率は[-c,c]の間の一様乱数でもx^{k}を中心としたガウス乱数でも何でもよい。なお、記事でVMCをする際にはガウス乱数を用いて計算することにする。

VMC

このようにして得られるメトロポリス法によるMCMCを用いてどのように変分計算を行うかを考えてみよう。
量子力学のおける変分原理によれば、試行関数を最適化するためには、

\langle E\rangle =\displaystyle\frac{\int dX\,\Psi^*(X,a)H\Psi(X,a)}{\int dX\,\Psi^*(X,a)\Psi(X,a)}

を最小化すればよいのであった。唐突だがこの式を少し変形してみよう。形を確率分布を挟んだ積分らしく変形すると、\langle E\rangleは局所エネルギーE_Lを用いて

\langle E\rangle =\displaystyle\frac{\int dX\,\Psi^*(X,a)\Psi(X,a)\displaystyle\frac{H\Psi(X,a)}{\Psi(X,a)}}{\int dX\,\Psi^*(X,a)\Psi(X,a)}\\ =\displaystyle\int dX\,\frac{|\Psi(X,a)|^2}{\int dX'\,|\Psi(X',a)|^2}\displaystyle\frac{H\Psi(X,a)}{\Psi(X,a)}\\ =\int dX\,\rho(X)E_L

と書くことができる。ここで\rho(X)とE_Lについて

\rho(X)= \displaystyle\int dX\,\frac{|\Psi(X,a)|^2}{\int dX'\,|\Psi(X',a)|^2}
E_L= \displaystyle\frac{H\Psi(X,a)}{\Psi(X,a)}

とした。
このように書き換えることでエネルギーの期待値\langle E\rangleは確率分布が\rho(X)である局所エネルギーE_Lの期待値と思うことができる。そこでMCMCを用いてサンプリング点を\rho(X)から選び出して平均を取ってやれば大数の法則からエネルギーの期待値\langle E\rangleを求めることができるであろう。これがVMCにおける多次元積分の基本的な方針である。
最後にここまでの話を踏まえてVMCでのメトロポリス法を整理しておこう。

VMC

**REPAT**
始めの変分パラメータaを決める。 ランダムな場所にN個のwalkerを初期配置として配置する。
**REPEAT**
**FOR 全walker**

1. k番目のサンプリング点をX^{(k)}とする。walkerの移動後の座標としてk+1番目の候補を提案分布に従ってX'=X^{(k)}+\Delta Xとして提案する。ただし、MCMCの条件を満たすために\Delta X:と-\Delta Xが同じ確率で現れるように候補の選び方を提案する。
2. メトロポリステスト:
[0,1]の間の一様乱数rを生成し、r< \rho(X')/\rho(X)=|\Psi(X')/\Psi(X)|^2であれば提案を受理してX^{(k+1)}=X'とする。そうでなければ候補を棄却してX^{(k+1)}=X^{(k)}とする。
3. サンプリング点での局所エネルギーE_Lを計算する。

**END FOR**
得られた局所エネルギーの平均を取ってエネルギーの期待値などを求める。
Newton法に従って変分パラメータaを更新する。
**END**

Newton法

最後に変分パラメータの調整の仕方を確認する。
エネルギーが最小になる\alphaでは\langle E\rangleのaによる微分が0になるのでdE/daのゼロ点を求めればよい。E自体を数値微分したものを使ってもよいが二階導関数を数値微分の数値微分で求めることになるので大きな誤差が出てしまう。そこで一階導関数は解析的に計算することにしよう。Eの表式に戻って微分を行うと、

\displaystyle\frac{d\langle E\rangle}{da} =\displaystyle\frac{d}{da}\frac{\langle\Psi(a)|H|\Psi(a)\rangle}{\langle\Psi(a)|\Psi(a)\rangle}\\ =\displaystyle\frac{d\langle \Psi(a)|}{da}\frac{H|\Psi(a)\rangle}{\langle \Psi(a)|\Psi(a)\rangle}+\frac{d}{da}\frac{\langle\Psi(a)|H}{\langle\Psi(a)|\Psi(a)\rangle}\frac{d|\Psi(a)\rangle}{da}\\ -\displaystyle\frac{\langle\Psi(a)|H|\Psi(a)\rangle}{\left(\langle \Psi(a)|\Psi(a)\rangle\right)^2}\left(\frac{d\langle \Psi(a)|}{da}|\Psi(a)\rangle+\langle \Psi(a)|\frac{d|\Psi(a)\rangle}{da}\right)

となる。ここで試行関数として実関数のみを用いることを仮定すると、計算結果として現れるのは実数であることに注意すると、

\displaystyle\langle \Psi(a)|H\frac{d|\Psi(a)\rangle}{da}=\left(\langle \Psi(a)|H\frac{d|\Psi(a)\rangle}{da}\right)^{*}=\frac{d\langle \Psi(a)}{da}|H\Psi(a)\rangle
\displaystyle\langle \Psi(a)|\frac{d|\Psi(a)\rangle}{da}=\left(\langle \Psi(a)|\frac{d|\Psi(a)\rangle}{da}\right)^{*}=\frac{d\langle \Psi(a)}{da}|\Psi(a)\rangle

と右側に落とした微分を左側に持ってくることができる。よってdE/daは、

\displaystyle\frac{d\langle E\rangle}{da}\\ =\displaystyle2\left(\frac{d\langle \Psi(a)|}{da}\frac{H|\Psi(a)\rangle}{\langle \Psi(a)|\Psi(a)\rangle}-\displaystyle\frac{\langle\Psi(a)|H|\Psi(a)\rangle}{\left(\langle \Psi(a)|\Psi(a)\rangle\right)^2}\frac{d\langle \Psi(a)|}{da}|\Psi(a)\rangle\right)\\ =\displaystyle 2\left(\int dX\frac{\Psi(X,a)H\Psi(X,a)}{\int dX'\Psi(X',a)^2}\frac{1}{\Psi(X,a)}\frac{d\Psi(X,a)}{da}\right.\\ \left.\displaystyle-\int dX\frac{\Psi(X,a)H\Psi(X,a)}{\int dX'\Psi(X',a)^2}\int dX\frac{\Psi(X,a)^2}{\int dX'\Psi(X',a)^2}\frac{1}{\Psi(X,a)}\frac{d\Psi(X,a)}{da}\right)\\ =\displaystyle 2\left(\int dX\frac{\Psi(X,a)H\Psi(X,a)}{\int dX'\Psi(X',a)^2}\frac{d\ln \Psi(X,a)}{da}\right.\\ \left.\displaystyle-\int dX\frac{\Psi(X,a)H\Psi(X,a)}{\int dX'\Psi(X',a)^2}\int dX\frac{\Psi(X,a)^2}{\int dX'\Psi(X',a)^2}\frac{d\ln \Psi(X,a)}{da}\right)\\ =\displaystyle 2\left(\left\langle E_L\frac{d\ln\Psi(a)}{da}\right\rangle-\left\langle E_L\right\rangle\left\langle\frac{d\ln \Psi(a)}{da}\right\rangle\right)

と計算できる。これをMCMCを用いて数値計算し、あとはNewton法の更新式

\displaystyle a_{\mathrm{new}}=a_{\mathrm{old}}-\left.\gamma\frac{d\langle E\rangle}{da}\right|_{a=a_{\mathrm{old}}}

に従って変分パラメータaを更新すればよい。ただし\gammaはdE/daの数値微分である。

さて、これで準備は終わりなので以上を踏まえて具体的な系で計算をしてみよう。

参考文献

[1]. J. M. Thijssen, 「計算物理学」,丸善出版, 2012.
[2]. 花田 政範, 松浦 壮, 「ゼロからできるMCMC マルコフ連鎖モンテカルロ法の実践的入門」, 講談社, 2020.

*1:nが大きくなる順にエネルギー固有値が並べられているとしよう。

*2:もちろんMCMCになるような一般条件があるのだが、本題ではないのでこの記事では述べない。

2022年度　大阪大学　大学院理学研究科　物理　解答

阪大の理学研究科　物理学・宇宙地球科学専攻の2022年度院試問題の解答です．間違いがあれば教えてください．
問題は，ここにあります．（9/10現在はまだ上がってませんが，そのうちアップロードされると思います．）

問題I

I

(1)

密度分布\rho(\textbf{r})はデルタ関数と次で定義される関数 f(x)

x = \left\{ \begin{array}{} 1 & 0\leq x\leq 2a\\ 0 & \mbox{otherwise} \end{array} \right.

を用いると，

\displaystyle \rho(\textbf{r})=\frac{M}{2a}f(x)\delta(y)\delta(z)

と書ける．よって棒の慣性モーメントI_Oは，回転軸をz軸とすると，

\displaystyle I_O=\int dV\,\rho(\textbf{r})(x^2+y^2)=\int_{0}^{2a}\frac{M}{2a}dx\, x^2\\ \displaystyle=\frac{4}{3}Ma^2

と求められる．

(2)

棒に働く力のモーメントが-aMg\sin\varphiなことから回転の運動方程式は，\varphiが十分小さいので\sin\varphi\sim\varphiと近似できることに注意すると，

\displaystyle I_O\frac{d^2\varphi}{dt^2}=-aMg\sin\varphi\sim -aMg\varphi\\ \displaystyle \rightarrow \frac{d^2\varphi}{dt^2}=-\frac{aMg}{I_O}\varphi

となる．よって棒の振り子運動の角振動数\omegaは，

\displaystyle \omega=\sqrt{\frac{aMg}{I_O}}

とわかる．

(3)

(2)で得た運動方程式の一般解は，

\varphi=A\sin(\omega t+\delta)

と表せる．ここでA,\deltaは積分定数である．
まずt=0で\varphi=0であることから\delta=0であることが従う．さらに角運動量と力積の関係から

\displaystyle \left.I_O\frac{d\varphi}{dt}\right|_{t=0}=aP

が成立するので，

I_O\omega A\cos(0)=aP

よりAは，

\displaystyle A=\frac{aP}{I_O\omega}=P\sqrt{\frac{a}{I_OMg}}

と求められる．以上から\varphiは，

\displaystyle \varphi=P\sqrt{\frac{a}{I_OMg}}\sin\left(\sqrt{\frac{aMg}{I_O}}t\right)

と決まる．
よって振れ幅の最大値\varphi_0は，

\displaystyle \varphi_0=P\sqrt{\frac{a}{I_OMg}}

と求められる．

II

(4)

重心の座標(x',y')=(x+a\sin\theta,-a\cos\theta)についてこの時間微分は

(\dot{x}',\dot{y}')=(\dot{x}+a\dot{\theta}\cos\theta,a\dot{\theta}\sin\theta)

であるから，重心の運動エネルギーは

\displaystyle \frac{1}{2}M((\dot{x}')^2+(\dot{y}')^2)=\frac{1}{2}M(\dot{x}^2+a^2\dot{\theta}^2+2a\dot{x}\dot{\theta}\cos\theta)

となる．よって台Aと棒からなる系の運動エネルギーはこれと棒の重心周りの回転エネルギー，そして台の運動エネルギーを足して

\displaystyle T=\frac{1}{2}I\dot{\theta}^2+\frac{1}{2}M(\dot{x}^2+a^2\dot{\theta}^2+2a\dot{x}\dot{\theta}\cos\theta)+\frac{1}{2}m\dot{x}^2\\ \displaystyle =\frac{1}{2}(M+m)\dot{x}^2+\frac{1}{2}(I+Ma^2)\dot{\theta}^2+Ma\cos\theta\dot{x}\dot{\theta}

となる．
一方で位置エネルギーについて基準面をy=0にとることと，y軸を鉛直上向きに取ることに注意すると，

U=Mgy'=-Mga\cos\theta

となる．

(5)

\thetaについて十分小さいとし，\displaystyle \cos\theta\sim1-\frac{1}{2}\theta^2と近似することにする．\theta,\dot{\theta}に関する3次以上の項及び定数項を無視することにすれば，LagrangianLは，

\displaystyle L=T-U \sim\frac{1}{2}(M+m)\dot{x}^2+\frac{1}{2}(I+Ma^2)\dot{\theta}^2+Ma\dot{x}\dot{\theta}-\frac{1}{2}Mga\theta^2

となる．よって，\mbox{（あ）}=\displaystyle\frac{1}{2}(M+m)，\mbox{（い）}=\displaystyle\frac{1}{2}(I+Ma^2)，\mbox{（う）}=\displaystyle Ma，\mbox{（え）}=\displaystyle\frac{1}{2}Mgaとわかる．

(6)

xについてのLagrange方程式を求めると

\displaystyle \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{x}}\right)=\frac{\partial L}{\partial x}\\ \rightarrow \displaystyle (M+m)\ddot{x}+Ma\ddot{\theta}=0

となる．一方で\thetaについてのLagrange方程式は，

\displaystyle \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}}\right)=\frac{\partial L}{\partial \theta}\\ \rightarrow \displaystyle (I+Ma^2)\ddot{\theta}+Ma\ddot{x}=-Mga\theta

となる．

(7)

m\to\inftyの時を考える．この時xについてのLagrange方程式において効くのはmがかかっている項のみであるから，この極限で運動方程式は

m\ddot{x}=0

となる．t=0において台がx=0にあり，そして速度が0であることを踏まえると台は動かない．
\ddot{x}=0を\thetaに関する運動方程式に代入すれば

(I+Ma^2)\ddot{\theta}=-Mga\theta

となる．よって振り子運動の角振動数\omegaは，

\displaystyle \omega=\sqrt{\frac{Mga}{I+Ma^2}}

とわかる．

(8)

m=0の時を考えよう．このときxについての運動方程式は

M\ddot{x}+Ma\ddot{\theta}=0\\ \rightarrow \ddot{x}+a\ddot{\theta}=0

となる．棒の重心に注目する．x'の二回微分について計算してみると，

\ddot{x}' =\frac{d}{dt}(\dot{x}+a\dot{\theta}\cos\theta)=\ddot{x}+a\ddot{\theta}\cos\theta-a\dot{\theta}^2\sin\theta\\ \sim \ddot{x}+a\ddot{\theta}\left(1-\frac{1}{2}\theta^2\right)-a\dot{\theta}^2\theta\\ \sim\ddot{x}+a\ddot{\theta}=0

となる．重心の初速度が0であることを考えると，台Aと棒は固定された重心のx'の周りで運動する．
最後に棒の振り子運動の角振動数\omegaを求める．xについての運動方程式から\ddot{x}=-a\ddot{\theta}であるから，これを\thetaに関する運動方程式に代入すると，

(I+Ma^2)\ddot{\theta}+Ma(-a\ddot{\theta})=-Mga\theta\\ \rightarrow I\ddot{\theta}=-Mga\theta

となる．よって角振動数\omegaは，

\displaystyle \omega=\sqrt{\frac{Mga}{I}}

である．

(9)

まず，xについてのの運動方程式から

\displaystyle \ddot{x}=-\frac{Ma}{M+m}\ddot{\theta}

となる．これを\thetaについての運動方程式に代入すると，

\displaystyle (I+Ma^2)\ddot{\theta}+Ma\left(-\frac{Ma}{M+m}\ddot{\theta}\right)=-Mga\theta\\ \displaystyle \rightarrow \frac{(M+m)I+mMa^2}{M+m}\ddot{\theta}=-Mga\theta

と変形できる．よって\thetaは角振動数

\displaystyle \omega=\sqrt{\frac{(M+m)Mga}{(M+m)I+mMa^2}}

の単振動をする．ゆえに，この方程式の一般解は積分定数A,\deltaを用いて

\theta=A\sin(\omega t+\delta)

と書ける．初期条件についてt=0において\theta=\theta_0，\dot{\theta}=0であることを踏まえれば，A=\theta_0，\delta=\frac{\pi}{2}と決まる．以上からまず\thetaのついて解が

\theta=\theta_0\cos\omega t

となることがわかる．
一方xについての運動方程式について両辺を積分すると，v_0，x_0を積分定数として

\displaystyle x=-\frac{Ma}{M+m}\theta+v_0t+x_0\\ \displaystyle =-\frac{Ma}{M+m}\theta_0\cos\omega t+v_0t+x_0

となる．t=0においてx=0，\dot{x}=0であることを考えると，v_0=0，x_0=\frac{Ma}{M+m}\theta_0であることがわかる．よってxは，

\displaystyle x=-\frac{Ma}{M+m}\theta_0\cos\omega t+\frac{Ma}{M+m}\theta_0

と求まる．

問題2

I

(1)

円筒座標系として(\rho,\theta,z)を取る．この時に電荷分布\rho(\textbf{r})は，

\rho(\textbf{r})=q\delta(\rho-a)\delta(z)

と表される．基底ベクトルを用いて\textbf{r}と\textbf{r}'を\textbf{r}=z\textbf{e}_z，\textbf{r}'=\rho'\textbf{e}_{\rho}+z'\textbf{e}_{z}と表すことにすれば，

|\textbf{r}-\textbf{r}'|=\sqrt{\rho'^2+(z-z')^2}

となる．よって静電ポテンシャル\phi(\textbf{r})は，

\displaystyle \phi(\textbf{r})=\int dV'\,\frac{1}{4\pi\epsilon_0}\frac{\rho(\textbf{r}')}{|\textbf{r}-\textbf{r}'|} =\frac{q}{4\pi\epsilon_0}\int_{0}^{\infty}d\rho'\int_{0}^{2\pi}d\theta'\int_{-\infty}^{\infty}dz'\frac{\rho'\delta(\rho'-a)\delta(z')}{\sqrt{\rho'^2+(z-z')^2}}\\ =\displaystyle \frac{q}{2\epsilon_0}\frac{a}{\sqrt{a^2+z^2}}

と求められる．

(2)

\textbf{E}(\textbf{r})=-\nabla\phi(\textbf{r})であるから，円筒座標におけるナブラ演算子が

\displaystyle \nabla=\textbf{e}_\rho\frac{\partial}{\partial \rho}+\frac{\textbf{e}_\theta}{\rho}\frac{\partial}{\partial \theta}+\textbf{e}_z\frac{\partial}{\partial z}

であることと，\phi(\textbf{r})がzにのみ依存することに注意すると，

\displaystyle \textbf{E}(\textbf{r})=-\nabla\phi(\textbf{r})=-\textbf{e}_z\frac{\partial}{\partial z}\frac{q}{2\epsilon_0}\frac{a}{\sqrt{a^2+z^2}}=\frac{q}{2\epsilon_0}\frac{az}{(a^2+z^2)^{\frac{3}{2}}}

と電場が求められる．原点において電場は明らかに0であるから電荷Q'が受ける力は\textbf{F}=\textbf{0}である．

(3)

\textbf{r}'=(a\cos\varphi,a\sin\varphi)であるからds'=ad\varphiである．

\textbf{r}-\textbf{r}'=(x-a\cos\varphi,y-a\sin\varphi,z)

なので，

\displaystyle \textbf{E}(\textbf{r})=\frac{1}{4\pi\epsilon_0}\int_{0}^{2\pi}d\varphi\,\frac{qa(\textbf{r}-\textbf{r}')}{|\textbf{r}-\textbf{r}'|^3} =\frac{qa}{4\pi\epsilon_0}\int_{0}^{2\pi}d\varphi\,\frac{1}{a^3}\left[1+\frac{3}{a}(x\cos\varphi+y\sin\varphi)\right]\left( \begin{array}{} x-a\cos\varphi\\ y-a\sin\varphi\\ z \end{array} \right)

を計算すればよい．ここで\cos\varphi，\sin\varphi，\cos\varphi\sin\varphiの積分は0になることを踏まえてnon zeroな項のみを抜き出せば，

\displaystyle \textbf{E}(\textbf{r})=\frac{q}{4\pi\epsilon_0a^2}\int_{0}^{2\pi}d\varphi\,\left( \begin{array}{} x-3x\cos^2\varphi\\ y-3y\sin^2\varphi\\ z \end{array} \right) =\displaystyle \frac{q}{4\pi\epsilon_0a^2}\left(\begin{array}{} -\pi x\\ -\pi y\\ 2\pi z \end{array} \right)\\ =\displaystyle\frac{q}{4\epsilon_0 a^2} \left(\begin{array}{} -x\\ -y\\ 2z \end{array} \right)

と求められる．念のために(2)の答えとの整合性を確認しておくと|z|\ll aの下で近似を行い最小次の項のみ残せば(3)での解と一致する．

II

(4)

同じく円筒座標を用いて考える．\textbf{r}=(x,y,z)，\textbf{r}'=(a\cos\varphi,a\sin\varphi)とすると，

\displaystyle |\textbf{r}-\textbf{r}'|^{-1}=\frac{1}{\sqrt{(x-a\cos\varphi)^2+(y-a\sin\varphi)^2+z^2}}=\frac{1}{\sqrt{r^2+a^2-2ax\cos\varphi-2ay\sin\varphi}}\\ \displaystyle =\frac{1}{a}\frac{1}{\sqrt{1-2\frac{x}{a}\cos\varphi-2\frac{y}{a}\sin\varphi+\left(\frac{r}{a}\right)^2}}\\ \displaystyle \sim\frac{1}{a}\left(1+\frac{1}{a}\left(x\cos\varphi+y\sin\varphi\right)\right)

と近似できるので，ベクトルポテンシャル\textbf{A}(\textbf{r})は，d\textbf{r}'=(-a\sin\varphi,a\cos\varphi,0)d\varphiであることから，

\displaystyle \textbf{A}(\textbf{r}) =\frac{\mu_0I}{4\pi}\int_{0}^{2\pi} \frac{ad\varphi}{|\textbf{r}-\textbf{r}'|}\left( \begin{array}{} -\sin\varphi\\ \cos\varphi\\ 0 \end{array} \right)\\ \displaystyle \sim\frac{\mu_0I}{4\pi}\int_{0}^{2\pi}ad\varphi \frac{1}{a}\left(1+\frac{1}{a}\left(x\cos\varphi+y\sin\varphi\right)\right) \left( \begin{array}{} -\sin\varphi\\ \cos\varphi\\ 0 \end{array} \right)

となる．再びここでもnon zeroな項のみ抜き出してやることにすれば，

\displaystyle \textbf{A}(\textbf{r}) =\frac{\mu_0I}{4\pi a}\int_{0}^{2\pi}d\varphi\left( \begin{array}{} -y\sin^2\varphi\\ x\cos^2\varphi\\ 0 \end{array} \right)\\ \displaystyle =\frac{\mu_0 I}{4a}\left(\begin{array}{} -y\\ x\\ 0 \end{array} \right)

とベクトルポテンシャルが計算できる．

(5)

\textbf{B}(\textbf{r})=\nabla\times\textbf{A}(\textbf{r})であるから，これを計算すればよい．磁束密度\textbf{B}のx，y成分については\textbf{A}がzに依存しないことやA_zが0であることが容易にわかるのでB_zのみ求めればよいであろう．実際に計算してみると，

\displaystyle B_z(\textbf{r})=\frac{\partial A_y}{\partial x}-\frac{\partial A_x}{\partial y} =\frac{\mu_0 I}{4a}\left(\frac{\partial x}{\partial x}-\frac{\partial-y}{\partial y}\right)\\ \displaystyle =\frac{\mu_0I}{2a}

となる．以上から磁束密度は，

\displaystyle \textbf{B}(\textbf{r})=\frac{\mu_0I}{2a}\textbf{e}_z

と計算できた．

(6)

円周上の電荷qは1秒間に\displaystyle\frac{v}{2\pi a}=\frac{a\omega}{2\pi a}=\frac{\omega}{2\pi}回，回転することになる．ここで\omegaは角速度である．また角運動量 \textbf{L}の大きさLについては，慣性モーメント\mathcal{I}_Cを用いてL=\mathcal{I}_C\omegaと書けるのであった．よって電流Iは，

\displaystyle I=q\frac{\omega}{2\pi}=q\frac{L}{\mathcal{I}_C}\frac{1}{2\pi}=\frac{qL}{2\pi\mathcal{I}_C}

と表せる．よって磁束密度についてこれを代入すれば

\displaystyle \textbf{B}(\textbf{r})=\frac{\mu_0I}{2a}\textbf{e}_z=\frac{\mu_0}{2a}\frac{qL}{2\pi\mathcal{I}_C}\textbf{e}_z=\frac{q}{4\pi a\mathcal{I}_C}\textbf{L}

と\textbf{B}と\textbf{L}の関係がわかる最後に磁気モーメント\textbf{M}と磁束密度\textbf{B}の相互作用エネルギーVについて，その定義V=-\textbf{M}\cdot\textbf{B}より，

\displaystyle V=-\textbf{M}\cdot\textbf{B}=-\frac{q}{4\pi a\mathcal{I}_C}\textbf{M}\cdot\textbf{L}

となる．

III

(7)

\textbf{B}=\nabla\times\textbf{A}であるから，磁束密度\textbf{B}は，

\displaystyle B_x=\frac{\partial A_z}{\partial y}-\frac{\partial A_y}{\partial z}=\frac{\mu_0Ia^2x}{4}\left(0-\left(-3\frac{z}{r^5}\right)\right)=\frac{3\mu_0Ia^2}{4r^5}xz
\displaystyle B_y=\frac{\partial A_x}{\partial z}-\frac{\partial A_z}{\partial x}=-\frac{\mu_0Ia^2y}{4}\left(-\left(3\frac{z}{r^5}\right)-0\right)=\frac{3\mu_0Ia^2}{4r^5}yz
\displaystyle B_z=\frac{\partial A_y}{\partial x}-\frac{\partial A_x}{\partial y}=\frac{\mu_0Ia^2}{4}\left(\frac{\partial}{\partial x}\left(\frac{x}{r^3}\right)-\frac{\partial }{\partial y}\left(\frac{-y}{r^3}\right)\right)\\ \displaystyle=\frac{\mu_0Ia^2}{4}\left(\frac{2}{r^3}-3\frac{x^2+y^2}{r^5}\right)\\ \displaystyle=\frac{\mu_0Ia^2}{4}\left(3\frac{z^2}{r^5}-\frac{1}{r^3}\right)

と求められる．さて，

\textbf{m}_1\cdot\textbf{r}=\pi Ia^2z

であるから

\displaystyle B_x=\frac{3\mu_0Ia^2}{4r^5}xz=3\frac{\mu_0}{4\pi r^5}\pi Ia^2zx=\frac{\mu_0}{4\pi r^5}3(\textbf{m}_1\cdot\textbf{r})x
\displaystyle B_y=\frac{3\mu_0Ia^2}{4r^5}yz=3\frac{\mu_0}{4\pi r^5}\pi Ia^2zy=\frac{\mu_0}{4\pi r^5}3(\textbf{m}_1\cdot\textbf{r})y
\displaystyle B_z=\frac{\mu_0Ia^2}{4}\left(3\frac{z^2}{r^5}-\frac{1}{r^3}\right)=\frac{\mu_0}{4\pi r^5}\left(3\pi Ia^2z^2-\pi Ia^2r^2\right)=\frac{\mu_0}{4\pi r^5}\left(3(\textbf{m}_1\cdot\textbf{r})z-\pi Ia^2r^2\right)

とも書ける．よって最終的に

\displaystyle \textbf{B}(\textbf{r})=B_x\textbf{e}_x+B_y\textbf{e}_y+B_z\textbf{e}_z =\frac{\mu_0}{4\pi r^5}\left[3(\textbf{m}_1\cdot\textbf{r})(x\textbf{e}_x+y\textbf{e}_y+z\textbf{e}_z)-\pi Ia^2\textbf{e}_zr^2\right]\\ \displaystyle =\frac{\mu_0}{4\pi r^5}\left[3(\textbf{m}_1\cdot\textbf{r})\textbf{r}-\textbf{m}_1r^2\right]

となる．

(8)

わからない．

問題3

I

(1)

\phi_1(0)=\phi_2(0)であるから

1+r=t
(2)

V(x)=-V_0\delta(x)をSchrödinger方程式に代入して区間[-\epsilon,\epsilon]で積分すると

\displaystyle -\frac{\hbar^2}{2m}\int_{-\epsilon}^{\epsilon}dx\,\frac{d^2\phi}{dx^2}-V_0\int_{-\epsilon}^{\epsilon}dx\,\delta(x)\phi(x)=E\int_{-\epsilon}^{\epsilon}dx\,\phi(x)\\ \displaystyle \rightarrow -\frac{\hbar^2}{2m}\left[\frac{d\phi_2(\epsilon)}{dx}-\frac{d\phi_1(-\epsilon)}{dx}\right]-V_0\phi_2(0)=E\int_{-\epsilon}^{\epsilon}dx\,\phi(x)

となる．ここでデルタ関数がかかった項の積分について\phi_1(0)でも\phi_2(0)でもどちらでも結果は変わらないが簡単な方を選んだ．ここで\epsilon\to0のlimitをとると右辺の積分は0になることから，上の式は

\displaystyle -\frac{\hbar^2}{2m}\left[\frac{d\phi_2(\epsilon)}{dx}-\frac{d\phi_1(-\epsilon)}{dx}\right]-V_0\phi_2(0)=0\\ \displaystyle \rightarrow -\frac{\hbar^2}{2m}\left[iqt-iq(1-r)\right]-V_0t=0\\ \displaystyle \rightarrow t+r=i\frac{2mV_0}{\hbar^2q}t+1

と変形される．\omega=\hbar^2q/mV_0を用いると，接続条件は

\displaystyle t+r=\frac{2i}{\omega}t+1

となる．

(3)

(1)と(2)の結果をtとrについて解くと

\displaystyle t=\frac{\omega}{\omega-i}=\frac{\omega(\omega+i)}{\omega^2+1}
\displaystyle r=\frac{i(\omega+i)}{\omega^2+1}

となる．よって，

\displaystyle |t|^2 =\frac{\omega^2(\omega^2+1)}{(\omega^2+1)^2}=\frac{\omega^2}{\omega^2+1}
\displaystyle |r|^2 =\frac{\omega^2+1}{(\omega^2+1)^2}=\frac{1}{\omega^2+1}

となる．

(4)

q=\sqrt{2mE}/\hbarであることと，\omega=\hbar^2q/mV_0なことからE\to0のとき\omega\to0であり，E\to\inftyであるとき\omega\to\inftyに対応する．
まず，E\to0のとき\omega\to0であったから

\displaystyle \lim_{\omega\to0}|t|^2=\lim_{\omega\to0}\frac{\omega^2}{\omega^2+1}=0
\displaystyle \lim_{\omega\to0}|r|^2 =\lim_{\omega\to0}\frac{1}{\omega^2+1}=1

となる．これは入射粒子のエネルギーが0であるときには全く透過しないであろうと予想されることから考えれば当然の結果である．
一方でE\to\inftyのときは\omega\to\inftyであったから，

\displaystyle \lim_{\omega\to\infty}|t|^2=\lim_{\omega\to\infty}\frac{\omega^2}{\omega^2+1}=1
\displaystyle \lim_{\omega\to\infty}|r|^2 =\lim_{\omega\to\infty}\frac{1}{\omega^2+1}=0

となる．この結果もエネルギーが\inftyであるときにはすべて壁を乗り越えうることを思えば妥当な結果になっている．
最後にE=\mathrm{const.}であるときにV_0\to-V_0とするとどうなるかを考える．この操作によっての変化は(2)の式を踏まえるとtとrの虚部の符号が反転することのみである．つまり，これは複素共役をとったことに対応するが，複素共役をとったことによって絶対値は変化しないので|r|^2と|t|^2はこの操作では変化しない．

II

(5)

x=0とx=aにおけるそれぞれのポテンシャル障壁でのrやtはこの場合でも変化しないと考えられる．
まずB_1についてはA_1が反射したときの寄与とD_1が透過したときの寄与が考えられ，それらを足し合わせたものがB_1になる．ゆえに，

B_1=rA_1+tD_1

となる．
続いてC_1については，A_1が透過したときの寄与とD_1が反射したときの寄与が考えることができる．よって，

C_1=tA_1+rD_1

が成立する．

(6)

まずD_1について

D_1=-\frac{r}{t}A_1+\frac{1}{t}B_1

とまとめる．これを先の二式目に代入すると，

\displaystyle C_1=tA_1+r\left(-\frac{r}{t}A_1+\frac{1}{t}B_1\right)=t\left(1-\left(\frac{r}{t}\right)^2\right)A_1+\frac{r}{t}B_1

となる．以上から行列Zは，

\displaystyle Z=\left(\begin{array}{} \displaystyle t\left(1-\left(\frac{r}{t}\right)^2\right) & \displaystyle \frac{r}{t}\\ \displaystyle -\frac{r}{t} & \displaystyle \frac{1}{t} \end{array}\right)

と求められる．

(7)

\psi_2(x)に対する2つの表式について等置すれば

A_2e^{iq(x-a)}+B_2e^{-iq(x-a)}=C_1e^{iqx}+D_1e^{-iqx}

となり，

C_1=A_2e^{-iqa}
D_1=B_2e^{iqa}

と求まる．よって，

\left(\begin{array}{} A_2\\ B_2 \end{array}\right) =\left(\begin{array}{} e^{iqa} & 0\\ 0 & e^{-iqa} \end{array}\right) \left(\begin{array}{} C_1\\ D_1 \end{array}\right)

と書くことができる．以上から

\displaystyle Y=\left(\begin{array}{} e^{iqa} & 0\\ 0 & e^{-iqa} \end{array}\right)

とわかった．

(8)

(7)における表式の\psi_2(x)と\psi_3(x)の関係についてこれらはどちらもx=aが中心の波であるから(5)と同様に考えよう．
まずB_2への寄与としては，A_2からの反射とD_2からの透過の寄与が考えられるので

B_2=rA_2+tD_2

となる．一方でC_2については，A_2からの透過とD_2からの反射が考えられるので

C_2=tA_2+rD_2

となる．これをC_2とD_2について解けばよいのだが，これは(5)の式と全く同一なので結果をそのまま流用すればよい．なので，

\left(\begin{array}{} C_2\\ D_2 \end{array}\right) =\left(\begin{array}{} \displaystyle t\left(1-\left(\frac{r}{t}\right)^2\right) & \displaystyle \frac{r}{t}\\ \displaystyle -\frac{r}{t} & \displaystyle \frac{1}{t} \end{array}\right) \left(\begin{array}{} A_2\\ B_2 \end{array}\right) =Z\left(\begin{array}{} A_2\\ B_2 \end{array}\right)

と書ける．これと(6),(7)での結果を踏まえれば

\left(\begin{array}{} C_2\\ D_2 \end{array}\right) =ZYZ\left(\begin{array}{} A_1\\ B_1 \end{array}\right)

となる．以上からXは，

X=ZYZ

とわかる．

(9)

a=2n\pi/qであるからe^{\pm iqa}=e^{\pm 2n\pi}=1となるのでYは単位行列に他ならない．よってXを求めるためにはZ^2を求めればよい．そこでまずZを\omegaを用いて表すことにしよう．(3)の結果からそれぞれ

\displaystyle \frac{r}{t}=\frac{i(\omega+i)}{\omega^2+1}\frac{\omega^2+1}{\omega(\omega+i)}=\frac{i}{\omega}
\displaystyle t\left(1-\left(\frac{r}{t}\right)^2\right)=\frac{\omega(\omega+i)}{\omega^2+1}\left(1+\frac{1}{\omega^2}\right)=1+\frac{i}{\omega}
\displaystyle \frac{1}{t}=\frac{\omega-i}{\omega}=1-\frac{i}{\omega}

となる．よってZは，

Z=\left(\begin{array}{} \displaystyle 1+\frac{i}{\omega} & \displaystyle \frac{i}{\omega}\\ \displaystyle -\frac{i}{\omega} & \displaystyle 1-\frac{i}{\omega} \end{array}\right)

と表せる．よって，Z^2は，

Z^2=\left(\begin{array}{} \displaystyle 1+\frac{i}{\omega} & \displaystyle \frac{i}{\omega}\\ \displaystyle -\frac{i}{\omega} & \displaystyle 1-\frac{i}{\omega} \end{array}\right)\left(\begin{array}{} \displaystyle 1+\frac{i}{\omega} & \displaystyle \frac{i}{\omega}\\ \displaystyle -\frac{i}{\omega} & \displaystyle 1-\frac{i}{\omega} \end{array}\right)\\ =\left(\begin{array}{} \displaystyle \left(1+\frac{i}{\omega}\right)^2+\frac{1}{\omega^2} & \displaystyle \frac{i}{\omega}\left(1+\frac{i}{\omega}\right)+\frac{i}{\omega}\left(1-\frac{i}{\omega}\right)\\ \displaystyle -\frac{i}{\omega}\left(1+\frac{i}{\omega}\right)-\frac{i}{\omega}\left(1-\frac{i}{\omega}\right) & \displaystyle \frac{1}{\omega^2}+\left(1-\frac{i}{\omega}\right)^2 \end{array}\right)\\ =\left(\begin{array}{} \displaystyle \frac{\omega+2i}{\omega} & \displaystyle \frac{2i}{\omega}\\ \displaystyle -\frac{2i}{\omega} & \displaystyle \frac{\omega-2i}{\omega} \end{array}\right)

と計算できる．以上から(8)において求めた式でC_2=t_2，D_2=0，A_1=1，B_1=r_2を代入すると

\displaystyle \left(\begin{array}{} t_2\\ 0 \end{array}\right) =\left(\begin{array}{} \displaystyle \frac{\omega+2i}{\omega} & \displaystyle \frac{2i}{\omega}\\ \displaystyle -\frac{2i}{\omega} & \displaystyle \frac{\omega-2i}{\omega} \end{array}\right) \left(\begin{array}{} 1\\ r_2 \end{array}\right)\\ =\left(\begin{array}{} \displaystyle \frac{\omega+2i}{\omega}+\displaystyle \frac{2i}{\omega}r_2\\ \displaystyle -\frac{2i}{\omega}+\displaystyle \frac{\omega-2i}{\omega}r_2 \end{array}\right)

となる．ゆえにまず，

\displaystyle r_2=\frac{2i}{\omega}\frac{\omega}{\omega-2i}=\frac{2i}{\omega-2i}=\frac{2i(\omega+2i)}{\omega^2+4}

と求まる．そしてt_2についても

\displaystyle t_2= \frac{\omega+2i}{\omega}+\frac{2i}{\omega}r_2=\frac{\omega+2i}{\omega}+\frac{2i}{\omega}\frac{2i}{\omega-2i}\\ \displaystyle =\frac{\omega}{\omega-2i}=\frac{\omega(\omega+2i)}{\omega^2+4}

となる．よって，|r_2|^2，|t_2|^2は，

\displaystyle |r_2|^2=\frac{4}{\omega^2+4}
\displaystyle |t_2|^2=\frac{\omega^2}{\omega^2+4}

と計算できる．最後にこれを(3)における結果と比較すると，

\displaystyle |r_2|^2-|r|^2=\frac{4}{\omega^2+4}-\frac{1}{\omega^2+1}=\frac{3\omega^2}{(\omega^2+4)(\omega^2+1)}>0

となる．つまり，II(9)における反射率の方がI(3)の反射率よりも大きいことがわかる．これはポテンシャル障壁が1つから2つになった結果により入射粒子がより反射されやすくなったからであると考えられる．

問題4

I

(1)

Gibbs free energyの全微分を計算し，Helmholtz free energyの全微分を代入すると，

dG=dF+VdP+PdV=-SdT-PdV+\mu dN+VdP+PdV\\ =-SdT+VdP+\mu dN

となる．よって，\mbox{（あ）}=-S，\mbox{（い）}=V，\mbox{（う）}=\muである．

(2)

(1)で求めたGibbs free energyの全微分から

\displaystyle \left(\frac{\partial G}{\partial N}\right)_{T,P}=\mu(T,P,N)

である．\mu(T,P,N)は今NによらないのでNで積分すると，

G=\mu N

と，Gを\muとNの2つの熱力学変数の積で表せる．そしてこの式の全微分を取り，それと(1)において求めたGibbs free energyの全微分を等置すると，

dG=Nd\mu+\mu dN=-SdT+VdP+\mu dN\\ \rightarrow Nd\mu=-SdT+VdP\\ \rightarrow d\mu=-sdT+vdP

となり，Gibbs-Duhemの関係式を得る．ただしここで，s=S/N，v=V/Nとした．

II

(3)

分配関数Z(T,V,N)をまず計算しよう：

\displaystyle Z(T,V,N)=\frac{V^N}{h^{3N}N!}\int d^3\textbf{p}_1\cdots d^3\textbf{p}_N\exp\left(-\frac{1}{k_\mathrm{B}T}\sum_{i=1}^{N}\frac{\textbf{p}_i^2}{2m}\right)\\ =\displaystyle\frac{V^N}{h^{3N}N!}\left(\int_{-\infty}^{\infty}dp\,\exp\left(-\frac{1}{k_{\mathrm{B}}T}\frac{p^2}{2m}\right)\right)^{3N}\\ =\displaystyle\frac{V^N}{h^{3N}N!}\left(\sqrt{2\pi mk_{\mathrm{B}}T}\right)^{3N} =\displaystyle \frac{V^N}{N!}\left(\frac{\sqrt{2\pi mk_{\mathrm{B}}T}}{h}\right)^{3N}\\ \displaystyle=\frac{V^N}{\lambda^{3N} N!}

さて，Helmholtz free energyについてF=-k_{\mathrm{B}}T\log Zであるから，

F(T,V,N) =-k_{\mathrm{B}}T\log Z(T,V,N) =-k_{\mathrm{B}}T\left[N\log V-3N\log\lambda-\log N!\right]\\ \sim-Nk_{\mathrm{B}}T\left[\log V-3\log\lambda-\log N+1\right]\\ =\displaystyle-Nk_{\mathrm{B}}T\log\left(\frac{eV}{\lambda^3N}\right)

とHelmholtz free energyが求まった．

(4)

Helmholtz free energyの全微分から化学ポテンシャル\muについて

\displaystyle \mu(T,V,N)=\left(\frac{\partial F(T,V,N)}{\partial N}\right)_{T,V}=\frac{\partial}{\partial N}\left(-Nk_{\mathrm{B}}T\log\left(\frac{eV}{\lambda^3N}\right)\right)\\ \displaystyle =-k_{\mathrm{B}}T\log\left(\frac{eV}{\lambda^3N}\right)-Nk_{\mathrm{B}}T\frac{\lambda^3N}{eV}\left(-\frac{eV}{\lambda^{3}N^{2}}\right)\\ \displaystyle =-k_{\mathrm{B}}T\log\left(\frac{V}{\lambda^3N}\right)

となる．ここで今考えているものは理想気体であるから状態方程式PV=Nk_{\mathrm{B}}Tを用いて表式からVを消去すると，

\displaystyle \mu(T,P,N)=-k_{\mathrm{B}}T\log\left(\frac{k_{\mathrm{B}}T}{\lambda^3P}\right)

と求まる．

(5)

一粒子状態密度（以下DOS）を求めるために，まずエネルギーが\epsilon以下であるような状態の数N(\epsilon)を求めよう．
まず自由粒子のSchrödinger方程式

\displaystyle -\frac{\hbar^2}{2m}\nabla^2\phi=\epsilon\phi

の解は

\phi=e^{\pm i\textbf{k}\cdot\textbf{r}}

となる．ここでk=\sqrt{2m\epsilon/\hbar^2}と置いた．今，周期境界条件\phi(x+L,y,z)=\phi(x,y,z)などを課すことにすれば（Lは箱の一辺の長さ），

k_iL=2\pi n_i\hspace{10pt}(i=x,y,z)

が得られる．
エネルギーが\epsilon以下の状態の総数を求めるためには半径がk=\sqrt{2m\epsilon/\hbar^2}である3次元球の中にある状態の数を数えればよいのだが，\textbf{k}-空間において1つの状態が占める体積は周期境界条件より得られた条件式から(2\pi/L)^3=(2\pi)^3/Vであるので，球の体積をこれで割ればよい．つまり，N(\epsilon)は，

\displaystyle N(\epsilon)=\frac{4\pi k^3}{\displaystyle\frac{(2\pi)^3}{V}}=\frac{4\pi}{3}\left(\sqrt{\frac{2m\epsilon}{\hbar^2}}\right)^3\frac{V}{8\pi^3}=\frac{4\pi V}{3h^3}\left(2m\epsilon\right)^{\frac{3}{2}}

と求められる．よってDOSD(\epsilon)は，

\displaystyle D(\epsilon)=\frac{dN(\epsilon)}{d\epsilon}=\frac{2\pi V(2m)^{\frac{3}{2}}}{h^3}\epsilon^{\frac{1}{2}}

となる．以上より，

\displaystyle A=\frac{2\pi(2m)^{\frac{3}{2}}}{h^3},\hspace{15pt}a=\frac{1}{2}

とわかる．

(6)

Maxwellの関係式

\displaystyle \frac{\partial P}{\partial \mu}=\frac{\partial N}{\partial V}

より両辺をVで積分すれば圧力Pが求められる．
さて，粒子数の平均Nについて

\displaystyle N=\int_{0}^{\infty}d\epsilon\,D(\epsilon)\frac{1}{e^{(\epsilon-\mu)/k_{\mathrm{B}}T}-1}=\int_{0}^{\infty}d\epsilon\,\frac{VA\epsilon^a}{e^{(\epsilon-\mu)/k_{\mathrm{B}}T}-1}

で与えられたから，これをVで微分し上式の右辺に代入すると次のようになる．

\displaystyle \frac{\partial P}{\partial \mu}=\frac{\partial N}{\partial V}=\int_{0}^{\infty}d\epsilon\,\frac{A\epsilon^a}{e^{(\epsilon-\mu)/k_{\mathrm{B}}T}-1}

この両辺を\muで積分しよう．積分範囲は-\infty,\muにとるが，\mu\to-\inftyでP\to0であるから積分の下限については省略することにする．

\displaystyle P(\mu,T)=\int^{\mu}d\mu'\int_{0}^{\infty}d\epsilon\,\frac{A\epsilon^a}{e^{(\epsilon-\mu')/k_{\mathrm{B}}T}-1}

さて，この積分を実行するために次の積分公式を用いよう．

\displaystyle \int^{x} dx'\,\frac{1}{Ae^{-x'}-1}=-x-\log\left(Ae^{-x}-1\right)

ここでAは1以上の定数であり，x\leq0である．
上の積分において\mu'/k_{\mathrm{B}}T=uと置換して積分公式を用いると，

\displaystyle P(\mu,T)=k_{\mathrm{B}}T\int_{0}^{\infty}d\epsilon\int^{\mu/k_{\mathrm{B}}T}du\,\frac{A\epsilon^a}{e^{\epsilon/k_{\mathrm{B}}T}e^{-u}-1}\\ \displaystyle =k_{\mathrm{B}}T\int_{0}^{\infty}d\epsilon\,A\epsilon^{a}\left(-\frac{\mu}{k_{\mathrm{B}}T}-\log\left(e^{(\epsilon-\mu)/k_{\mathrm{B}}T}-1\right)\right)\\ \displaystyle =-k_{\mathrm{B}}T\int_{0}^{\infty}d\epsilon\,A\epsilon^{a}\left(\frac{\mu}{k_{\mathrm{B}}T}+\log\left(e^{(\epsilon-\mu)/k_{\mathrm{B}}T}-1\right)\right)

となり，圧力の表式を得ることができる．

(7)

Helmholtz free energyの全微分から圧力Pは，

\displaystyle P=-\frac{\partial F}{\partial V}

で求められるのであった．まず，風船の中の圧力P_{\mathrm{in}}について求めよう．風船内のHelmholtz free energyは，

F_{\mathrm{tot}}(V)=F(V)+\gamma V{2/3}

であったからP_{\mathrm{in}}は，

\displaystyle P_{\mathrm{in}}=-\frac{\partial F(V)}{\partial V}-\frac{2}{3}\gamma V^{-1/3}

と計算できる．一方，風船外部の圧力P_{\mathrm{ex}}については，風船内部の気体と風船外部の気体は同一のものであるからそのHelmholtz free energyは風船内の気体のHelmholtz free energyと同一であると考えられる．ゆえにP_{\mathrm{ex}}は，

\displaystyle P_{\mathrm{ex}}=-\frac{\partial F(V)}{\partial V}

となる．よって圧力差\Delta Pは，

\displaystyle \Delta P=P_{\mathrm{in}}-P_{\mathrm{ex}}=-\frac{2}{3}\gamma V^{-1/3}

である．よって，Bとbは，

\displaystyle B=-\frac{2}{3}\gamma,\hspace{10pt}b=-\frac{1}{3}

とわかる．

(8)

\mu_L=\mu_Gであるから，d\mu_L=d\mu_Gが成立する．ここにGibbs-Duhemの関係式から得た式

\begin{array}{} d\mu_L=v_LdP_L\\ d\mu_G=v_GdP_G \end{array}

を代入すると，

\displaystyle v_LdP_L=v_GdP_G\\ \displaystyle \rightarrow dP_G=\frac{v_L}{v_G}dP_L

という式を得ることができる．
さらに\Delta P=P_L-P_Gの左辺に(7)の結果を代入し，そのうえで両辺の全微分を取ると

dP_L-dP_G=BbV^{b-1}dV\\ \rightarrow dP_L=BbV^{b-1}dV+dP_G

という式を得ることができるので，これを先に得た式のdP_Lに代入すると，

\displaystyle dP_G=\frac{v_L}{v_G}dP_L=\frac{v_L}{v_G}\left(BbV^{b-1}dV+dP_G\right)\\ \displaystyle \rightarrow \frac{dP_G}{dV}=\frac{v_L}{v_G-v_L}BbV^{b-1}

となり，P_Gが満たす微分方程式が得られた．
最後に微分方程式をv_L/v_Gの一次のオーダーで近似する．x\ll1のとき(1+x)^n\sim 1+nxと近似できることを用いると，

\displaystyle \frac{dP_G}{dV}=\frac{v_L}{v_G-v_L}BbV^{b-1}=\frac{v_L}{v_G}\frac{1}{\displaystyle 1-\frac{v_L}{v_G}}BbV^{b-1}\\ \displaystyle \sim\frac{v_L}{v_G}\left(1+\frac{v_L}{v_G}\right)BbV^{b-1}\\ \displaystyle \sim \frac{v_L}{v_G}BbV^{b-1}

を得られる．

(9)

v_G=k_{\mathrm{B}}T/P_Gを(8)で得られた近似された微分方程式に代入すると，

\displaystyle \frac{dP_G}{dV}=\frac{v_L}{v_G}BbV^{b-1}=\frac{v_L}{k_{\mathrm{B}}T}P_GBbV^{b-1}\\ \displaystyle \sim\frac{dP_G}{P_G}=\frac{v_L}{k_{\mathrm{B}}T}BbV^{b-1}dV

となる．両辺を積分すると，

\displaystyle \log P_G=\frac{v_L}{k_{\mathrm{B}}T}BV^{b}+C\hspace{10pt}(C:\mathrm{arbitrary\hspace{2pt}constant})\\ \displaystyle P_G=A\exp\left(\frac{v_L}{k_{\mathrm{B}}T}BV^{b}\right)\hspace{10pt}(A:\mathrm{arbitrary\hspace{2pt}constant})

を得る．ここで(7)よりB=-2/3\gamma，b=-1/3であるからこれを代入すると，

\displaystyle P_G(V)=A\exp\left(-\frac{2}{3}\frac{v_L\gamma}{k_{\mathrm{B}}T}V^{\displaystyle-\frac{1}{3}}\right)

と解を求めることができた．最後に積分定数Aを境界条件から定めよう．与えられた境界条件はP_{\infty}=\lim_{V\to\infty}P_G(V)であるから，

\displaystyle P_{\infty}=\lim_{V\to\infty}P_G(V)=\lim_{V\to\infty}A\exp\left(-\frac{2}{3}\frac{v_L\gamma}{k_{\mathrm{B}}T}V^{\displaystyle-\frac{1}{3}}\right) =A\exp\left(-\frac{2}{3}\frac{v_L\gamma}{k_{\mathrm{B}}T}\cdot0\right)=A

となり，A=P_{\infty}とわかる．以上よりP_G(V)は，

\displaystyle P_G(V)=P_{\infty}\exp\left(-\frac{2}{3}\frac{v_L\gamma}{k_{\mathrm{B}}T}V^{\displaystyle-\frac{1}{3}}\right)

と求められた．

二重ポテンシャル障壁による共鳴透過

量子力学

一次元量子系のポテンシャル障壁による散乱において、E>V_0のとき、Lを障壁の長さとすると
\displaystyle \frac{\sqrt{2m(E-V_0)}}{\hbar}L=n\pi
を満たす場合、透過率TがT=1となり障壁を反射されることなく完全に通り抜けることが知られていて、これは共鳴透過と呼ばれている。
ポテンシャル障壁が1つしかないときは、E>V_0の場合にしか共鳴透過は起こらないのだが、ポテンシャル障壁が2つある場合にはE\lt V_0の場合にもこの共鳴透過が起こる*1。このことを実際に散乱問題を解くことで確認してみよう。

問題設定

f:id:MiyanTarumi:20210502205910j:plain
二重ポテンシャル障壁
上のようなポテンシャルの散乱問題を考えよう。そして、z\lt 0における波動関数を\varphi_{\mathrm{inc}}、0\lt z\lt \lt Lの波動関数を\varphi_{\mathrm{I}}、L\lt z\lt L+\ellの波動関数を\varphi_{\mathrm{II}}、L+\ell\lt z\lt 2L+\ellの波動関数を\varphi_{\mathrm{III}}、2L+\ell\lt zの波動関数を\varphi_{\mathrm{sc}}と置くことにする。

Schrödinger方程式とその接続条件

各区間におけるSchrödinger方程式を並べてやると次のようになる。

\displaystyle -\frac{\hbar^2}{2m} \frac{d^2\varphi_{\mathrm{inc}}}{dx^2}=E\varphi_{\mathrm{inc}}
\displaystyle \left(-\frac{\hbar^2}{2m} \frac{d^2}{dx^2}+V\_0\right)\varphi_{\mathrm{I}}=E\varphi_{\mathrm{I}}
\displaystyle -\frac{\hbar^2}{2m} \frac{d^2\varphi_{\mathrm{II}}}{dx^2}=E\varphi_{\mathrm{II}}
\displaystyle \left(-\frac{\hbar^2}{2m} \frac{d^2}{dx^2}+V\_0\right)\varphi_{\mathrm{III}}=E\varphi_{\mathrm{III}}
\displaystyle -\frac{\hbar^2}{2m} \frac{d^2\varphi_{\mathrm{sc}}}{dx^2}=E\varphi_{\mathrm{sc}}

これらを解けばA,B,C,D,E,F,G,H,I,Jを積分定数として、
\varphi_{\mathrm{inc}}=Ae^{ikz}+Be^{-ikz}

\varphi_{\mathrm{I}}=Ce^{i\kappa z}+De^{-i\kappa z}

\varphi_{\mathrm{II}}=Ee^{ik(z-L)}+Fe^{-ik(z-L)}

\varphi_{\mathrm{III}}=Ge^{ik(z-L-\ell)}+He^{-ik(z-L-\ell)}

\varphi_{\mathrm{sc}}=Ie^{ik(z-2L-\ell)}+Je^{-ik(z-2L-\ell)}
ここで、
\displaystyle k=\frac{\sqrt{2mE}}{\hbar},\hspace{10pt}\kappa=\frac{\sqrt{2m(V_0-E)}}{\hbar}
である。
接続条件は、ポテンシャルが有限なことを踏まえると各波動関数とその微分が各区間端において連続であればよい。ゆえに、
A+B=C+D,\hspace{10pt}ik(A-B)=\kappa(C-D)

Ce^{\kappa L}+De^{-\kappa L}=E+F,\hspace{10pt}\kappa(Ce^{\kappa L}-De^{-\kappa L})=ik(E-F)

Ee^{ik\ell}+Fe^{-ik\ell}=G+H,\hspace{10pt}ik(Ee^{ik\ell}-Fe^{-ik\ell})=\kappa(G-H)

Ge^{\kappa L}+He^{-\kappa L}=I+J,\hspace{10pt}\kappa(Ge^{\kappa L}-He^{-\kappa L})=ik(I-J)
であり、これらの方程式から積分定数を決めればよい。しかし、これを強引に解くのはかなり骨が折れるので行列を用いて解くことにする。そのために、それぞれを整理して次のように書いてみよう。面倒なので\displaystyle\beta=\frac{\kappa}{ik}と書くことにして、

\begin{pmatrix} 1&1\\ 1&-1 \end{pmatrix} \begin{pmatrix} C\\ D \end{pmatrix} = \begin{pmatrix} 1&1\\ \frac{1}{\beta}&-\frac{1}{\beta} \end{pmatrix} \begin{pmatrix} A\\ B \end{pmatrix}
\begin{pmatrix} 1&1\\ 1&-1 \end{pmatrix} \begin{pmatrix} E\\ F \end{pmatrix} = \begin{pmatrix} e^{\kappa L}&e^{-\kappa L}\\ \beta e^{\kappa L}&-\beta e^{-\kappa L} \end{pmatrix} \begin{pmatrix} C\\ D \end{pmatrix}
\begin{pmatrix} 1&1\\ 1&-1 \end{pmatrix} \begin{pmatrix} G\\ H \end{pmatrix} = \begin{pmatrix} e^{ik\ell}&e^{-ik\ell}\\ \frac{1}{\beta}e^{ik\ell}&-\frac{1}{\beta}e^{-ik\ell} \end{pmatrix} \begin{pmatrix} E\\ F \end{pmatrix}
\begin{pmatrix} 1&1\\ 1&-1 \end{pmatrix} \begin{pmatrix} I\\ J \end{pmatrix} = \begin{pmatrix} e^{\kappa L}&e^{-\kappa L}\\ \beta e^{\kappa L}&-\beta e^{-\kappa L} \end{pmatrix} \begin{pmatrix} G\\ H \end{pmatrix}

と書き直せる。積分定数を求めるにはこれらの連立方程式を解けばよい。

積分定数の決定1

さて、上で求めた連立方程式を解くためにまず準備をいくつか行っておく。
行列について

E_2=\begin{pmatrix} 1&1\\ 1&-1 \end{pmatrix}
M(x,y)=\begin{pmatrix} e^{x}&e^{-x}\\ ye^{x}&-ye^{-x} \end{pmatrix}

とE_2とM(x,y)を定義しよう。
また、連立方程式を解くためにE_2の逆行列を求めると、その逆行列E_2^{-1}は、

E_2^{-1}=\displaystyle\frac{1}{2} \begin{pmatrix} 1&1\\ 1&-1 \end{pmatrix} =\frac{1}{2}E_2

となっている。
以上の下で上の連立方程式を改めて書き直すと、

\begin{pmatrix} C\\ D \end{pmatrix} =\displaystyle\frac{1}{2}E_2M\left(0,\frac{1}{\beta}\right) \begin{pmatrix} A\\ B \end{pmatrix}
\begin{pmatrix} E\\ F \end{pmatrix} =\displaystyle\frac{1}{2}E_2M\left(\kappa L,\beta\right) \begin{pmatrix} C\\ D \end{pmatrix}
\begin{pmatrix} G\\ H \end{pmatrix} =\displaystyle\frac{1}{2}E_2M\left(ik\ell,\frac{1}{\beta}\right) \begin{pmatrix} E\\ F \end{pmatrix}
\begin{pmatrix} I\\ J \end{pmatrix} =\displaystyle\frac{1}{2}E_2M\left(\kappa L,\beta\right) \begin{pmatrix} G\\ H \end{pmatrix}

となる。ここで、共通して出てくるのはE_2M(x,y)という行列なのでこれをN(x,y)と書くことにしよう。そしてこれを具体的に計算してみると、

N_2(x,y)=E_2M(x,y)= \begin{pmatrix} 1&1\\ 1&-1 \end{pmatrix} \begin{pmatrix} e^{x}&e^{-x}\\ ye^{x}&-ye^{-x} \end{pmatrix} = \begin{pmatrix} (1+y)e^{x}&(1-y)e^{-x}\\ (1-y)e^{x}&(1+y)e^{-x} \end{pmatrix}

となる。これを用いて連立方程式を解いてみる。上の連立方程式を順に代入していくと、

\begin{pmatrix} I\\ J \end{pmatrix} =\displaystyle\frac{1}{16}N(\kappa L,\beta)N\left(ik\ell,\frac{1}{\beta}\right)N(\kappa L,\beta)N\left(0,\frac{1}{\beta}\right) \begin{pmatrix} A\\ B \end{pmatrix}

となり、行列の計算を行えばこの連立方程式を解けたことになる。
具体的に行列を計算していけばよいのだが、少々面倒なのでかけていく行列の規則性に注目してみよう。上の式を見ると、N(x,y)はy=\betaとy=\frac{1}{\beta}のペアの繰り返しで出てきている。そこでこの二つをまとめてT\left(x,\beta;x',\frac{1}{\beta}\right)=N(x,\beta)N\left(x',\frac{1}{\beta}\right)と表すことにし、これを転送行列と呼ぶことにする。なお、連立方程式を見直せばこの転送行列一つ一つ一つがポテンシャル障壁を一つ通り抜けることに対応していることがわかるだろう。
転送行列の成分を具体的に計算してみよう。定義から

T\left(x,\beta;x',\frac{1}{\beta}\right)=N(x,\beta)N\left(x',\frac{1}{\beta}\right)\\ =\begin{pmatrix} (1+\beta)e^{x}&(1-\beta)e^{-x}\\ (1-\beta)e^{x}&(1+\beta)e^{-x} \end{pmatrix} \begin{pmatrix} \left(1+\frac{1}{\beta}\right)e^{x'}&\left(1-\frac{1}{\beta}\right)e^{-x'}\\ \left(1-\frac{1}{\beta}\right)e^{x'}&\left(1+\frac{1}{\beta}\right)e^{-x'} \end{pmatrix}\\ =\begin{pmatrix} e^{x'}\left((1+\beta)\left(1+\frac{1}{\beta}\right)e^{x}+(1-\beta)\left(1-\frac{1}{\beta}\right)e^{-x}\right)&e^{-x'}\left((1+\beta)\left(1-\frac{1}{\beta}\right)e^{x}+(1-\beta)\left(1+\frac{1}{\beta}\right)e^{-x}\right)\\ e^{x'}\left((1-\beta)\left(1+\frac{1}{\beta}\right)e^{x}+(1+\beta)\left(1-\frac{1}{\beta}\right)e^{-x}\right)&e^{-x'}\left((1-\beta)\left(1-\frac{1}{\beta}\right)e^{x}+(1+\beta)\left(1+\frac{1}{\beta}\right)e^{-x}\right) \end{pmatrix}

とわかる。
今の場合の転送行列2つにおいてx=\kappa Lは共通しているので、T\left(x,\beta;x',\frac{1}{\beta}\right)の成分のうちxのみによる部分を文字で置くことにしよう。

a(x)=(1+\beta)\left(1+\frac{1}{\beta}\right)e^{x}+(1-\beta)\left(1-\frac{1}{\beta}\right)e^{-x}\\ b(x)=(1+\beta)\left(1-\frac{1}{\beta}\right)e^{x}+(1-\beta)\left(1+\frac{1}{\beta}\right)e^{-x}\\ c(x)=-b(x)=(1-\beta)\left(1+\frac{1}{\beta}\right)e^{x}+(1+\beta)\left(1-\frac{1}{\beta}\right)e^{-x}\\ d(x)=(1-\beta)\left(1-\frac{1}{\beta}\right)e^{x}+(1+\beta)\left(1+\frac{1}{\beta}\right)e^{-x}

このように置いたうえで連立方程式の解をこれらを用いて表すと次のようになる。

\begin{pmatrix} I\\ J \end{pmatrix} =\displaystyle\frac{1}{16}T\left(\kappa L,\beta;ik\ell,\frac{1}{\beta}\right)T\left(\kappa L,\beta;0,\frac{1}{\beta}\right) \begin{pmatrix} A\\ B \end{pmatrix}\\ =\displaystyle\frac{1}{16} \begin{pmatrix} e^{ik\ell}a(\kappa L)&e^{-ik\ell}b(\kappa L)\\ e^{ik\ell}c(\kappa L)&e^{-ik\ell}d(\kappa L) \end{pmatrix} \begin{pmatrix} a(\kappa L)&b(\kappa L)\\ c(\kappa L)&d(\kappa L) \end{pmatrix} \begin{pmatrix} A\\ B \end{pmatrix}\\ =\displaystyle\frac{1}{16}\begin{pmatrix} e^{ik\ell}a^2+e^{-ik\ell}bc&e^{ik\ell}ab+e^{-ik\ell}bd\\ e^{ik\ell}ac+e^{-ik\ell}cd&e^{ik\ell}bc+e^{-ik\ell}d^2 \end{pmatrix} \begin{pmatrix} A\\ B \end{pmatrix}

ここで最後の表式についてa(\kappa L)などの引数について省略した。
以上から、連立方程式を完全に解くことができてA,B,I,Jは、

I=\displaystyle\frac{1}{16}\left(\left(e^{ik\ell}a^2+e^{-ik\ell}bc\right)A+\left(e^{ik\ell}ab+e^{-ik\ell}bd\right)B\right)\\ J=\displaystyle\frac{1}{16}\left(\left(e^{ik\ell}ac+e^{-ik\ell}cd\right)A+\left(e^{ik\ell}bc+e^{-ik\ell}d^2\right)B\right)

とわかる。
ここで物理的な要請から変数を減らす。粒子はx=-\inftyから入射してくるとしよう。すると2L+\ell\lt zの領域では反射波がないのでJ=0と置ける。また、簡単のためA=1としよう。これらから残りの変数について完全に決定することができる。
まず、Bについて２つ目の式から

B=-\displaystyle \frac{c\left(e^{ik\ell}a+e^{-ik\ell}d\right)}{e^{ik\ell}bc+e^{-ik\ell}d^2}

となる。そして、Iについても

I=\displaystyle\frac{1}{16}\left(\left(e^{ik\ell}a^2+e^{-ik\ell}bc\right)+\left(e^{ik\ell}ab+e^{-ik\ell}bd\right)B\right)\\ =\displaystyle\frac{1}{16}\left(\left(e^{ik\ell}a^2+e^{-ik\ell}bc\right)+\left(e^{ik\ell}ab+e^{-ik\ell}bd\right)\left(-\frac{c\left(e^{ik\ell}a+e^{-ik\ell}d\right)}{e^{ik\ell}bc+e^{-ik\ell}d^2}\right)\right)\\ =\displaystyle\frac{1}{16}\left(\frac{\left(e^{ik\ell}a^2+e^{-ik\ell}bc\right)\left(e^{ik\ell}bc+e^{-ik\ell}d^2\right)-bc\left(e^{ik\ell}a+e^{-ik\ell}d\right)^2}{e^{ik\ell}bc+e^{-ik\ell}d^2}\right)\\ =\displaystyle\frac{1}{16}\frac{a^2d^2+b^2c^2-2adbc}{e^{ik\ell}bc+e^{-ik\ell}d^2}\\ =\displaystyle\frac{1}{16}\frac{\left(ad-bc\right)^2}{e^{ik\ell}bc+e^{-ik\ell}d^2}

と計算できる。

積分定数の決定2

あとはa,b,c,dを代入すればよいのだが、先にa,b,c,dについてある程度式変形を行っておこう：

\displaystyle a(\kappa L)=(1+\beta)\left(1+\frac{1}{\beta}\right)e^{\kappa L}+(1-\beta)\left(1-\frac{1}{\beta}\right)e^{-\kappa L}\\ \displaystyle=2\left(e^{\kappa L}+e^{-\kappa L}\right)+\left(\beta+\frac{1}{\beta}\right)\left(e^{\kappa L}-e^{-\kappa L}\right)\\ \displaystyle=2\left(2\cosh(\kappa L)+\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\right)
\displaystyle b(\kappa L)=(1+\beta)\left(1-\frac{1}{\beta}\right)e^{\kappa L}+(1-\beta)\left(1+\frac{1}{\beta}\right)e^{-\kappa L}\\ \displaystyle=\left(\beta-\frac{1}{\beta}\right)\left(e^{\kappa L}-e^{-\kappa L}\right)\\ \displaystyle=2\left(\beta-\frac{1}{\beta}\right)\sinh(\kappa L)
\displaystyle c(\kappa L)=-b(x)=-2\left(\beta-\frac{1}{\beta}\right)\sinh(\kappa L)
\displaystyle d(\kappa L)=(1-\beta)\left(1-\frac{1}{\beta}\right)e^{\kappa L}+(1+\beta)\left(1+\frac{1}{\beta}\right)e^{-\kappa L}\\ \displaystyle=2\left(e^{\kappa L}+e^{-\kappa L}\right)-\left(\beta+\frac{1}{\beta}\right)\left(e^{\kappa L}-e^{-\kappa L}\right)\\ \displaystyle=2\left(2\cosh(\kappa L)-\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\right)

これらを踏まえたうえで計算を行おう。
まず、Bの分子の部分について計算を行う。

-c\left(e^{ik\ell}a+e^{-ik\ell}d\right)\\ =\displaystyle2\left(\beta-\frac{1}{\beta}\right)\sinh(\kappa L)\left[e^{ik\ell}2\left(2\cosh(\kappa L)+\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\right)\right.\\ \displaystyle\hspace{10pt}\left.+e^{-ik\ell}2\left(2\cosh(\kappa L)-\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\right)\right]\\ \displaystyle=2\left(\beta-\frac{1}{\beta}\right)\sinh(\kappa L)\left[4\cosh(\kappa L)\left(e^{ik\ell}+e^{-ik\ell}\right)+2\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\left(e^{ik\ell}-e^{-ik\ell}\right)\right]\\ \displaystyle=16\left(\beta-\frac{1}{\beta}\right)\sinh(\kappa L)\left[\cosh(\kappa L)\cos(k\ell)+\frac{i}{2}\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\sin(k\ell)\right]\\ \displaystyle=16\left(\beta-\frac{1}{\beta}\right)\sinh(\kappa L)R\cos(k\ell+\phi)

ここで

\displaystyle R^2=\cosh^2(\kappa L)-\frac{1}{4}\left(\beta+\frac{1}{\beta}\right)^2\sinh^2(\kappa L)\\ \displaystyle\cos(\phi)=\frac{\cosh(\kappa L)}{R}\\ \displaystyle\sin(\phi)=-\frac{i}{2}\frac{\displaystyle\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)}{R}

である。
次にIの分子の部分の計算を行う。まず、

ad-bc\\ \displaystyle=\left(2\left(2\cosh(\kappa L)+\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\right)\right)\left(2\left(2\cosh(\kappa L)-\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\right)\right)\\ \displaystyle\hspace{10pt}+4\left(\beta-\frac{1}{\beta}\right)^2\sinh^2(\kappa L)\\ \displaystyle=4\left(4\cosh^2(\kappa L)-\left(\beta+\frac{1}{\beta}\right)^2\sinh(\kappa L)^2+\left(\beta-\frac{1}{\beta}\right)^2\sinh^2(\kappa L)\right)\\ \displaystyle=4\left(4\cosh^2(\kappa L)-4\sinh^2(\kappa L)\right)\\ =16

である。なお最後の式変形において\cosh^2(x)-\sinh^2(x)=1を用いた。
よって、

(ad-bc)^2=16^2

とわかった。
最後にBとIに共通する分母の部分の計算を行おう。

e^{ik\ell}bc+e^{-ik\ell}d^2\\ \displaystyle=-4e^{ik\ell}\left(\beta-\frac{1}{\beta}\right)^2\sinh^2(\kappa L)+4e^{-ik\ell}\left(2\cosh(\kappa L)-\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\right)^2\\ \displaystyle=-4e^{ik\ell}\left(\left(\beta+\frac{1}{\beta}\right)^2-4\right)\sinh^2(\kappa L)+4e^{-ik\ell}\left(2\cosh(\kappa L)-\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\right)^2\\ \displaystyle=4e^{-ik\ell}\left(4\cosh^2(\kappa L)-4\left(\beta+\frac{1}{\beta}\right)\cosh(\kappa L)\sinh(\kappa L)\right)\\ \displaystyle\hspace{10pt}-4\sinh^2(\kappa L)\left(2i\left(\beta+\frac{1}{\beta}\right)^2\sin(kl)-4e^{ik\ell}\right)\\ \displaystyle=16\cos(k\ell)\left(\cosh^2(\kappa L)-\left(\beta+\frac{1}{\beta}\right)\cosh(\kappa L)\sinh(\kappa L)\right)\\ \displaystyle\hspace{10pt}-16i\sin(k\ell)\left(\cosh^2(\kappa L)-\left(\beta+\frac{1}{\beta}\right)\cosh(\kappa L)\sinh(\kappa L)\right)\\ \displaystyle\hspace{10pt}-8i\left(\beta+\frac{1}{\beta}\right)^2\sin(kl)\sinh^2(\kappa L)+16\left(\cos(k\ell)+i\sin(k\ell)\right)\sinh^2(\kappa L)\\ \displaystyle=16\cos(k\ell)\left(\cosh^2(\kappa L)-\left(\beta+\frac{1}{\beta}\right)\cosh(\kappa L)\sinh(\kappa L)\right)\\ \displaystyle\hspace{10pt}+16i\sin(k\ell)\left(\beta+\frac{1}{\beta}\right)\cosh(\kappa L)\sinh(\kappa L)-16i\sin(k\ell)\\ \displaystyle\hspace{10pt}-8i\left(\beta+\frac{1}{\beta}\right)^2\sin(kl)\sinh^2(\kappa L)+16\cos(k\ell)\sinh^2(\kappa L)\\

ここで、\betaがその定義から純虚数であったことを思い出そう。このことに注意したうえで実部と虚部に上の式を分けると、実部について

\mathrm{Re}\left[e^{ik\ell}bc+e^{-ik\ell}d^2\right]\\ \displaystyle=16\cos(k\ell)\cosh^2(\kappa L)+16i\sin(k\ell)\left(\beta+\frac{1}{\beta}\right)\cosh(\kappa L)\sinh(\kappa L)+16\cos(k\ell)\sinh^2(\kappa L)\\ \displaystyle=32\cos(k\ell)\cosh^2(\kappa L)+16i\sin(k\ell)\left(\beta+\frac{1}{\beta}\right)\cosh(\kappa L)\sinh(\kappa L)-16\cos(k\ell)\\ \displaystyle=32\cosh(\kappa L)\left(\cosh(\kappa L)\cos(k\ell)+\frac{i}{2}\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\sin(k\ell)\right)-16\cos(k\ell)\\ \displaystyle=16\left[2R\cos(k\ell+\phi)\cosh(\kappa L)-\cos(k\ell)\right]

と整理できる。
続いて虚部については

\mathrm{Im}\left[e^{ik\ell}bc+e^{-ik\ell}d^2\right]\\ \displaystyle=16i\cos(k\ell)\left(\beta+\frac{1}{\beta}\right)\cosh(\kappa L)\sinh(\kappa L)-8\left(\beta+\frac{1}{\beta}\right)^2\sin(kl)\sinh^2(\kappa L)-16\sin(k\ell)\\ \displaystyle=16i\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\left(\cosh(\kappa L)\cos(k\ell)+\frac{i}{2}\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\sin(k\ell)\right)-16\sin(k\ell)\\ \displaystyle=16\left[i\left(\beta+\frac{1}{\beta}\right)R\cos(k\ell+\phi)\sinh(\kappa L)-\sin(k\ell)\right]

となる。よって、分母の部分は、

e^{ik\ell}bc+e^{-ik\ell}d^2\\ \displaystyle=16\left[2R\cos(k\ell+\phi)\cosh(\kappa L)-\cos(k\ell)\right]+16i\left[i\left(\beta+\frac{1}{\beta}\right)R\cos(k\ell+\phi)\sinh(\kappa L)-\sin(k\ell)\right]

積分定数の決定3

以上からBとIについて

\displaystyle B=\frac{\displaystyle\left(\beta-\frac{1}{\beta}\right)\sinh(\kappa L)R\cos(k\ell+\phi)}{\displaystyle\left[2R\cos(k\ell+\phi)\cosh(\kappa L)-\cos(k\ell)\right]+i\left[i\left(\beta+\frac{1}{\beta}\right)R\cos(k\ell+\phi)\sinh(\kappa L)-\sin(k\ell)\right]}
\displaystyle I=\frac{\displaystyle 1}{\displaystyle\left[2R\cos(k\ell+\phi)\cosh(\kappa L)-\cos(k\ell)\right]+i\left[i\left(\beta+\frac{1}{\beta}\right)R\cos(k\ell+\phi)\sinh(\kappa L)-\sin(k\ell)\right]}

と求められた。

反射率Rと透過率T

今入射波の絶対値を1と置いているので反射率Rと透過率Tはそれぞれ係数の絶対値二乗で与えられることに注意しよう。
面倒なのは分母の部分なので分母のみを先に計算しておく。

\displaystyle\left|\left[2R\cos(k\ell+\phi)\cosh(\kappa L)-\cos(k\ell)\right]+i\left[i\left(\beta+\frac{1}{\beta}\right)R\cos(k\ell+\phi)\sinh(\kappa L)-\sin(k\ell)\right]\right|^2\\ \displaystyle=\left[2R\cos(k\ell+\phi)\cosh(\kappa L)-\cos(k\ell)\right]^2+\left[i\left(\beta+\frac{1}{\beta}\right)R\cos(k\ell+\phi)\sinh(\kappa L)-\sin(k\ell)\right]^2\\ \displaystyle=4R^2\cos^2(k\ell+\phi)\cosh^2(\kappa L)-4R\cos(k\ell+\phi)\cosh(\kappa L)\cos(k\ell)+\cos^2(k\ell)\\ \displaystyle\hspace{10pt}-\left(\beta+\frac{1}{\beta}\right)^2R^2\cos^2(k\ell+\phi)\sinh^2(\kappa L)-2i\left(\beta+\frac{1}{\beta}\right)R\cos(k\ell+\phi)\sinh(\kappa L)\sin(k\ell)+\sin^2(k\ell)\\ \displaystyle=1+R^2\cos^2(k\ell+\phi)\left(4\cosh^2(\kappa L)-\left(\beta+\frac{1}{\beta}\right)^2\sinh^2(\kappa L)\right)\\ \displaystyle\hspace{10pt}-4R\cos(k\ell+\phi)\left(\cosh(\kappa L)\cos(k\ell)+\frac{i}{2}\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\sin(k\ell)\right)\\ \displaystyle=1+R^2\cos^2(k\ell+\phi)\left(4\cosh^2(\kappa L)-\left(\beta+\frac{1}{\beta}\right)^2\sinh^2(\kappa L)-4\right)\\ \displaystyle=1+R^2\cos^2(k\ell+\phi)\left(4-\left(\beta+\frac{1}{\beta}\right)^2\right)\sinh^2(\kappa L)

となる。
よって反射率Rと透過率Tは、\betaが純虚数なことに注意して

R=|B|^2\\ \displaystyle\hspace{10pt}=\frac{\displaystyle-\left(\beta-\frac{1}{\beta}\right)^2\sinh^2(\kappa L)R^2\cos^2(k\ell+\phi)}{\displaystyle1+R^2\cos^2(k\ell+\phi)\left(4-\left(\beta+\frac{1}{\beta}\right)^2\right)\sinh^2(\kappa L)}
T=|I|^2\\ \displaystyle\hspace{10pt}=\frac{\displaystyle1}{\displaystyle1+R^2\cos^2(k\ell+\phi)\left(4-\left(\beta+\frac{1}{\beta}\right)^2\right)\sinh^2(\kappa L)}

であることがわかった。これは明らかにR+T=1を満たす。
こうして得られた透過率Tについて、\cos(k\ell+\phi)=0であるときにT=1、つまりエネルギーがポテンシャルよりも小さいのにも関わらず粒子が完全に透過してしまうことがわかるだろう。これが二重ポテンシャル障壁における共鳴透過と呼ばれる現象であり、江崎玲於奈によるトンネル電流の最も簡単なモデルである。
共鳴透過が起こる条件を調べよう。\cos(k\ell+\phi)の定義より、\cos(k\ell+\phi)=0は次のように書ける。

\displaystyle\cosh(\kappa L)\cos(k\ell)+\frac{i}{2}\left(\beta+\frac{1}{\beta}\right)\sinh(\kappa L)\sin(k\ell)=0\\ \displaystyle\Leftrightarrow \frac{2i}{\displaystyle\left(\beta+\frac{1}{\beta}\right)}=\tanh(\kappa L)\tan(k\ell)

そして、左辺について\beta=\frac{\kappa}{ik}、\kappa=\sqrt{\frac{2m(V_0-E)}{\hbar^2}}、k=\sqrt{\frac{2mE}{\hbar^2}}を代入して計算すると、

\displaystyle\frac{2i}{\displaystyle\left(\beta+\frac{1}{\beta}\right)}=\frac{2i}{\displaystyle\frac{\kappa}{ik}+\frac{ik}{\kappa}}=\frac{-2k\kappa}{\kappa^2-k^2}\\ \displaystyle=\frac{-2\sqrt{E(V_0-E)}}{V_0-2E}

となるので、結局、共鳴透過の条件は、

\displaystyle\frac{2\sqrt{E(V_0-E)}}{2E-V_0}=\tanh(\kappa L)\tan(k\ell)

となる。
そこで、\hbar=1の単位系において、横軸に入射エネルギーEを取り、V_0=10、m=\frac{1}{2}、L=1、\ell=5としてプロットしたのが以下の図である。

f:id:MiyanTarumi:20210504025020p:plain
透過率
上の図において、赤線は共鳴透過の条件の左辺、青線は条件式の右辺であり、これらの交点が共鳴透過が起こるようなエネルギーに対応している。そして、緑線は透過率Tの対数を取ったものであり、\log T=0である点で共鳴透過が起こっている。

共鳴透過の条件と井戸型ポテンシャルの固有状態

なぜこのような不思議な現象が起きているかを考えてみよう。
そのヒントは、二重ポテンシャル障壁において障壁と障壁の間が有限井戸型ポテンシャルのようになっていることにある。
有限井戸型ポテンシャルの固有値を決める方程式は、

\displaystyle \frac{2k\kappa}{\kappa^2-k^2}=\tan(k\ell)\Leftrightarrow \frac{2\sqrt{E(V_0-E)}}{2E-V_0}=\tan(k\ell)

であった。なお、ここで用いた井戸型ポテンシャルのモデルは一般に使われるような軸対象の問題設定とは少し違うことに注意せよ。具体的には、井戸の始まりはz=Lであり、そこから[L,L+\ell]まで井戸の底が続き、井戸の終わりがz=L+\ellにあるような問題設定である。
さて、これと先の共鳴透過の条件式を見比べると、\tanh(\kappa L)が右辺にかかっているかどうかの違いしかないことがわかる。\tanh(x)がx\rightarrow\inftyで1になることを考えれば、有限井戸型ポテンシャルは二重ポテンシャル障壁においてV_0\rightarrow\inftyとしたものであり、そしてこれは数値計算すればわかることなのであるが、共鳴透過が起こるエネルギーと有限井戸型ポテンシャルのエネルギー固有値は、右辺に\tanh(\kappa L)がかかるかどうかの違いしかないのでほぼ一致する。つまり、入射粒子のエネルギーがポテンシャル障壁のはざまによって作られた疑似的な井戸型ポテンシャルが作っているであろうエネルギー固有状態のエネルギーと一致しているときに共鳴が起こり、入射粒子が減衰されることなく透過するのである。

終わりに

二重ポテンシャル障壁において共鳴透過が起こることの原因は、ポテンシャル障壁の間にできた疑似的な井戸型ポテンシャルの中でエネルギー固有状態のようなものが形成されていることであったことを考えると、ポテンシャル障壁の数を3つ、4つと増やしていっても同じような現象が起きそうである。けど、転送行列の計算がめっちゃ面倒そうなので、こういう行列が簡単にできる方法があれば誰か教えてください。
量子力学の演義でやらされるポテンシャル障壁の散乱にもう一つ障壁を付け加えるだけで変な現象が起きたりするのは結構面白いとおもう。その分計算も結構面倒になるけど。。。

*1:直観的にはかなり気持ち悪い現象だと思う。

統計力学における形式的なアンサンブルの作り方

統計力学

統計力学においてカノニカルアンサンブルやグランドカノニカルアンサンブル、あるいはT \mbox{-} pアンサンブルなどを構成する方法を紹介する。
なお、この記事は統計力学のはじめの一歩 2020年版４章などの内容をまとめ直したものである。

統計力学の復習

Boltzmann原理によれば系のマクロなエントロピーS(E)とエネルギーが[E,E+\Delta E]の間にあるような状態の数W(E,\Delta E)の間には

S(E)=k_{\mathrm{B}}\log W(E,\Delta E)

の関係があるのであった。
またこのBoltzman原理は状態密度\Omega(E)を用いれば

S(E)=k_{\mathrm{B}}\log\Omega(E)\Delta E

と書くこともできる。
ミクロカノニカルアンサンブルから話を始めよう。
ミクロカノニカルアンサンブルにおいて、エネルギーがE,E+\Delta Eの間にある状態である確率は、

P(E)\Delta E=\displaystyle{\frac{1}{W(E,\Delta E)}}

となるのであった。
そして、導出は省略する*1が温度Tの熱浴と対象系を全系としてミクロカノニカルアンサンブルを適用してやれば、対象系のエネルギーがEであるような確率P(E)は、

P(E)=\displaystyle\frac{\Omega(E)e^{-\beta E}}{Z(\beta)}

で与えられる。ここにZ(\beta)は分配関数と呼ばれ、規格化定数のようなものである。そしてそれは、

Z(\beta)=\int dE\, \Omega(E)e^{-\beta E}

で計算される。
この他にもグランドカノニカルアンサンブルなども同様に化学ポテンシャル\muと温度Tの熱浴を考え、それと対象系を含めたものにミクロカノニカルアンサンブルを適用すれば得ることができる。
各アンサンブルについてこのように物理的な状況を考え、それに対してミクロカノニカルアンサンブルを適用し新たなアンサンブルを得てもよいが、もっと形式的に新たなアンサンブルを得る方法がないか考えてみよう。

分配関数について

カノニカルアンサンブルにおいて規格化定数のような役割を果たす分配関数についてもう少し考えてみる。分配関数とは、

Z(\beta)=\int dE\, \Omega(E)e^{-\beta E}

であった。この積分を計算することを考える。まず、Boltzmann原理を用いて状態密度\Omega (E)をエントロピーS(E)を用いて書き換えると、

Z(\beta)=\displaystyle\int dE\, \exp \left[\frac{1}{k_{\mathrm{B}}}\left(S(E)-\frac{E}{T}\right)\right]

となる。そして、被積分関数が指数関数であることから鞍点法を用いて積分を評価しよう。

\displaystyle\frac{d}{dE}\left(\frac{1}{k_{\mathrm{B}}}\left(S(E)-\frac{E}{T}\right)\right)=0

であるような点をE^*と書くことにすると、被積分関数はこのE^*の近くでのみ大きな値を持つのでこの付近でのみ積分を考えればよい。すると2次までの近似で被積分関数は、

\displaystyle\exp\left[\frac{1}{k_{\mathrm{B}}}\left(S(E^*)-\frac{E^*}{T}+\frac{1}{2}\frac{d^2 S(E^*)}{dE^2}(E-E^*)^2\right)\right]

と近似される。積分されるのは2次の項であり、それはGauss積分で計算でき、結局

Z(\beta)=\sqrt{\frac{k_{\mathrm{B}}\pi}{2\left|\frac{d^2 S(E^*)}{dE^2}\right|}}\exp\left[\frac{1}{k_{\mathrm{B}}}\left(S(E^*)-\frac{E^*}{T}\right)\right]

となる。
さて、この分配関数の対数をとれば係数の部分はO(\sqrt{N})なので無視できることを考えると、

\displaystyle\log Z(\beta)=\frac{1}{k_{\mathrm{B}}}\left(S(E^*)-\frac{E^*}{T}\right)

となり、Mathieu関数\Phi(T)あるいは、Helmholtz自由エネルギーF(T)が得られる。

\Phi(T)=k_{\mathrm{B}}\log Z(\beta)
F(T)=\displaystyle-\frac{1}{\beta}\log Z(\beta)

これらは、ミクロカノニカルアンサンブルにおけるBoltzman原理のカノニカルアンサンブルにおける対応物と考えられる。
つまり、ミクロカノニカルアンサンブルにおいて規格化定数のような役割を果たしていた状態密度\Omega(E)がBoltzman原理によってマクロ系のエントロピーと結びついていたのと同様に、カノニカルアンサンブルでは規格化定数にあたる分配関数Z(\beta)がマクロ系のMathieu関数(あるいはHelmholtz自由エネルギー)と結びついていることがわかる。
では、Mathieu関数（あるいはHelmholtz自由エネルギー）とは熱力学においてどのような量であっただろうか。
エントロピーの微分を思い出すと、

\displaystyle dS=\frac{1}{T}dU+\frac{p}{T}dV

なのであった。そして、これをUについてLugendre変換したものがMathieu関数（あるいはHelmholtz自由エネルギー）である。

\Phi\left[\frac{1}{T};V\right]=\displaystyle S-\frac{U}{T}
F[T;V]=-T\Phi

ここで、分配関数の定義を見直せば、これは状態密度\Omega(E)の\betaを変数とするLaplace変換になっている。即ち、熱力学においてLugendre変換によって各熱力学関数が結びついていたのと対応して、統計力学では各アンサンブルがLpalace変換によって結びつくことになる。

形式的にミクロカノニカルアンサンブルからカノニカルアンサンブルを得る方法

先の過程を逆にたどってみよう。
まず第一にエントロピーからLugendre変換で結びつく何らかの熱力学関数(今の場合はMathieu関数)がミクロな系の量とBoltzman原理のように結びついてほしい。つまり、

\displaystyle\Phi(T)=k_{\mathrm{B}}\log Z(T)

という関係にあってほしい。ここでZ(T)はまだわからないミクロカノニカルアンサンブルにおける\Omega(E)にあたる量であり、\Phi(T)はMathieu関数である。
次にこのようなZ(T)を構成することを考える。上の式を変形すると、

Z(T)=\displaystyle\exp\left[\frac{\Phi}{k_\mathrm{B}}\right]=\exp\left[\frac{1}{k_{\mathrm{B}}}\left(S-\frac{E}{T}\right)\right]

となる。
さて、Mathieu関数においてEはLugendre変換から\displaystyle \frac{dS}{dE}=\frac{1}{T}によって定まっている*2ので、上の式におけるEもこの関係を満たしていなければならない。すると、

Z(T)=\displaystyle\int dE\,\exp\left[\frac{\Phi(E)}{k_{\mathrm{B}}}\right]=\int dE\,\exp\left[\frac{1}{k_{\mathrm{B}}}\left(S(E)-\frac{E}{T}\right)\right]

というZ(T)ならば鞍点法を用いればLugendre変換から決まるEに関する条件を満たしてなおかつ対数を取った時に熱力学関数が出るような規格化定数になる。そして、最後にBoltzman原理を用いれば

Z(T)=\displaystyle\int dE\,\Omega(E)e^{-\beta E}

となり、これによってエネルギーがEの時の相対確率が\Omega(E)e^{-\beta E}と決まる。
つまり、分配関数Z(T)を構成したければ、エントロピーからLugendre変換したい変数EについてLugendre変換後の自然な変数\displaystyle\frac{1}{T}をBoltzman定数k_{\mathrm{B}}で割ったもの（つまり逆温度\beta）を変数とするLaplace変換すればよいことがわかる：

\Omega(E)\overset{\mathrm{Laplace}}{\longrightarrow} Z(T)=\displaystyle\int dE\,\Omega(E)e^{-\beta E}

こうしてMathieu関数とBoltzman原理からカノニカルアンサンブルを構成できた。
以上の手続きを形式的にまとめると次のようにまとめられる。

  1. エントロピーS(E)から使いたい、即ちBoltzman原理のようにそれによってミクロ系とマクロ系の対応を与える熱力学関数をLugendre変換によって求める。
  2. 状態密度\Omega(E)をLugendre変換によって変えた後の自然な変数（をBoltzman定数k_{\mathrm{B}}で割ったもの）を用いてLaplace変換し、それを新たなアンサンブルにおける分配関数として定義する。
  3. 等重率の仮定より分配関数からある状態の確率を決める。

グランドカノニカルアンサンブル

形式的にアンサンブルを構成する方法の応用としてミクロカノニカルアンサンブルからグランドカノニカルアンサンブルを求めてみる。
エントロピーS[U,V,N]について、その微分は、

\displaystyle dS=\frac{1}{T}dU+\frac{p}{T}dV-\frac{\mu}{T}dN

であった。
まず、エントロピーの変数のUとNについてLugendre変換し、新たな熱力学関数\Phi_{\mu}\left[\frac{1}{T},-\frac{\mu}{T};V\right]を作る。

\Phi_{\mu}\left[\frac{1}{T},-\frac{\mu}{T};V\right]=\displaystyle S-\frac{U}{T}+\frac{\mu N}{T}

ここで、UとNについては、Lugendre変換から\displaystyle \frac{dS}{dU}=\frac{1}{T}と\displaystyle \frac{dS}{dN}=-\frac{\mu}{T}によって決まっている。
この熱力学関数 \Phi_{\mu}\left[\frac{1}{T},-\frac{\mu}{T};V\right]のLugendre変換後の自然な変数について、Uが\displaystyle\frac{1}{T}に、Nが\displaystyle-\frac{\mu}{T}になっているので、変換後の自然な変数をBoltzman定数 k_{\mathrm{B}}で割ったもの、つまり\betaと-\beta\muを変数とするLaplace変換を行えばよい。つまり、

\Omega(E)\overset{\mathrm{Laplace}}{\longrightarrow}\Xi\left(\beta,\mu\right)=\int dE\int dN\,\Omega(E)e^{-\beta E+\beta\mu N}

であり、この\Xi\left(\beta,\mu\right)がグランドカノニカルアンサンブルにおける大分配関数である。
そして、\Xi\left(\beta,\mu\right)を大分配関数としたことから、エネルギーがEであり、粒子数がNであるような状態の出現確率P(E,N)が、

P(E,N)=\displaystyle \frac{e^{-\beta E+\beta\mu N}}{\Xi\left(\beta,\mu\right)}

と決まる。
こうしてBoltzman原理とミクロカノニカルアンサンブルからグランドカノニカルアンサンブルが構成することができた。
なお、構成方法からわかるようにグランドカノニカルアンサンブルにおけるBoltzman原理の対応物は、

\Phi_{\mu}\left[\frac{1}{T},-\frac{\mu}{T};V\right]=k_{\mathrm{B}}\log\Xi\left(\beta,\mu\right)
J\left[T,\mu;V\right]=\displaystyle-\frac{1}{\beta}\log\Xi\left(\beta,\mu\right)

である。ここで、J\left[T,\mu;V\right]はグランドポテンシャル。

T-pアンサンブル

T\mbox{-}pアンサンブルについても同様に導出してみよう。
まず、エントロピーS[U,V,N]について、UとVについてLugendre変換し、新たな熱力学関数*3\Psi\left[\frac{1}{T},\frac{p}{T}\right]を構成する：

\Psi\left[\frac{1}{T},\frac{p}{T}\right]=\displaystyle S-\frac{E}{T}-\frac{pV}{T}

ここで、UとVについては、Lugendre変換から\displaystyle \frac{dS}{dU}=\frac{1}{T}と\displaystyle \frac{dS}{dV}=\frac{p}{T}によって決まっている。
この熱力学関数\Psi\left[\frac{1}{T},\frac{p}{T}\right]のLugendre変換後の自然な変数について、Uが\displaystyle\frac{1}{T}に、Vが\displaystyle\frac{p}{T}になっているので、変換後の自然な変数をBoltzman定数k_{\mathrm{B}}で割ったもの、つまり\betaと\beta pを変数とするLaplace変換を行えばよい。つまり、

\Omega(E)\overset{\mathrm{Laplace}}{\longrightarrow}Y\left(\beta,p\right)=\int dE\int dV\,\Omega(E)e^{-\beta E-\beta pV}

であり、このY\left(\beta,p\right)がT\mbox{-}pアンサンブルにおける分配関数である。
そして、Y\left(\beta,p\right)を分配関数としたことから、エネルギーがEであり、体積がVであるような状態の出現確率P(E,V)が、

P(E,V)=\displaystyle \frac{e^{-\beta E-\beta pV}}{Y\left(\beta,p\right)}

と決まる。
こうしてBoltzman原理とミクロカノニカルアンサンブルからT\mbox{-}pアンサンブルを構成することができた。
なお、構成方法からわかるようにT\mbox{-}pアンサンブルにおけるBoltzman原理の対応物は、

\Psi\left[\frac{1}{T},\frac{p}{T}\right]=k_{\mathrm{B}}\log Y\left(\beta,p\right)

である。

その他のアンサンブル

同様のことを行えば物理的に使えるかどうかはともかくとして様々なアンサンブルを作ることが可能である。
また、ゴムの熱力学などで使えるT\mbox{-}fアンサンブルについてもこの手法を用いて、それに対応するアンサンブルを構成することができる。この場合エントロピーの微分は、fを張力、lを変位として

dS=\displaystyle\frac{1}{T}dU-\frac{f}{T}dl

なので、UとlについてLugendre変換した熱力学関数\Phi_{f}\left[\frac{1}{T},-\frac{f}{T}\right]を考え、それに対応する分配関数Y(\beta,f)を構成してやれば、

Y(\beta,f)=\int dE\int dl\,\Omega(E)e^{-\beta E+\beta fl}

となる。そして、エネルギーがE、変位がlの状態の出現確率P(E,l)は、

P(E,l)=\displaystyle\frac{e^{-\beta E+\beta fl}}{Y(\beta,f)}

となる。

*1:例えば、キャレン「熱力学および統計物理入門（下）」第16章などを参照。

*2:これは熱浴によって対象系の温度がTに保たれていることを意味している。

*3:Planck関数、あるいはMathieu関数と呼ぶらしい。