量子カオスの定義を検索してもどれも曖昧なものばかりで、定義の明言を避けている印象があります。この記事の目的は、特に Bohigas, Giannoni, and Schmit (1984) による予想について、具体的にシナイ・ビリヤードと呼ばれる物理系の実装を添えて量子カオスとはなんなのか、その輪郭をざっくり掴むことにあります。数式を限りなく減らしたので、あまり身構えず気軽に読んでください。

(注：この記事は、物工/計数 Advent Calendar 2022のために書かれたものです。私は2022年度現在学部3年であり、この分野の専門家でもなんでも無いので、間違い・思い込みなどたくさんあるかと思います。見つけていただいた方はぜひご連絡いただければ幸いに思います。)

古典力学系におけるカオス

量子カオスの前に、古典カオスについて紹介しましょう。

カオスの定義

カオスという言葉は日常会話でも比較的使われることが多く、多くの場合「混沌」の意味で使われているかと思います。しかし、力学系ではそんなに曖昧な定義ではありません。今回は Strogatz (2018) に準じた定義を採用します。以下、Strogatz (2018) からの引用です：

今のところカオスの定義として万能なものは知られていないが、ほとんどの人は、以下の暫定的な定義に含まれている3つの構成要素に同意するだろう。
「 カオスは決定論的なシステムにおける非周期的な長時間挙動であり、初期条件への鋭敏な依存性を示す。」
(1) 「非周期的な長時間挙動」とは、$t\to \infty$ で固定点、周期軌道、あるいは準周期軌道のいずれにも落ち着かない軌道が存在することを意味している。現実的な要請として、そのような軌道が、あまりにまれなものではないことが必要となる。例えば、この比周期軌道へ至る初期条件が開集合として存在することや、ランダムな初期条件に対してゼロではない確率でそのような軌道が生じることを要求する必要があるだろう。
(2) 「決定論的」とは、系がランダムな、あるいはノイズ的な入力やパラメータを持たないことを意味している。不規則な挙動は、ノイズ的な外力からではなく、系の非線形性から生じる。
(3) 「初期条件への鋭敏な依存性」とは、近接軌道が指数関数的に速く分離していくこと、つまり正のリアプノフ指数を持つことを意味する。

なんだか仰々しい定義ですね。ここまで定義が長いのに万能な定義ではないとはどういうことだと思われるかもしれませんが、例えば $\dot{x}=x$ のような、一見カオスに見えるがそうではない一部の例外を強く意識しているからだと思われます。

この記事ではあまり細かいことは考えません。上の定義の気持ちだけ書けば、以下のようになると思います（個人の感想です）。

カオスとは、性質が同じ2つの物体を力学系の微妙に違う初期条件からスタートさせたとき、そのうち軌道が全く異なるものになる様子。

イメージとしてはこんな感じです。（図はリアプノフ指数のWikipediaより引用。Yapparina, CC0, via Wikimedia Commons）

コンピュータが扱える精度の限界より遥かに小さい初期値のズレからでもカオスは起こるので、どれだけ正確に初期条件を決めても長時間の挙動予測が不可能になります。たとえ「秩序のある」初期条件からダイナミクスが始まっても、長時間経った後では秩序が失われてしまうことになります。

日常会話で「混沌」を意味する語として使われているカオスですが、力学では上述のような定義や性質を持っていて、まさに「混沌」そのもの、あるいはその原因といえそうですね。

カオスを示す力学系の例

有名なものには、 二重振り子 や 天体の運動(三体問題) などが挙げられます。カオス的な挙動を示す系の多くは 非可積分系 と呼ばれる、自由度に対して保存量が足りない系 となっています。多体問題と呼ばれる、3以上の自由度が相互作用を起こす系では多くが非可積分であることが知られています。

一方で、上述のような複雑な系でなくてもカオス的な挙動を起こすことが知られています。この記事で主に注目するシナイ・ビリヤードはそのような例の1つです。

シナイ・ビリヤード

