動径方向に対する密度構造とポテンシャルの振動成分は、ルジャンドル陪多項式と呼ばれる関数が関わっており、複数のモード（m、l、n）が存在します。今回はm=2、l=2で計算を行っています。下に動径方向の密度の振動成分とポテンシャルの振動成分のグラフを掲載します。
図１ 横軸は距離、縦軸が密度の振動成分
図２ 横軸は距離、縦軸がポテンシャルの振動成分
n=1～6に対して、blue、green、purple、navy、orange、redが対応しています。この密度振動の分布がポテンシャルの振動によってどう変化するかという方程式が存在するのですが、様々なモードに対して方程式が存在するため、matrix（行列）として表記することができます。matrixは回転の効果を含まないmatrix（図３左）と回転の効果を含むmatrix（図３右）が存在します。これはKalnajsさんという方が考えた理論です。この2つのmatrixを足し合わせたmatrixの固有ベクトルを（動径方向の密度振動×球面調和関数）にかけて、nについて足し合わせて、実数部分（複素数のreal part）をプロットしたものが最初のページに掲載した3次元プロット図になっています。
図３ Kalnajs matrix
なお回転の効果を含むmatrixに係数αがついており、このαの値と速度分散の非等方性を導入する変数qの値によって密度振動のspiral patternがどの程度強く表れるかが決まります。最初のページで示したspiral patternはα=0.9、q=2としています。次のページでαの値によってどのようなspiral patternが生じるかが確認できます。
今回はpythonを使ってプログラミングを行い、matplotlibと呼ばれるモジュールを利用してグラフを作成しました。pythonでは、scipyと呼ばれる超幾何関数、ルジャンドル陪関数を扱えるモジュールが存在し、numpyと呼ばれるモジュールでは行列の固有値と固有ベクトルを簡単に計算することができます。複素数の取り扱いも簡単にできます。今回のシミュレーションを行う上では、複素関数、解析力学における、作用・角変数の知識が必要であり、これらの知識を組み合わせてシミュレーションを行っているのは、他では例があまりないです。