二次元の箱の中心に円盤状の障壁を置き、その中で質点を運動させる事を考えます。ハミルトン力学系として捉えるためには、円盤の内側と箱の外側ではポテンシャルが無限に大きく、それ以外の場所ではポテンシャルが0であるような系を考えればよいです。

質点はどのような運動を行うでしょうか。Processingを使える環境をお持ちの方は、ぜひ以下のサンプルコードをまるまるコピペして実行してみてください。

sinai.pde

int Nparticle = 1000;
float LX=400;
float LY=400;
float R=min(LX/2,LY/2);
float v0=700;
float theta=PI/9;

float[] x =new float[Nparticle];
float[] y =new float[Nparticle];
float[] cx=new float[Nparticle];
float[] cy=new float[Nparticle];
float[] vx=new float[Nparticle];
float[] vy=new float[Nparticle];

float m=1;
float dt=0.01;

int count=0;
void setup() {
  size(900,900);
  background(0);
  for (int i=0;i<Nparticle;i++){
    x[i]=-R+random(1);
    y[i]=-R+random(1);
    vx[i]=v0*cos(theta);
    vy[i]=v0*sin(theta);
  }
}

void draw() {
  translate(width/2,height/2);
  count+=1;
  background(0);
  stroke(255);
  strokeWeight(1);
  line(LX,LY,LX,-LY);
  line(LX,-LY,-LX,-LY);
  line(-LX,-LY,-LX,LY);
  line(-LX,LY,LX,LY);
  ellipse(0,0,R*2,R*2);
  fill(255);
  for (int i=0; i<Nparticle; i++){
    float[] hold=new float[2];
    x[i] += vx[i]*dt;
    y[i] += vy[i]*dt;
    if (x[i]<-LX){
      x[i]=-2*LX-x[i];
      vx[i]=-vx[i];
    }else if(x[i]>LX){
      x[i]=2*LX-x[i];
      vx[i]=-vx[i];
    }
    if (y[i]<-LY){
      y[i]=-2*LY-y[i];
      vy[i]=-vy[i];
    }else if(y[i]>LY){
      y[i]=2*LY-y[i];
      vy[i]=-vy[i];
    }
    if (x[i]*x[i]+y[i]*y[i]<R*R){
      float T=optimaltime(x[i],y[i],vx[i],vy[i],R);
      cx[i]=x[i]-vx[i]*T;
      cy[i]=y[i]-vy[i]*T;
      hold=mirror(x[i],y[i],cx[i],cy[i],vx[i],vy[i],T);
      x[i]=hold[0];
      y[i]=hold[1];
      vx[i]=(x[i]-cx[i])/sqrt((x[i]-cx[i])*(x[i]-cx[i])+(y[i]-cy[i])*(y[i]-cy[i]))*v0;
      vy[i]=(y[i]-cy[i])/sqrt((x[i]-cx[i])*(x[i]-cx[i])+(y[i]-cy[i])*(y[i]-cy[i]))*v0;
    }
    fill(255);
    ellipse(x[i],y[i],1,1);
  }
}

float optimaltime(float x,float y,float vx,float vy,float R){
  float V2=vx*vx+vy+vy;
  float RV=x*vx+y*vy;
  float T=(RV+sqrt(RV*RV+V2*(R*R-x*x-y*y)))/V2;
  return T;
}

float[] mirror(float x,float y,float cx,float cy,float vx, float vy, float t){
  float theta=findtheta(cx,cy);
  float phi=findtheta(x-cx,y-cy);
  float[] ret= new float[2];
  ret[0]=cx+(-cos(2*(theta-phi))*vx+sin(2*(theta-phi))*vy)*t;
  ret[1]=cy+(-sin(2*(theta-phi))*vx-cos(2*(theta-phi))*vy)*t;
  return ret;
}


float findtheta(float x,float y){
  float theta=0;
  if (y>0){
    theta=acos(x/sqrt(x*x+y*y));
  } else if (y<0){
    theta=PI+acos(-x/sqrt(x*x+y*y));
  }
  return theta;
}

このシミュレーションでは、相互作用しない1000個の質点を用いて、全く同じ初期速度で、ほんの少しだけ異なる初期位置から計算を開始しています。詳しい実装法などを解説すると長くなりすぎるので省きますが、基本的には障壁に衝突した点に引かれる接線に反射の法則を適用しただけです。

Processingは環境設定がほとんど不要で強力な機能が豊富なC++ベースのプログラミング環境であり、ジェネラティブアートの作成などに広く使われています。簡単なシミュレーションを書くのにとても重宝しますので、皆さんもぜひ使ってみてくださいね。

さて、もしエラーなくコードが動けば、以下のような挙動が見られたはずです：

確かに、上述したカオスの定義のとても良い例になっています。この系はシナイ・ビリヤードと呼ばれる系で、次に述べる量子カオスの特徴づけを行う上で重要な役割を果たすことになりました。

量子力学系におけるカオス

ここでは、ランダム行列が量子カオスと関わりを持つに至る経緯を説明し、量子力学的なシナイ・ビリヤードのエネルギー準位をPythonを用いて数値計算します。筆者も勉強中の部分が多く、ほとんどが文献頼りになりますがご了承ください。また、記述の少なくない部分で Ullmo (2016) を参考にさせていただきました。

Wigner (1951) による原子物理へのランダム行列の導入

水素のRydbergエネルギースペクトル。摂動が加わると準位が交差することがある。Michael Courtney氏の著作、CC BY-SA 3.0, via Wikimedia Commons

反磁性リチウムのRydbergエネルギースペクトル。摂動が加わっても準位が交差しないことがわかる。Michael Courtney氏の著作、CC BY-SA 3.0, via Wikimedia Commons

高励起状態において電場を加えていくと、水素のエネルギースペクトルは交差することがあるのに対して、反磁性リチウムのエネルギースペクトルは交差しない準位反発という現象を示します。古典的に可積分な系の多くは水素のエネルギースペクトルに類似した現象を示すのに対し、古典的に非可積分な系はリチウムのような準位反発を示すことが知られています。

Wigner (1951) は、このような不思議な現象には系の複雑性やランダムネスが背後にあると考え、十分に複雑な系では、系の対称性以外の情報は全く寄与しないと仮定し、適切に設定したランダムなハミルトニアンで近似できると仮定しました。このようなランダムなハミルトニアンは、行列として捉えると、ランダム行列と呼ばれる分野で研究されています。

結果として、ランダム行列のスペクトルで準位反発を記述することに成功し、Wignerが提唱したランダム行列による準位反発の記述は理論と実験で精度良く一致しました。

準位統計

準位統計は、スペクトルの統計的性質を明らかにするための手段です。次節での展開をよりよく理解するために、準位反発が起きているかそうでないかを記述するための準位統計の一つの方法である、Nearest Neighbor Spacingについて説明します。

Nearest Neighbor Spacingとは、最も近いほかのスペクトルとの距離のことです。上図でn番目のエネルギー固有値に対応するNearest Neighbor Spacingは、$E_n-E_{n-1}$ということになります。注目する固有値が縮退している場合、この値は0になります。準位反発が起きる場合、この値は0ではなくなりそうです。

ハミルトニアンの全てのエネルギー準位を計算し、横軸にNearest Neighbor Spacing、縦軸にその値の出現回数をプロットした分布は、準位反発が起きる系とそうでない系で、0近傍で異なる振る舞いを示すと予想されます。

具体的には、準位反発なければ縮退が容易に起きうるので、Nearest Neighbor Spacingが0となる場合が多くなると考えられます。その場合は0で極大を取る分布になるはずですが、準位反発の効果が大きくなるにつれ、極大が移動すると考えられます。下図はその様子を表したもので、(a)から(h)にかけて分布が変化していく様子を表しています(Michael Courtney氏の著作、CC BY-SA 3.0, via Wikimedia Commons)。

おそらくランダム行列の文脈から導入されたものだと思います（筆者はランダム行列をちゃんと勉強していないので不安）が、あるスペクトルに注目して、その上下のスペクトルとの間隔を計算し、大きいスペクトル間隔で小さいスペクトル間隔の値を割ったもの（下図では$r_n$と表記）を指標とする準位統計も有力な方法です。

Atas et al. (2013) にはその期待値とランダム行列の分類の関係が細かく論じられていて、現在も研究でよく使われています。

Bohigas, Giannoni, and Schmit (1984) の予想

Wignerによって導入されたランダム行列による準位反発の記述は、その後Wignerの意図に反して（？）系の複雑性というよりも系の量子カオスの結果だと認識されることになります。

まず、Berry and Tabor (1977) により、可積分な系では準位反発が生じないことが示されました。一方で、Bohigas, Giannoni, and Schmit (1984) により、以下の予想が打ち立てられました(論文の文をそのまま日本語訳したもの)：

古典的な類似系がK-Systemである時間反転不変系のスペクトルは、Gaussian Orthogonal Ensembleが予言したのと同じ揺らぎ特性を示す。

これがBohigas, Giannoni, and Schmit (1984) の予想と呼ばれるものです。難解な用語が多くて意味を汲み取りづらいですが（筆者もエルゴード理論はあまり勉強したことがありません）、大体の意味は以下の通りになると思います（間違っていたらすいません）。

古典的なカオスを示す系の量子的なスペクトルの性質は、Wignerらによって導入されたランダム行列の一種のスペクトルの性質と似ている。

ということになると思います。Wignerの主張では系の複雑性がランダム行列の適用の根拠でしたが、Bohigasらはカオスが本質だと指摘しているのだと思います。

この根拠として、Bohigasらは、対称性を落としたシナイ・ビリヤード（や他にも様々なダイナミカルビリヤード）のエネルギー準位からNearest Neigbor Spacingの分布を計算しています。

量子シナイ・ビリヤード

せっかくですので、私達もBohigasらの計算結果が正しいのか確認してみましょう。以下、Pythonを用いた実装です。IPython Notebookを用いることを前提にしています。一つ一つのコードブロックは.ipynb形式のコードブロックと対応させているので、コピペを行う際は順番を守るように気をつけてください。

Bohigasらはビリヤードの対称性を落としていましたが、簡単のためにProcessingで実装した古典系と同じポテンシャルから始めてみましょう。まずはビリヤードをポテンシャルを関数として定義します。$V_0$ はハードコアポテンシャルを近似するために正の巨大な数を代入することにします。

sinai.ipynb

def potential(x,y,LX,LY,R,V0):
    if x<-LX or x>LX:
        return V0
    elif y<-LY or y>LY:
        return V0
    elif x**2+y**2<R**2:
        return V0
    else:
        return 0

Nearest Neighbor Spacingを計算するための関数を定義しておきます（もしかすると筆者の不勉強のせいでトリビアルな部分の定義が間違っているかもしれませんが、結果に大きく影響することはないでしょう。あまり考えずコピペして動かすのが良いと思います）。

import numpy as np
import matplotlib.pyplot as plt

def NNSD(eigenvalues):
    N=len(eigenvalues)
    eigenvalues=np.array(sorted(np.real(eigenvalues)))
    NNSD_list=np.zeros(N)
    for i in range(1,N-2):
        NNSD_list[i]=min(eigenvalues[i+1]-eigenvalues[i], eigenvalues[i]-eigenvalues[i-1])
    NNSD_list[0]=eigenvalues[1]-eigenvalues[0]
    NNSD_list[N-1]=eigenvalues[N-1]-eigenvalues[N-2]
    Delta=(eigenvalues[-1]-eigenvalues[0])/N
    NNSD_list=NNSD_list/Delta
    return NNSD_list

def ERatio(eigenvalues):
    N=len(eigenvalues)
    eigenvalues=np.array(sorted(np.real(eigenvalues)))
    ERatio_list=np.zeros(N-1)
    for i in range(1,N-1):
        di=eigenvalues[i]-eigenvalues[i-1]
        di1=eigenvalues[i+1]-eigenvalues[i]
        ERatio_list[i]=(min(di,di1))/(max(di,di1))
    return ERatio_list

定義した関数が本当にシナイ・ビリヤードになっているのか確認しましょう。まずはシステムサイズの設定をします。

LX=100
LY=100
V0=100000
R=min(LX/2,LY/2)

ポテンシャルをヒートマップにプロットしてみます。シナイ・ビリヤードと同じ形状であることが確認できるはずです。

pot_test=np.zeros((4*LX+1,4*LY+1))
for x in range(4*LX+1):
    for y in range(4*LY+1):
        pot_test[x,y]=potential(x-2*LX,y-2*LY,LX,LY,R,V0)
plot1=plt.subplot()
plot1.pcolor(pot_test)
plot1.set_xlim(0,LX*4)
plot1.set_ylim(0,LY*4)
plot1.set_title("potential heatmap")
plot1.set_aspect('equal')

では量子系の対応するハミルトニアンを作っていきます。一応システムサイズを直感的な値に再設定し、系の空間方向の分割数を設定します。divx, divyに対応した数だけ線形独立な状態 $|x\rangle,\,|y\rangle$ が存在するということです。

LX=1
LY=1
R=min(LX/2,LY/2)
m=1
divx=40
divy=40

ハミルトニアンの行列表示を与える関数を定義します。ラプラシアンを数値微分に置き換えているだけですのであまり気にせずに実行してみてください。

def Hamiltonian(LX,LY,divx,divy,m,V0,sinai):
    dx=LX/divx
    dy=LY/divy
    A=(1/dx**2+1/dy**2)/(2*m)
    mAh=-A/2
    Hmat=np.zeros(((2*divx+1)*(2*divy+1),(2*divx+1)*(2*divy+1)))
    for indx in range(2*divx+1):
        x=(indx-divx)/(2*divx+1)
        for indy in range(2*divy+1):
            y=(indy-divy)/(2*divy+1)
            if sinai==True:
                Hmat[indx*(2*divy+1)+indy,indx*(2*divy+1)+indy]=A+potential(x,y,LX,LY,R,V0)
            else:
                Hmat[indx*(2*divy+1)+indy,indx*(2*divy+1)+indy]=A
            if indx!=2*divx:
                Hmat[(indx*(2*divy+1)+indy)  , (indx+1)*(2*divy+1)+indy]=mAh
            if indx!=0:
                Hmat[(indx*(2*divy+1)+indy)  , (indx-1)*(2*divy+1)+indy]=mAh
            if indy!=0:
                Hmat[(indx*(2*divy+1)+indy)  , indx*(2*divy+1)+indy-1  ]=mAh
            if indy!=2*divy:
                Hmat[(indx*(2*divy+1)+indy)  , indx*(2*divy+1)+indy+1  ]=mAh
    return Hmat

実際に行列表示を見てみましょう。大規模な疎行列が得られるはずです。

Hmat=Hamiltonian(LX,LY,divx,divy,m,V0,True)
print(Hmat)

得たのは疎行列なので、疎行列用のライブラリを用いて固有値を計算したいところですが、全てのエネルギー準位を計算しなければ準位統計を計算することはできないので、全ての固有値を求めます。マシンスペックや系の空間分割数にもよりますが、筆者のラップトップではdivx=divy=40程度で計算に30秒程度かかりました。もしあまりにも計算に時間がかかるようでしたら、 divxやdivyの値を減らして前のセルから実行し直してみてください。

eps,psi=np.linalg.eigh(Hmat)

固有値を小さいものから見てみましょう。

for i in range(28):
    print(eps[i])

実行がうまく行けば、ほとんど4重縮退しているように見えるはずです。これは系が離散的な回転対称性を持っているからだと考えられます。というのは、ビリヤードの中心を原点として全体を $0,\frac{\pi}{2},\pi,\frac{3\pi}{2}$ ラジアンだけ回転させるともとのビリヤードに一致します。このような対称性を4回回転軸をもつと表現したりします。

ハミルトニアンに対称性がある場合、縮退が起きうることはよく知られています。 Bohigasらが計算したシナイ・ビリヤードは対称性を落としたものだったので、おそらく縮退は起きなかったでしょう。

Bohigasらのビリヤードのスペクトルと似たものを得るために、ここでは4重縮退した固有値の中から1つだけを取ってきて、その固有値だけで準位統計を計算することにします。

Nndeg=int(np.floor((2*divx+1)*(2*divy+1)/4))
epsq=np.zeros(Nndeg)
for i in range(Nndeg):
    epsq[i]=eps[4*i]

では準位統計を計算してみましょう。

plt.clf()
plt.hist(NNSD(epsq),bins=50)
plt.xlabel("Nearest neighbor energy spacing")
plt.ylabel("Number of occurrences")
plt.title(f"Nearest Neighbor Energy Spacing - Distributuion")
plt.show()

うまくいけば、分布の極大は0ではない場所にあるはずです。このことは、準位反発が起きてエネルギースペクトルのNearest Neighbor Spacingが平均的に0よりも大きくなることを意味します。

可積分な系では、このように極大が0からずれることは起こらないはずです。ハミルトニアン構築の際に、最後のパラメータをFalseにするとビリヤード内部のポテンシャルを0にする（つまり2次元の井戸型ポテンシャルを作る）ことができます。すると可積分なので準位反発は起きず、準位統計は0付近で巨大な値を示すことがわかりますので、是非試してみてください。

また、スペクトル間隔の比率の分布と期待値を計算してみましょう。 Atas et al. (2013) によれば、期待値が $4-2\sqrt{3}$ となれば、ランダム行列の一種になるようです。

Nbins=50
X=plt.hist(ERatio(epsq),bins=Nbins)
plt.xlabel("r")
plt.ylabel("Number of occurrences")
plt.title(f"Nearest Neighbor Energy Spacing - Distributuion")
plt.show()
print("Theory:",4-2*np.sqrt(3))
print("Sample:",np.dot(X[0],X[1][0:Nbins])/len(epsq))

非常に近い値になっていて、ランダム行列の一種とスペクトルの性質がよく似ているとわかるはずです。 Bohigasらの主張がなんとなく理解できたかと思います。

参考文献

Atas, Y. Y., E. Bogomolny, O. Giraud, and G. Roux. 2013. “Distribution of the Ratio of Consecutive Level Spacings in Random Matrix Ensembles.” Phys. Rev. Lett. 110 (February): 084101. https://doi.org/10.1103/PhysRevLett.110.084101.

Berry, Michael Victor, and Michael Tabor. 1977. “Level Clustering in the Regular Spectrum.” Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 356 (1686): 375–94. https://royalsocietypublishing.org/doi/10.1098/rspa.1977.0140

Bohigas, O., M. J. Giannoni, and C. Schmit. 1984. “Characterization of Chaotic Quantum Spectra and Universality of Level Fluctuation Laws.” Phys. Rev. Lett. 52 (January): 1–4. https://doi.org/10.1103/PhysRevLett.52.1.

Martı́nez-Peña, Rodrigo, Gian Luca Giorgi, Johannes Nokkala, Miguel C. Soriano, and Roberta Zambrini. 2021. “Dynamical Phase Transitions in Quantum Reservoir Computing.” Phys. Rev. Lett. 127 (August): 100502. https://doi.org/10.1103/PhysRevLett.127.100502.

Strogatz, Steven H. 2018. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. CRC press.

Ullmo, Denis. 2016. “Bohigas-Giannoni-Schmit Conjecture.” Scholarpedia 11 (9): 31721. http://www.scholarpedia.org/article/Bohigas-Giannoni-Schmit_conjecture.

Wigner, Eugene P. 1951. “On the Statistical Distribution of the Widths and Spacings of Nuclear Resonance Levels.” Mathematical Proceedings of the Cambridge Philosophical Society 47 (4): 790–98. https://doi.org/10.1017/S0305004100027237.

量子カオスって何だろう？