調和代数幾何解析(調和解析と不変式論)

GoodmanとWallachの本

Symmetry, Representations, and Invariants
https://www.amazon.co.jp/dp/1441927298

の12章を眺めていて、ふと、ここに書いてあるのは、調和解析+代数幾何 = 調和代数幾何解析だなぁと思った(勝手に命名)。まぁ、この本の範囲では、"解析学"は、全然出てこないのだけど


【参考文献】同じく、GoodmanとWallachの本

Representations and Invariants of the Classical Groups
https://www.amazon.co.jp/Representations-Invariants-Encyclopedia-Mathematics-Applications/dp/0521663482

の12章あたりにも同等の内容があるっぽい。以下のPDFにも、同様の内容がある

An Algebraic Group Approach to Compact Symmetric Spaces
https://sites.math.rutgers.edu/~goodman/pub/symspace

Harmonic Analysis on Compact Symmetric Spaces : the Legacy of Élie Cartan and Hermann Weyl
https://sites.math.rutgers.edu/~goodman/pub/weyl_goodman.pdf



コンパクト対称空間G/K上の二乗可積分関数の空間L^2(G/K)をG作用で既約分解するという問題は、Cartanが1929年に解いたらしい。この問題が、GoodmanとWallachの本の12章のそもそもの主題である。Cartanの議論には瑕疵があったらしく、1970年にHelgasonがギャップを埋めたということで、Cartan-Helgasonの定理と呼ばれていることがある。古典的な例として、L^2(S^1)のSO(2)作用による分解はFourier級数で、L^2(S^2)のSO(3)作用による分解は球面調和関数展開で記述できる。この定理は、対称空間を舞台とする調和解析の基本と言える。

Cartan-Helgasonの定理には、Peter-Weylの定理と同様に、解析的なバージョンと代数的なバージョンが存在する。

Peter-Weylの定理の場合は、コンパクトLie群Gに対して、L^2(G)の既約分解を記述するのが解析的バージョンで、一方、コンパクトLie群には、実代数群の構造が一意に入る(※)ので、複素化して、群上のアフィン座標環を(Gの複素化群の作用で)既約分解するのが代数的バージョン。

※)代数群側から始めるなら、複素簡約(代数)群G_{\mathbf{C}} \subset GL(n,\mathbf{C})には、G=G_{\mathbf{C}} \cap U(n)として、コンパクト実型が取れると言ってもいい。

個人の好みは別として、Peter-Weylの定理の解析バージョンは、関数空間がHilbert空間だから積構造は捨てられているのに対して、代数バージョンの方は、可換環なので、よりrichな情報を持っている。球面調和関数の積とか議論することを考えると、代数的な方が、面倒事がない。


簡単な例として、SU(2)の場合を考える。SU(2)は、幾何学的にはS^3と同一視出来るので、Peter-Weylの定理の適用対象であると同時に、Cartan-Helgasonの定理の適用対象にもなる。Peter-Weylの定理で考える場合は、SU(2)の作用で分解し、Cartan-Helgasonの定理で考える場合は、例えば、SO(4)などの作用で既約分解する。L^2(SU(2))の基底として、Wigner D-matrixを取ることができ、L^2(S^3)の基底として、hyperspherical harmonicsを取ることができる。

SU(2)を複素化した群は、SL(2,\mathbf{C})で、座標環は\mathbf{C}[a,b,c,d]/(ad-bc-1)と書ける。それで、
\left(\begin{array}{cc} a & b \\ c & d \end{array} \right) = \left( \begin{array}{cc} x_1 + i x_2 & -x_3 + i x_4 \\ x_3+i x_4 & x_1 - i x_2 \end{array} \right)
と変数変換を行えば、
\mathbf{C} [a,b,c,d]/(ad-bc-1)\simeq \mathbf{C}[x_1,x_2,x_3,x_4]/(x_1^2+x_2^2+x_3^2+x_4^2-1)
となって、後者は、3次元球面上の座標環の複素化と思える。

\mathbf{C}[a,b,c,d]/(ad-bc-1)は代数群の座標環なので、Peter-Weylの定理を考えることができる。一次の元は、a,b,c,dと4つあって、これは、二次元既約表現の行列要素を与える。二次の元は、a^2,b^2,c^2,d^2,ab,ac,ad,bc,bd,cdの10個あるけど、関係式ad-bc=1によって、独立な元の数は9個。適当な基底を取って、三次元既約表現の行列要素を見れば、これらの元が全部出ることが分かる。


GoodmanとWallachの本には、複素簡約群とその閉部分群(それぞれG,Kとする)に対して、以下のような環が定義されている
\mathcal{R}(G/K) = \{f \in \mathbf{C}[G] : f(gk) = f(g) \verb|, for all | k \in K,g \in G\}
Goodman-Wallachは、Kの自明表現の誘導表現として、この環を定義している。誘導表現という視点は、GやKが有限群であるような状況で議論を拡張するには有用であるけど、ここでは、\mathcal{R}(G/K)がK作用による不変式環であるという点に着目してみる。

また簡単な例として、\mathcal{R}(SL_{2}/GL_{1})を考える。\mathbf{C}[SL_{2}] = \mathbf{C}[a,b,c,d]/(ad-bc-1)と座標を取って、
\left(\begin{array}{cc} a & b \\ c & d \end{array} \right) \left( \begin{array}{cc} t & 0 \\ 0 & t^{-1} \end{array} \right) = \left( \begin{array}{cc} at & bt^{-1} \\ ct & dt^{-1} \end{array} \right)
なので、ab ,ad, bc, bdは、K=GL_{1}作用で不変である。K=GL_{1}作用で不変な元は、これらの元で生成されることは明らか。4つあるけど、条件ad-bc=1によって、独立な元は3つ以下。

更にad+bc , ab-cd , ab+cdに対して
(ad+bc)^2 + (ab-cd)^2 - (ab+cd)^2 = (ad - bc)^2 = 1
という関係式があって、独立な元は2つとなる。この二次の関係式と球面の関係は明らかだけど、上でやったのと同様に
(a,b,c,d) = (x_1 + i x_2 , -x_3 + i x_4 , x_3+i x_4 , x_1 - i x_2)
と変数変換を行うと
ad+bc = x_1^2+x_2^2-x_3^2-x_4^2
ab-cd = 2(x_1x_3+x_2x_4)
ab+cd = 2i(x_1x_4-x_2x_3)
なので、x_1,x_2,x_3,x_4を実数に制限すると、SU(2) \simeq S^3からS^2へのHopf写像を与えていることが分かる

cf)Hopf fibration
https://vertexoperator.github.io/2018/11/20/elementary_Hopf_fibration.html

以上のような計算によって、
\mathcal{R}(SL_{2}/GL_{1}) \simeq \mathbf{C}[y_0,y_1,y_2]/(y_0^2+y_1^2+y_2^2-1)
が分かり、これはSU(2)/U(1) \simeq S^2という同型の複素化に相当する。


【超余談】
三角形の合同条件と不変式論
https://m-a-o.hatenablog.com/entry/20130104/p1
の計算は、\mathcal{R}(GL_{2}/SO(2))を決定する問題になっている。三角形のモジュライ空間自体は、鏡映で移るものを区別するなら、GL(2,\mathbf{R})/SO(2)で、同一視する場合は、GL(2,\mathbf{R})/O(2)になり、いずれにせよコンパクトではないんだけど、複素化すると、上のHopf fibrationの計算と殆ど同じことになる。

\mathbf{C}[GL_{2}] \simeq \mathbf{C}[a,b,c,d,1/(ad-bc)]で、ad-bcは、SO(2)作用で不変(三角形のモジュライという観点からは、面積が回転不変ということ)なので、結局、主要な仕事は、\mathbf{C}[a,b,c,d]^{SO_2}の決定になる。この不変式環の生成元は、a^2+b^2,c^2+d^2,ac+bdとなり、(a^2+b^2)(c^2+d^2) = (ac+bd)^2 +(ac-bd)^2という関係式が成立する。

ところで、平面上の三角形は平面上の並進自由度も持つことを考えると、三角形のモジュライ空間は、GA(2)/SE(2)あるいは、GA(2)/E(2)でもある((GA(n)はアフィン変換群で、SE(n),E(n)はEuclid変換群)けど、GA(n)/SE(n) \simeq GL(n,\mathbf{R})/SO(n)、GA(n)/E(n) \simeq GL(n,\mathbf{R})/O(n)なので、簡約群の計算で済む。Riemann計量は、各点でGL(n,\mathbf{R})/O(n)に値を取る場であるが、GA(n)/E(n)に値を取る場と思うこともできる。一般相対論の場合、直交群はローレンツ群に、Euclid群はPoinare群に変更すれば、同じ同型があって、Einsteinによる定式化は前者に基づくけど、後者に基づく理論は、Einstein-Cartan理論とか呼ばれる。このCartanは、Cartan-Helgasonの定理のCartanと同じ人。


複素簡約群GとKに対して、G/Kがアフィン代数多様体であるというのは、松島の定理の簡単な部分である。考えてみると、コンパクト対称空間を複素化すると、アフィンになるというのは、ちょっと面白いことだと思う。複素射影空間や複素グラスマン多様体を、複素化するというのは、ぱっと見奇妙な感じがしないでもないけど、仕方ない。

【松島の定理】松島の定理は、複素簡約群Gに対して、G/Kがアフィンであることと、Kが複素簡約群であることが同値という定理。この定理は、1960年に発表され、色んな証明があるようである。2000年以降でも、以下の論文を見つけた
Invariant Ideals and Matsushima's Criterion
https://arxiv.org/abs/math/0506430

The Kempf-Ness theorem and Invariant Theory
https://arxiv.org/abs/math/0605756


別に定理を持ちださなくても、その座標環の計算は、上で見たように、不変式の計算で原理的には出来る。SL_{n}/SO_{n}やSL_{2n}/Sp_{n}のような場合は、極分解や、そのvariationで、行列の分解に帰着して計算できるけど、一般的には、それだけで済むのか分からない。21世紀になって、DerksenやKemperによって、具体的なアルゴリズムも与えられている(※)。群が少し大きくなると、計算が終わらないと思うけど。

※)計算不変式論の使い方:合同変換群の多項式不変量の計算事例
https://m-a-o.hatenablog.com/entry/20131227/p2

【行列の分解による不変式の計算】例えば、non-oriented Grassmann多様体の場合でも、極分解と似たような方法で、不変式を構成できる。実Grassmann多様体O(p+q)/(O(p)\times O(q))を例に取ると、行列をブロック化して
\begin{pmatrix} A & B\\ C & D \end{pmatrix} \begin{pmatrix} P & 0 \\ 0 & Q \end{pmatrix} = \begin{pmatrix} AP & BQ \\ CP & DQ \end{pmatrix}
なので、AA^{T}やCA^{T},CC^{T},\cdotsなんかの行列成分に不変式が出る。問題は、これらで、不変式の生成元が尽きるのかということだけど、ちゃんと数えてないので不明。まぁ、何となく、他にないやろうという気はするけど


そういうわけで、\mathcal{R}(G/K)は、自明でない環構造を持ってるのだけど、Cartan-Helgasonの定理の証明には、環構造を具体的に知る必要は全然ないし、環構造があることすら重要ではない。これは、Cartan-Helgasonの定理の解析的バージョンがL^2関数空間の分解として書かれることを思えば、当然ではある。GoodmanとWallachの本にある戦略は、\mathbf{C}[G]に対する(代数的)Peter-Weylの定理から、Gの有限次元既約表現をKに制限して、Kの一次元表現が何個出るか見ればいいよね的なものである

まぁしかし、せっかく環構造があるのだから、もうちょっと、そこに注目してみてもいいのじゃないかと漠然と思う。別に、それで何か新しい物が出てくるとか強気に主張できないのが、苦しいとこではあるけれど。

逆に、Cartan-Helgasonの定理による既約分解から、何とかして環構造を復元できるかということを考えると、よく分からない。球面調和関数の空間は、対称トレースレステンソルの空間と同一視できるけど、EastwoodのCartan productという(結合的かつ可換な)積構造を入れることができ、この環は、(複素化した)球面上の座標環とは異なる。

The Cartan Product
https://doi.org/10.36045/bbms/1110205624

尤も、この場合、違いは、r^2=1とするか、r^2=0とするかなので、この2つの環は一方から他方へ連続変形できる。他の場合は、どうなってるのかというのは謎



ところで、Cartanは、Schoutenと共同で、絶対平行性を持つRiemann多様体の分類を調べる中で、曲率テンソルの共変微分が0になる空間(今では局所対称空間と呼ばれる)に関心を持って、対称空間を調べ始めたようだ。例えば、定曲率空間は、局所対称空間なので、定曲率空間を大幅に一般化したと思うことができる。Schoutenとの論文が1926年出版で、この時点では調和解析の片鱗もない。

元々の動機は、今となっては割とどうでもいいものになってしまったけれど、その後、調和解析をやるだけなら、もうちょっと条件の緩い空間でもいいんじゃねと思ったのか何なのか、20世紀中盤〜後半にかけて、
・弱対称空間(Selberg 1956)
・Gelfand pair(Gelfand 1960s?)
・spherical pair(Vinberg&Kimel'feld 1978, Krämer 1979)
・spherical homogeneous space(Brion 1986)
などの概念が現れた。これらは同じようなもので、別に、互いの仕事を認識してなかったというわけでもないようだけど、なんか色々な名前が付いててめんどくさい。

Gelfand pairの心は、不変微分作用素環が可換というものだと思う。Selbergは、弱対称空間が、この性質を持つことを示したそうだ(論文は見てないけど)。Selbergの論文は、Poisson和公式の一般化/Selberg trace formulaにあったそうだ。多分、Selbergの仕事が先っぽいけど、Gelfand pairは、有限群とかも対象になる点で、より一般的なものではある。このへんは、解析学らしさが感じられる話。

一方、spherical pairやspherical homogeneous spaceとかは、なんか関数空間を、G作用で既約分解すると、multiplicity-freeになるという感じのもので、Cartan-Helgasonの定理から、対称空間が、この性質を持つことが分かる。


spherical homogeneous spaceという名前を、最初に使ったのは、Brionという人だと思われる。
Quelques proprietes des espaces homogenes spheriques
https://link.springer.com/article/10.1007/BF01168684

もうちょっと一般化したspherical varietyという用語も、この人が命名したのが、多分、最初っぽい
Spherical Varieties an Introduction
https://link.springer.com/chapter/10.1007/978-1-4612-3702-0_3

私は全然知らなかったけど、spherical varietyの研究は(論文が沢山出てくるので)結構、流行ってるっぽい。このへんの話は、調和解析も表現論も離れて、また違う動機の元で進んでるように見えるので、よく分からない。


で、不変微分作用素環の可換性と、関数空間をG作用で既約分解したらmultiplicity-freeになるのが同値という事実は、色んなバリエーションがある。今の話に関係しそうなとこだと
Multiplicity-free spaces
https://doi.org/10.4310/jdg/1214438422

Compact weakly symmetric spaces and spherical pairs
https://arxiv.org/abs/math/9808039

Weakly symmetric spaces and spherical varieties
https://link.springer.com/article/10.1007/BF01236659

とりあえず、不変式環を見ると、(原理的には)計算できる！と思って安心する

Siegelモジュラー形式環と計算機代数

asir

Rungeという人の1993年の論文

On Siegel modular forms. Part I.
https://doi.org/10.1515/crll.1993.436.57

に、Siegelモジュラー形式環の構造が、有限群の不変式の計算と、環の商体内での正規化(整閉包)によって計算できるという結果が書いてある(Theorem3.12)。有限群が作用する環は、テータ関数(テータ零値)の多項式たちで生成され、単純な多項式環というわけではない。偶数ウェイトのSiegelモジュラー形式環は、もう少し、簡単な計算でいける(Corollary3.17)。


Siegelモジュラー形式は1930年代終わり頃、Siegelが導入した。Siegelモジュラー形式に関する概略が、以下の論説の1〜2節などに書かれてる。
コンパクト化の今昔
https://www.jstage.jst.go.jp/article/sugaku1947/51/2/51_2_129/_article/-char/ja/

歴史的に見ると、(19世紀頃の)初期のモジュラー形式論は、個々のモジュラー形式を理解しようという特殊関数論的(というか古典解析的)な動機が強かったように見える。一方で、Siegelモジュラー形式環を調べる動機の一つは、離散群による商空間という扱い難い対象を理解しようというところにあって、解析的側面は薄れてて、代数幾何学的な話になっている。

【モジュライ空間史】On the early history of moduli and Teichmüller spaces
https://arxiv.org/abs/1602.07208
を読んだところ、Siegelは、1935年の論文で、Torelli空間とSiegel上半平面、Siegel基本領域を定義し、Torelli空間から、Siegel上半平面への複素構造を保つ写像を定義したということが書いてあった。Torelli空間は、おそらく、Weilによる命名で、今では、タイヒミュラー空間のTorelli群による商と定義されることが多い。タイヒミュラー空間の論文は、この少し後に書かれたらしいけど、TeichmüllerがSiegelの仕事を知ってたのかは謎。論文は読んでたとしても、当時、Teichmüllerは、ナチの反ユダヤ主義を支持し、Siegelは、ユダヤ人だったのかは知らないけど、1940年にアメリカに亡命してるので、仲良く交流したとは思えない。それより前の時代だと、Schottkyなんかは、純粋に解析的なアプローチを取っていたようだし、1913年にTorelliの定理を示したTorelliも、周期行列の集合に良い幾何構造が入るという認識は持ってなかったっぽい。1920年代に、対称空間の分類とかが行われてて、Siegelに、その影響があったかは不明。モジュライ空間を厳密に定式化するという点で、Siegelは大きな貢献があったようである



私は、コンパクト化の議論とか読んでも、そんなんええから、Siegelモジュラー形式環の生成元と関係式を書いてくれへんかなという気持ちになる。現在、種数3までのSiegelモジュラー形式環の構造が決定できているっぽい。種数1のSiegelモジュラー形式環は、19世紀には、既に知られていた通り、4次と6次のEisenstein級数で生成される。この環と、有限群の不変式環との間に自然な同型が存在するというのは、誰が最初に気付いたのか分からないけど、有限群の不変式環とある種のモジュラー形式環に自然な同型があるという現象は、1970年頃には知られてたっぽい(もっと古い起源があるのかは不明)。

例えば、1972年のフランス語の論文

Polynômes des poids de certains codes et fonctions thêta de certains réseaux
http://www.numdam.org/item/ASENS_1972_4_5_1_157_0/

には、位数192の群の不変式環が、\mathbf{C}[E_4,\Delta]と自然に対応するということが書かれてるように見える(フランス語は読めない)。

そういう話とは関係なく、Siegelモジュラー形式環を頑張って決めようとした人もおり、種数2のSiegelモジュラー形式環は、井草(1962)
On Siegel modular forms of genus two
https://www.jstor.org/stable/2372812

種数3のSiegelモジュラー形式環は、露峰(1986)の結果がある。
On Siegel Modular Forms of Degree Three
https://www.jstor.org/stable/2374517?seq=1

種数3の場合は、計算が大変なようで、生成元がminimal setかどうか分かってなかったらしいけど、2019年に、minimal setを決めたと主張している論文が出ているよう(内容は未精査)。19個の生成元と55個の関係式があるらしいので、結果をトレースするだけでも、手計算は辛そう

Siegel modular forms of degree three and invariants of ternary quartics
https://arxiv.org/abs/1907.07431


Rungeの結果を使って、モジュラー形式環を決定するには、テータ零値の間に存在する代数的関係式を決める必要があるけど、一般的な方法は知られてないっぽい。m=(m',m'')\in \{0,1\}^g \times \{0,1\}^gに対して
\theta_m(z,\tau) = \displaystyle \sum_{x \in \mathbf{Z}^g} \exp 2\pi i\left(\dfrac{1}{2} {}^{t}[x+\dfrac{1}{2}m'] \tau [x+\dfrac{1}{2}m'] + \langle x+ \dfrac{1}{2}m',z+\dfrac{1}{2}m'' \rangle \right)
とすると、zについて偶関数であるものが、2^{g-1}(2^g+1)個ある。z=0の時しか考えないので、oddな方は、どうでもいい。Rungeの論文では、テータ関数の二乗が、
f_a(\tau) = \theta_{(a,0)}(0,2\tau) = \displaystyle \sum_{x \in \mathbf{z}^g} \exp 2\pi i \left( {}^{t}[x+\dfrac{a}{2}] \tau[x+\dfrac{a}{2}] \right)
を使って書けることから、この形のテータ零値の作る環を中心に考える。この形のテータ零値は、どれも0でなく、2^g個の線形独立な関数となっている(Rungeの論文p59)

テータ関数やテータ零値には、沢山の恒等式があって、関係式の間の関係式が必要なほど、ややこしい。必要なテータ零値の間の関係式も、知られている関係式から得られるものもある。Rungeの論文§6では、種数3の場合に、既知のテータ関係式から、f_aの関係式を得て、Siegelモジュラー形式環の構造を決定している(Theorem6.2)

cf)A remark on a theorem of Runge
https://link.springer.com/article/10.1007/s000130050220
[PDF]https://www.mathi.uni-heidelberg.de/~freitag/preprints/runge.pdf
の式(2)も参照


もっと高い種数では、既知のテータ恒等式に帰着できない関係式が出るのか、よく分からない。

A theta relation in genus 4
https://projecteuclid.org/euclid.nmj/1114631553
では、種数4の場合に、一つの関係式を見つけて、"It is a natural question, whether our relationis a consequenceof Riemann's theta relations. We could not decide this."と書いてる(2001年付けだから、少し前だけど、凄く古いわけではない)。全部が既知のテータ恒等式に帰着できるなら、原理的には、変数消去の計算で、全ての関係式を見つけることができる。実際上は、evenなテータ関数が2^{g-1}(2^g+1)個あって、(τではなく)2τの関数であるf_aが、2^g個あるので、単純に考えると、g=4で152変数、g=5で560変数の多項式の変数消去計算をやらないといけないことになって、現在のコンピュータとアルゴリズムで計算できるかは怪しいけど

【参考文献】arXiv:1401.5368やarXiv:1510.02699には、色々なテータ恒等式の発見に、誰が寄与したかという話が少し書いてある

【テータ恒等式とPlucker関係式】2x4行列の2x2小行列式D_{ij} (1 \leq i \lt j \leq 4)で生成される環は、多項式環を、一個のPlucker関係式で生成されるイデアルで割ったものになる。
\mathbf{C}[D_{12},D_{13},D_{14},D_{23},D_{24},D_{34}] \simeq \mathbf{C}[X_{12},X_{13},X_{14},X_{23},X_{24},X_{34}]/(X_{12}X_{34}-X_{13}X_{24}+X_{14}X_{23})
小行列式D_{ij}に対して、4x4対角行列h \in H_4 \subset GL(4)の作用を
\rho(h)(D_{ij})(x) =D_{ij}(xh)
で定める。\rho(h)(D_{12}D_{34}) = det(h) D_{12}D_{34}etc.などが成立するけど、h \in H_4の作用で、det(h)の多項式倍だけ変化する元の全体は環をなし、
\mathbf{C}[D_{12}D_{34} , D_{13}D_{24} , D_{14}D_{23}] \simeq \mathbf{C}[y_0,y_1,y_2]/(y_0-y_1+y_2)
のようになる。一般の群と指標に対して、この手の環は、半不変式環と呼ばれてることがある。一方、種数1のテータ零値には
\theta_{00}^4 = \theta_{10}^4 + \theta_{01}^4
という関係式があるので、環同型\mathbf{C}[\theta_{00}^4,\theta_{10}^4,\theta_{01}^4] \simeq \mathbf{C}[D_{12}D_{34} , D_{13}D_{24} , D_{14}D_{23}]が成立する。

代数的に書くと、大した意味はなさそうだけど、これは偶然でなく、19世紀から知られる超幾何関数と楕円関数の不思議な繋がりが背景にある。\mathbf{C}[\theta_{00}^4,\theta_{01}^4,\theta_{10}^4]は、モジュラー群\Gamma(2)のモジュラー形式環M(\Gamma(2))になっている。このモジュラー形式環は、楕円曲線の単位元以外の等分点e_1,e_2,e_3を使うと、
M(\Gamma(2)) = \mathbf{C}[e_1 , e_2 , e_3]
のようにも書ける。テータ零値で書くと
e_1=-\dfrac{1}{3}(\theta_{00}^4 + \theta_{10}^4)
e_2=\dfrac{1}{3}(\theta_{00}^4 + \theta_{01}^4)
e_3=-\dfrac{1}{3}(\theta_{01}^4 - \theta_{10}^4)
で、e_1+e_2+e_3=0が成立する。SL(2,\mathbf{Z})/\Gamma(2) \simeq S_3は、e_1,e_2,e_3の置換として作用し、M(SL(2,\mathbf{Z})) \simeq M(\Gamma(2))^{S_3}で、M(SL(2,\mathbf{Z}))の生成元として、e_1e_2+e_2e_3+e_3e_1とe_1e_2e_3が取れて、それぞれ、定数倍を除いてEisenstein級数と一致する(E_4=-3(e_1e_2+e_2e_3+e_3e_1) , E_6 = \dfrac{27}{2}e_1e_2e_3)(逆に、M(\Gamma(2))はM(SL(2,\mathbf{Z}))上の自由加群になってて、この構造を具体的に書いてみるのも面白い)。一方、\mathbf{C}[D_{12},D_{13},D_{14},D_{23},D_{24},D_{34}]のProjは、Grassmann多様体Gr(2,4)で、これをH_4で"簡約"した空間が、Gauss超幾何関数の住む空間と思える。

"私説超幾何関数"(ISBN:978-4-320-01576-0)という本を読むと、種数2でも類似の同型が作れると書いてある。今度は、Grassmann多様体Gr(3,6)を考える。3x6行列の3x3小行列式D_{ijk}(1 \leq i < j < k \leq 6)は、20個あって、代数的トーラスのdetに関する半不変式環Z(3,6)は、10個の生成元と5個の独立な関係式からなる。種数2のテータ零値は10個あるけど、テータ零値の2次式を10個うまく選ぶと、それらの二乗から生成される環(適当なモジュラー群のモジュラー形式環となるらしい)と半不変式環Z(3,6)とが同型になるようだ。種数3以上に一般化する話は知られてないっぽいし、根源的な理由も分からないので、モジュラー形式というより、テータ零値の恒等式(の一部)に、何か良い構造があるという可能性もあるんじゃないかと勝手に思ってる



それ以外のステップ。有限群の不変式環を計算するアルゴリズムは以前に実装した

有限群の不変式環の生成元を計算するアルゴリズム
https://m-a-o.hatenablog.com/entry/20130316/p2

asir_misc/finvar.rr
https://github.com/vertexoperator/asir_misc/blob/master/finvar.rr

にも、Risa/Asirで実装したものが置いてある。しかしまぁ、過去の私の記述によれば、"有限群の不変式の計算は、SingularやMagamaにも実装されてるらしいので、計算するだけなら、それ使うのが一番確実かもしれない"らしいので、Singular使えばいい気がする。

正規化の方は、Risa/Asirには、ないので、使いたければ自分で実装するしかない(と思う)けど、Singularには、既に実装されている

D.4.18 normal_lib
https://www.singular.uni-kl.de/Manual/4-0-3/sing_1269.htm#SEC1344

つまり、総合的に見て、Singular使えばいいのでは。一応、Singularで使ってる正規化アルゴリズムは、以下の論文に書かれているようだ。

Normalization of Rings
https://arxiv.org/abs/0904.3561

以前測定した限りでは、グレブナー基底の計算は、Singularの方が、Risa/Asirより、速い(数倍くらい違うこともあった気がする)。ただ、Risa/Asirは、何故か、代数的数のサポートが割と手厚い気がするので捨てがたい。あと、Singularは、(Mac以外)インストールめんどくさい



種数1の場合は、容易に計算できるし、Rungeの論文では、Remark3.6の直後に、結果が書いてある。種数1の場合は、群がShephard-Toddの(既約擬鏡映)群No.8というものになる

cf)List of irreducible complex reflection groups
https://en.wikipedia.org/wiki/Complex_reflection_group#List_of_irreducible_complex_reflection_groups

【用語について】正方行列Aが擬鏡映とは、Aが有限位数で、rank(I-A)が1であることを指す。英語だと、pseudo-reflection。位数が2の時は、ある超平面に関する鏡映になる。同じ概念に対して、complex reflection,unitary reflectionとかいう呼び方をされることもある。擬鏡映は、条件から、1でない固有値を一つだけ持ち、それは、１のべき根になってるから、こういう名前が付いてるのだと思う。+1でない固有値が実であれば、-1しかありえず、Aは鏡映になり、普通の鏡映は、実鏡映ということになる。個人的には、擬鏡映という名前が一番わかり易い。擬鏡映で生成される群が、擬鏡映群


私の実装でも

Z2 = newalg(x^2+1);
T=newmat(2,2,[[1,1],[1,-1]])*(1+Z2)/2;
D=newmat(2,2,[[1,0],[0,Z2]]);

finvar(grp([T,D]),[f0,f1]);

とかやって、一瞬で正しい答えが出た。
f_0 = \theta_{00}(2\tau) , f_1 = \theta_{10}(2\tau)
で、
f_0^2 = \dfrac{\theta_{00}^2 + \theta_{01}^2}{2}
f_1^2 = \dfrac{\theta_{00}^2 - \theta_{01}^2}{2}
が成立し、 2つの独立な不変式は
g_4 = f_0^8+14f_1^4 f_0^4+f_1^8
g_6 = f_0^{12} - 33f_0^4f_1^4(f_0^4+f_1^4)+ f_1^{12}
となる。Eisenstein級数は
E_4 = \dfrac{1}{2} ( \theta_{00}^8 + \theta_{10}^8 + \theta_{01}^8)
E_6 = \dfrac{1}{2} (-3\theta_{10}^8(\theta_{00}^4+\theta_{01}^4) + \theta_{00}^{12} + \theta_{01}^{12})
なので、両者が一致することは、単なる計算問題


種数2の場合。群H_2の生成元は、整係数Symplectic群の生成元から来ていて、Jから来る方はいいけど、整係数対称行列は、無限個ある。しかし
S = \left( \begin{matrix} a & b \\ b & c \end{matrix} \right)
と置くと
D_S = diag(1 , i^c , i^a , i^{a+2b+c})
になって、a,b,cは、0~3の範囲を舐めれば十分。有限個の生成元が決まったので、群の要素を列挙しようと思ったけど、pythonのsetに相当するものがないRisa/Asirの実装では、まる一日かけても、46080個の元を列挙するのは無理だった。仕方ないので、pythonで列挙して、ファイルから群の要素を読み込むという方法で対応(計算時間は、夜寝ている間に終わる程度)。

Rungeの論文には、これは擬鏡映群だと書いてあるけど、群に含まれる擬鏡映をコード書いて列挙してみると60個ある。位数は全て2で、どうやら、Shephard-Toddの群No.31っぽい(ちゃんとチェックはしてないけど、擬鏡映を列挙できてるので、比較すれば、力ずくで確認できるはず)。Molienの公式でHilbert級数を計算する(これはRisa/Asirで実装してある。昼寝してたら終わってた)と、
\dfrac{1}{(1-t^8)(1-t^{12})(1-t^{20})(1-t^{24})}
になっている。なので、この群を、4変数多項式環に作用させた不変式環は、8,12,20,24次の代数的に独立な不変式で生成される(8*12*20*24=46080で、群の位数が得られる)。具体的な式は、Rungeの論文p77冒頭に記載されているけど、Reynolds作用素の計算を、コンピュータにやってもらえばできる。論文にも書いてあるけど、20次と24次の不変式の空間は、一次元ではないので、定数倍を除いても、論文と同一の式にならない可能性はある。というか、私はならなかった。

モジュラー形式としては、ウェイト4,6,10,12に対応し、(Corollary3.17によって)種数2、偶数ウェイトのSiegelモジュラー形式環は、この4つの元で生成される。奇数ウェイトのモジュラー形式として、ウェイト35のものがあり、この5つの元が、種数2のSiegelモジュラー形式環を生成する


種数3の場合は、Rungeの群の位数が約3.7億なので、全部要素を列挙して計算するというのは望みがない。多分、私のPCでは、メモリに載らない。詳細は、
On Siegel modular forms. II
https://projecteuclid.org/euclid.nmj/1118775400
の方にあって、群の構造を調べて何とかしてるようだ。

現在のアルゴリズムは、Molienの公式やReynolds作用素の直接計算に依存してるので、位数の大きな群に対して使うには、どうしても限界がある。けど、種数が大きくなるごとに試行錯誤してたらキリがないので、超でかい有限群の不変式と戦えるアルゴリズムが欲しい。次元が上がると、今度はグレブナー基底の計算が辛くなるという二重苦が待っている

解析力学から見た古典/量子Kepler系の解法

古典Kepler系の解析解
https://vertexoperator.github.io/2016/12/01/kepler2d_classical.html

2次元量子Kepler系の束縛状態
https://vertexoperator.github.io/2016/12/01/kepler2d_quantum.html

で、2次元のKepler系を、古典力学と量子力学で解いている。多少見かけは違ってるけど、同じ方法で空間変数の取り直しを行って、調和振動子に帰着(全く同じなので、書いてないけど、正エネルギー/散乱状態の場合は、inverted harmonic oscillatorというものに帰着する)、周波数の形も古典/量子で同じ形になる

ところで、古典Kepler系の方は、時間変数も取り直しているのに対して、量子Kepler系では、時間変数の変換は出てきてない。量子系の方は、時間非依存のシュレディンガー方程式を考えているのだから当然だけど、量子Kepler系のHamiltonianに対して、時間依存シュレディンガー方程式を考えて、時間変数の取り換えを行ってから、(新しい時間変数に対する)定常シュレディンガー方程式を考えると、結局、Hamiltonianは、左からrを掛けるという変換を受けることになる。

"量子Kepler系の束縛状態"の解法2の方で、「量子Kepler系のシュレディンガー方程式を
r(\hat{H} - E)\Psi = 0
と置くと・・・」と書いてあるのは、自明な式変形を書いてるのではなく、シュレディンガー方程式をこのように見ることが、古典系に於ける時間変数の取り換えと同等の操作を行うことに相当している。また、spectrum generating algebraに含まれるのは、Hamiltonianそのものではなく、r(\hat{H} - E)であることからも、この変換の必要性が理解できる。

ただ、これだけだと、ややad hocな印象もある



Duru-Kleinertによる、量子Kepler系の経路積分に関する以下の論文には、"pseudo Hamiltonian"というものが書いてある(論文中の式(10))。彼らは、これ以前にも同じような論文を書いてるけど、そっちには、このような説明はないので全然参考にならない

Quantum Mechanics of H-Atom from Path Integrals
https://doi.org/10.1002/prop.19820300802
[PDF]http://users.physik.fu-berlin.de/~kleinert/83/83.pdf

r(\hat{H} -E)は、Duru-Kleinertが"pseudo Hamiltonian"と呼んでいるものに相当する。


で、Duru-Kleinertの論文の式(15)を見ると、このpseudo Hamiltonianが住んでるのは、拡大相空間であることが読み取れ、pseudo Hamiltonianに対して正準方程式を考えると、普通の運動方程式+時間変数の取り換えの式が出ている。Duru-Kleinertは、拡大相空間とは言ってないけど、昔の私が調べた所によると、1949年にLanczosが導入した比較的新しい概念らしい。昔の私が書いてることを信じるなら、H-Eはextended Hamiltonianと呼ばれることもあるそう

cf)extended phase space
https://formalgroup.tumblr.com/post/105356647155/extended-phase-space

古典力学では、H-E=0は定義みたいなものであり、量子力学では、H-Eを波動関数に作用させると0になるというのは、シュレディンガー方程式そのもの。


拡大相空間を考えるなら、(明らかに、拡大相空間も、自然なsymplectic構造を持ってるので)正準変換も"拡大"されるべきと思うのだけど、

Stackel系の全ての保存量を保つ離散化
https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/42747

とかいう関係なさそうな文献を見ると、そのような変換を考えた人はいるらしい。Duru-Kleinertのpseudo Hamiltonianは、extended Hamiltonianに、一般化した正準変換を施したものであると理解できる。

DuruとKleinertは、経路積分のために、pseudo Hamiltonianを考えたけど、経路積分自体は、どうでもよく、背景に、こういう構造があることが重要なのだと思う



古典力学では、ハミルトニアンは関数なので、右から関数を掛けるか左から関数を掛けるか気にする必要はないけど、量子力学では、事情が違う。Kepler系では、r(\hat{H}-E)としたのだけど、(\hat{H}-E)rでない理由があるだろうかという疑問が残っている。一般には
f_{l}(\hat{H}-E)f_{r}
という形の変換が考えられる。この形の変換を取り扱った文献で、一番古いのは、調べた限り

On the generalization of the Duru-Kleinert-propagator transformations
https://link.springer.com/article/10.1007/BF01318171

じゃないかと思う。

Kleinertの本が、Kleinert自身によって公開されているけど、

New Path Integral Formula forSingular Potentials
http://users.physik.fu-berlin.de/~kleinert/b5/psfiles/pthic12.pdf

この形の変換が扱われている。式(12.26)とか。本を読めばわかるけど、新しい時間変数と、時間変数の関係は、f_lとf_rの積のみで決まる。式(12.37)など。


この形の変換で一番良く見るのは、相似変換で、その場合は、時間変数を変える必要はない。相似変換については、当たり前すぎて、注意を払われることもないけど、それより、ずっと一般的な変換を行っても、ツケを時間変数に押し付けてしまえるというのは、言われれば当たり前だけど、何かに使えることもあるかもしれない。Kepler系の場合は、rを掛けることによって、Coulombポテンシャルの特異性が見かけ上消えるという見方もできる


【おまけ】Kleinertは、1967年に、A.O. Barutと共に
Transition Probabilities of the Hydrogen Atom from Noncompact Dynamical Groups
https://doi.org/10.1103/PhysRev.156.1541

という論文を書いていて、3次元量子Kepler系の束縛状態の空間が、so(4,2)の既約表現空間になってることを、(私が知る限り)初めて明確に書いた人でもある。この論文には、まだ"pseudo Hamiltonian"は書かれていない。もう少し後のBarutの論文では、Hamiltonianそのものより、"pseudo Hamiltonian"(とは呼んでないけど)を重視する書き方をしているものがあったと記憶している

Group Theory and Orbital Fluctuations of the Hydrogen Atom
http://users.physik.fu-berlin.de/~kleinert/kleiner_re222/barut92.html

のAbstractで、Kleinertは、以下のように書いている。

The group theoretic development was triggered by A.O. Barut who suggested to me the search for a dynamical group larger than SO(4). In this way he became partly responsible also for the later and most recent path integral development.

Chapman-Jouguet理論

物理で定義とか真面目に考え出すと、辛いことになるけど、"法線方向の速度が不連続に変化する波面の伝播"を衝撃波と総称する立場からすると、発熱反応を伴う圧縮衝撃波がデトネーション波ということになる。発熱反応を伴う衝撃波では、膨張波も存在することができ、こっちは、デフラグレーション波と呼ぶ。日本語では、デトネーションは爆轟と訳され、デフラグレーションは爆燃と訳されている。

※)接線方向の速度が不連続に変化する不連続面は、渦層と呼ばれ、衝撃波面とは異なる。同じ不連続面とはいえ、使われる手法とかも全然別のよう。

普通の衝撃波は、化学反応とかは伴わないので、そのへんがデトネーション/デフラグレーションとの違い。デフラグレーションは、不連続波面を持つけど、伝播速度は音速より遅いので、衝撃波のイメージに、そぐわない。まぁ、あまり深く考えず、気相爆轟=超音速で伝わるガス爆発くらいの認識でいいっぽい。

気相爆轟は、19世紀の終わり頃から研究されるようになり、理論的な研究は、Mikelson(1890),Chapman(1899),Jouguet(1905)などが独立に行ったらしい。これらは、同じ結果を与えるらしく、現在では、Chapman-Jouguet理論(以下、CJ理論)として知られる。現在、CJ理論は、衝撃波のRankine-Hugoniot理論を拡張した形で展開される。

どうでもいいけど、CJ理論のChapmanとChapman-Enskog展開のChapmanは別人。Chapman-Enskog展開の人は、Sydney Chapman(1888-1970)という名前で、CJ理論の人は、David Chapman(1869-1958)らしいけど、David Chapmanには、Sydney John Chapman(1871-1951)という弟がいて経済学者だったらしい。ややこしいわ

CJ理論では、波面の前後が、それぞれ未燃ガス、既燃ガスになっていて、反応は一瞬で終了するものとされている。未燃ガス、既燃ガスは共に、理想気体として記述するのが一般的。基本的に、代数的な計算しか出てこないので、難しいわけでないけど、既燃ガスの比熱比は常温のものと違うことに気付かず、計算してみると、値が合わないという罠にはまった

後になって気付いたけど、デトネーションのことは、ランダウ・リフシッツの流体力学の教科書に載っていて、CJ理論が説明されてる。この本も、測定値と理論値の比較とかはないので、助けにはならなかったろうけど。

Fluid Mechanics §121
https://archive.org/details/FluidMechanics/page/n491



罠にはまった時に、Chapmanの原論文を読んだ。

On the rate of explosion in gases
https://www.tandfonline.com/doi/abs/10.1080/14786449908621243

Chapmanの議論は、現在の定式化とは違っていて分かりにくいけど、TABLE IIIに、"Velocity of Explosion"の式として

V^2 = \dfrac{2RJ}{\mu C_{v}^2} \left( {(m-n)C_p + m C_{v} }C_{p} t_0 + (C_p+C_v)\mathrm{h} \right)

が書かれていて、様々な混合気体の燃焼反応に対して、Vを計算し、測定値と比較している。単位とかが書かれてないけど、私が解読した限りでは、多分、以下のようになっている。

Chapmanの論文には、多くの混合気体に対する測定値と計算値が載ってるけど、全部やるのは大変なので、水素と酸素がモル比で1:1に混合した気体の燃焼反応(2モルの水素と2モルの酸素を燃焼させて、1モルの水蒸気が発生し、1モルの酸素が反応せず残るという状況)
2\mathrm{H}_2 + 2\mathrm{O}_2 = 2 \mathrm{H}_{2}\mathrm{O} + \mathrm{O}_2
を例に取ると、

RJ = 8.314462618(J/mol/K) : 気体定数
μ = 68(g) : 単位量あたりの気体の質量(反応前後で保存されるはずなので、水素2モル+酸素2モルの質量であり、また、水2モルと酸素1モルの質量でもある)
m = 3(mol) : 反応前のモル量(水2モルと酸素1モル)
n = 4(mol) : 反応後のモル量(水素2モルと酸素2モル)
h = 117000(cal) : 反応熱=1molのH2O(g)につき、245kJ弱の計算
t0 = 273.15(K) : 基準温度
Cv = 11.63(cal/mol/K) : 生成ガスの定積モル比熱(表に記載されてる値)
Cp = 13.62(cal/mol/K) : 生成ガスの定圧モル比熱(Cvから、マイヤーの関係式を使って計算)
V = 2341.8(m/s) : Chapman-Jouguet速度(論文では2340と計算されている)

だと思う。


生成ガスの比熱について、論文には、以下のようなexcuseが長々と書いてあって、時代を感じる

A few words are necessary regarding explosions in which water is formed. If the specific heat of steam is taken as 3/2 × specific heat of the diatomic gases, the found rates of explosion fall below the calculated rates when the dilution with inert gas is great, and vice versd when the dilution is small.

It is possible to account for this by two theories. The first theory is that at high temperatures the water is dissociated, whereas at low temperatures the combination of hydrogen and oxygen is complete. The second theory is that the specific heat of steam rises more rapidly with the temperature than the specific heat of the diatomic gas.
(中略)
We are therefore encouraged to test the first theory, i.e. that the specific heat of steam rises more rapidly with the temperature than that of the diatomic gases.

水蒸気の比熱を、2原子分子の1.5倍とすると、うまくいかないので、高温では、水蒸気の比熱は、2原子分子より急速に上昇するのだろうとか。Chapmanが書いている1.5倍というのは、多分、エネルギー等分配則の誤った適用による(比熱が一分子あたり原子数に比例すると考えたんだろう)ものだろう。Chapmanは、もう一つの可能性として、現在でいうところの弱電離プラズマのようなものを考えて却下している。



反応直後のガスは、高圧で、温度は数千度に達するらしいので、常温の場合と違って、水や酸素でも、分子振動の比熱への寄与が無視できなくなる。分子内振動を調和振動子で近似して、独立なモードに分離すると基準振動が得られる。基準振動波数は、赤外分光やラマン分光で測定でき、水素分子の基準振動波数は、ν=4160(/cm)程度。分光業界(?)では、波数の記号は、νの上にチルダを付けたものを使うことが多い気がするけど、以下では、普通のνを使う。回転スペクトルを測るマイクロ波分光はGHzを使ってるっぽいので、赤外/ラマンも、THz単位に揃えろやと思うけど。水の場合、適当に検索した数値によると、1595(/cm),3756(/cm),3657(/cm)。

世の中には合成法が確立してなかったり合成難度が高い化合物もあるけど、その場合でも(低分子なら)割と高速かつ、それなりの精度で、基準振動波数を数値計算することもできる。例えば、現時点で、合成法が確立してない低分子の一つに、テトラヘドランC4H4があり、分子量などからすると、常温常圧で気体の可能性がある(構造異性体の1,2−シクロブタジエンや、1,3-シクロブタジエンの沸点が40度前後という情報がある)

cf)Tetrahedrane
https://en.wikipedia.org/wiki/Tetrahedrane

振動モードが一個の場合、分子振動の分配関数は
q_{\mathit{vib}} = \left( \dfrac{e^{-\beta c \mathit{h} \nu/2}}{1 - e^{-\beta c \mathit{h} \nu}}　\right)^N
なので、一個の振動自由度による定積比熱への寄与は
C_{v,\mathit{vib}} = N k_{B} \dfrac{x^2 e^{x}}{(e^{x}-1)^2}
と計算できる。但し、x = \dfrac{c h \nu}{k_{B}T}

ν=4160(/cm)、T=300(K)だと、x=19.95で、y=\dfrac{ C_{v,\mathit{vib}} }{N k_{B}}は、8.62e-7とほぼ無視できる。T=3000(K)だと、x=1.995で、y=0.725になって、無視できない寄与があることが分かる。

十分高い温度では、振動の自由度1つにつき、モル比熱は気体定数Rだけ増加する。振動自由度を高温近似で扱うと、高温の二原子分子理想気体では、定積比熱が3.5R,定圧比熱が4.5Rになり、高温の非直線型三原子分子理想気体では、定積比熱は6R,定圧比熱は7Rになる(はず)。

6R=50(J/mol/K)=12(cal/mol/K)程度で、上のChapmanの計算に出ている11.63(cal/mol/K)に近い値になっている。実際に上の計算で扱ってるのは、水蒸気2モルと酸素1モルの混合気体の定積比熱なので、それぞれ理想気体として、高温近似が成立するならば、定積比熱は、10.27(cal/mol/K)程度という計算。完全に一致とは言えないけど、それなりに合ってそうではある


100年以上前の論文の数値と理論値だけを頼りにするのもどうかと思うので、もうちょっと新しいデータを見ておく(数千度とかで比熱を測定できるのか知らないけど、以下のデータは、測定値というわけではない)

Water@NIST Chemistry Webbook
https://webbook.nist.gov/cgi/cbook.cgi?ID=C7732185&Units=SI&Mask=1&Type=JANAFG&Table=on#JANAFG

を見ると、水の6000(K)に於けるCpの値は、60.59(J/mol/K)となっていて、これは、7R=58.20(J/mol/K)と近い。4100(K)では、丁度、58.20(J/mol/K)となっている。二原子分子気体の例として、H2(g)を見ると

Hydrogen@NIST Chemistry Webbook
https://webbook.nist.gov/cgi/cbook.cgi?ID=C1333740&Units=SI&Mask=1&Type=JANAFG&Table=on#JANAFG

に、3200(K)で、Cp=37.53(J/mol/K)とある。これは、4.5R=37.415(J/mol/K)と、ほぼ一致している。上の基準振動の数値に基づく計算を考慮すると、3200(K)で、この数値に達するのちょっと早くね？という気もするけど。6000(K)に於けるCpの値は、41.97(J/mol/K)となっている。

この手の小さな多原子分子では、3000〜4000(K)で、比熱への振動自由度の寄与を高温近似で扱うことが許されそうにも思える。幸い、デトネーションに於ける生成ガスの温度は、この水準にある

これらのデータが、どれくらい信用できるものなのか分からないけど、並進、回転、分子内振動を単純に考慮しただけのモデルとは少しズレがあるけど、実験値と比較しないことには、何も言えそうにない。


【おまけ】大学一年の時、エネルギー等分配則による気体の比熱の説明を習った際、結構ずれてる気体もあって胡散臭いと思った記憶がある(条件とか分子は忘れた)。ファインマンの教科書の

The Feynmann Lectures on Physics　Vol.I 40-6 The failure of classical physics
http://www.feynmanlectures.caltech.edu/I_40.html#Ch40-S5

を見ると、古典統計力学で気体の比熱を説明できないことは、丁度Chapmanの論文が書かれた時代の物理学者も不可解に思ってたと書いてある。そんな説明を21世紀にもなって教えないで欲しかった。ファインマンは、分子振動を調和近似して量子統計力学で扱えばいいと言ってるけど、計算例は一つもないし、他で見たこともないので、やってみたのが、以下の表。

分子 Cp(J/mol/K) Fr Cp/R-2.5-Fr/2 基準振動(/cm) X R*(2.5+Fr/2+X)
N2(g) 29.124 2 0.0028 2358.57 0.001477 29.113
O2(g) 29.376 2 0.0331 1580.19 0.028395 29.337
HCl(g) 29.136 2 0.0043 2990.95 0.000112 29.102
HI(g) 29.156 2 0.0067 2309.01 0.001798 29.116
Br2(g) 36.048 2 0.8356 325.32 0.8177 35.899
CO2(g) 37.129 2 0.9656 [667,667,1388,2349] 0.9565 37.053
N2O(g) 38.617 2 1.1446 [589,589,1285,2224] 1.1428 38.603
NO2(g) 36.974 3 0.4470 [750,1320,1617] 0.4653 37.127
O3(g) 39.238 3 0.7192 [719.6 ,1062.1 ,1157.8] 0.6745 38.866
H2S(g) 34.192 3 0.1124 [1183,2615,2626] 0.1098 34.171
NH3(g) 35.652 3 0.2880 [950,1627,1627,3337,3444,3444] 0.2671 35.479
ND3(g) 38.225 3 0.5974 [793,1225,1225,2495,2652,2652] 0.5254 37.626
PH3(g) 37.102 3 0.4623 [992.1,1118.3,1118.3,2322.9,2327.7,2327.7] 0.4656 37.130
BF3(g) 50.446 3 2.0673 [480.7,480.7,696.7,888,1463.3,1463.3] 2.0673 50.446
C2H2(g) 44.095 3 1.3034 [612,729,729,1974,3289,3374] 1.2929 44.077
CH4(g) 35.639 3 0.2864 [1306,1306,1306,1534,1534,2917,3019,3019] 0.2864 35.639
CH3Cl(g) 40.731 3 0.8988 [732,1017,1017,1355,1452,1452,2937,3039,3039] 0.8998 40.739
SF6(g) 96.960 3 7.6616 [346,346,346,524,524,524,...](略) 7.6648 96.987
C3H8(g) 73.60 3 4.8520 [216, 268, 369, 748,...](略) 4.524 70.874

定圧モル比熱Cpの出典は、JANAF熱化学データ表かNIST Chemistry WebBookで、298.15(K) , 0.1(MPa)の値。Frは回転の自由度で、Rは気体定数。ND3は、アンモニアの水素を全部重水素に置換したもの。

表の値Xは、全基準振動波数について、以下の和を計算したもの(温度Tは、298.15K)
X = \displaystyle \sum_{i} f(\dfrac{c h \nu_i}{k_{B}T})
f(x) = \dfrac{x^2 e^{x}}{(e^{x}-1)^2}
並進と回転運動からのみ比熱への寄与がある理想気体を高温近似で考えると、Cp=(5＋Fr)R/2になる。分子内振動からの比熱への寄与を量子統計力学に基づいて計算し、加えると、Cp=(2.5+Fr/2+X)*Rになる。これを測定値と比較すればいい

基準振動の出典は、以下のどれか
(1)JANAF熱化学データ表
(2)Vibrational Modes in HITRAN2012 Database
https://www.cfa.harvard.edu/hitran/vibrational.html
(3)Tables of molecular vibrational frequencies, consolidated volume I (Takehiko Shimanouchi)
https://nvlpubs.nist.gov/nistpubs/Legacy/NSRDS/nbsnsrds39.pdf
https://archive.org/details/tablesofmolecula39shim

(2)のアセチレンのデータは、縮退度にミスがあるので修正してある。また、六フッ化硫黄の基準振動は(2)から、
[346,346,346,524,524,524,615,615,615,643,643,775,948,948,948]
を使い、27個あるプロパンの基準振動は、(3)に従って、以下の数値を使った
[216,268,369,748,869,922,940,1054,1158,1192,1278,1338,1378,1392,1451,1462,1464,1472,1476,2887,2887,2962,2967,2968,2968,2973,2977]

計算値は、これなら、過去の自分も納得しただろうと思う程度には合ってる。どっちみち、調和近似が完全に正確な結果を与えるはずはないので、多少のずれは、あって然るべきもの



Chapmanの論文にあるCJ速度の式は、実際の数値を見ると、大抵の場合、第二項が卓越しており、更にマイヤーの関係式を使うと
V^2 \approx \dfrac{2(C_p - C_v)}{\mu C_v^2} (C_p + C_v) h = \dfrac{2(C_p^2-C_v^2)}{C_v^2} \dfrac{h}{\mu}
と書ける。この近似式は、現代の標準的なCJ理論でも導くことができる。冒頭のランダウ・リフシッツの教科書だと式(121.9)

水素2モルと酸素1モルの爆発の場合、生成ガスは、水蒸気2モルのみで、
h = 242(kJ/mol)
μ = 18(g/mol)
を使うことにする。更に、高温の非直線型三原子分子理想気体に対するCp=7R,Cv=6Rを使うと、CJ速度は
V = 3116(m/s)
となる。

デトネーション伝播の基礎
https://www.jstage.jst.go.jp/article/jcombsj/55/174/55_317/_pdf/-char/ja

の表1を見ると、CJ速度は2841(m/s)で、"既燃ガスの実効的な比熱比"は、1.129となっている。観測されるデトネーション速度も、こんなものらしい


上で使用した反応熱は298K,1atmのもので、0K換算の値を使うべきと思う。0Kに於ける生成エンタルピーは、定圧比熱を、0(K)から298(K)まで積分した値を引けば得られる。面倒なので、(振動の影響を無視した常温での)定圧比熱と基準温度の積で見積もることにすると、H2(g)の場合、定圧比熱が3.5Rとすれば、8.672(kJ/mol)くらい引けばいい。O2(g)とH2O(g)に対して、同様に熱補正を考慮すると、0(K)に於ける反応エンタルピーは、上の値から3.1(kJ/mol)ほど引いたものになる。まぁ影響は小さい


もうひとつ、上の近似式は、反応ガス、生成ガスが共に理想気体と仮定されていて、定圧比熱が温度に依らず一定であることは重要。もうちょっとよく考えると、CJ理論のRankine-Hugoniotの式を使う定式化では、波面前後でのエネルギー保存則が、
h_1(T_1) + \dfrac{u_1^2}{2} = h_2(T_2) + \dfrac{u_2^2}{2}
と書かれる。添字の1と2は、反応ガスと生成ガスの物理量を表し、h_i,u_i,T_iは、比エンタルピーと流速、温度。q = h_1(0)-h_2(0)が、0Kに於ける(単位質量当たりの)反応熱(Chapmanの論文の記号に従うと、h/μに相当する量)。で、生成ガスが理想気体と仮定すれば
h_2(T_2) = h_2(0) + \dfrac{\gamma_2}{\gamma_2-1} \dfrac{N k_B T}{\mu}
となる。\gamma_2は生成ガスの比熱比。\mu,Nは、反応前後の質量と、反応後の分子数。

T_2は高温であるが、単原子分子気体の場合を除けば、分子内振動の影響で、比熱は全温度領域で一定でありえない。単原子分子気体は反応性が低く、デトネーションの文脈では多原子分気体が、いつでも重要だろう。多原子分気体では、上のエンタルピーの式を、そのままで使うのは、どう考えても良くない

考えられる処方として
C = \lim_{T_2 \to \infty} \left( h_2(T_2) - h_2(0) - \dfrac{\gamma_2}{\gamma_2-1} \dfrac{N k_B T}{\mu} \right)
を定義して、高温領域で
h_2(T_2) \approx (h_2(0) + C) + \dfrac{\gamma_2}{\gamma_2-1} \dfrac{N k_B T}{\mu}
と近似すれば、q = h_1(0)-h_2(0)-Cと取りなおすだけで、Chapman-Jouguet速度に関する解析的な式を使うことも許されるはず

少し手を抜いて、並進と回転成分については、全領域で高温近似が妥当だとみなし、振動成分からくる比熱のみ、量子統計力学と古典統計力学で、差があると近似する。基準振動モードが一つで、波数νの場合は
C = \dfrac{N k_B T}{\mu} \displaystyle \int_{0}^{\infty}(f(\dfrac{c h \nu}{k_B T}) - 1)dT
f(x)=\dfrac{x^2 e^{x}}{(e^{x}-1)^2}
で、積分を実行すると
C = -\dfrac{N}{\mu} \dfrac{c h \nu}{2}
となる。基準振動モードが複数ある場合は足すだけ。H2O(g)の場合、1595(/cm),3756(/cm),3657(/cm)を使うと

Cμ = -6.02214076e23 * 2.99792458e8 * 6.62607015e-34 *(1595e2+3756e2+3657e2)*0.5 / 1000.0 = -53.88(kJ/mol)

となる。符号が負なので、CJ速度は、より大きく評価されるようになり、ますます測定値と合わなくなる(ダメ)けど、少なくとも、Chapmanの原論文の意図からすると、このような数値を使うのが正しい気がする。Cの値は熱力学データから決まらない(比熱の温度依存性からフィッティングを行えば、基準振動を知ることができるかもしれないけど、普通、そんなことはしない)せいか、普通、教科書を読んでも、こうは書いてない

まだ、生成ガスの比熱比の値が、考慮に入れてない効果によって、7/6より小さいかもしれないという希望は残ってるけれど


Lasing action and the relative populations of vibrationally excited carbon monoxide produced in pulse-discharged carbon disulfide-oxygen-helium mixtures
https://pubs.acs.org/doi/abs/10.1021/j100639a019
では、CS2-O2-He混合ガス火炎中に生成するCOの振動準位の分布が、点火後の経過時間で、どう変化するか計測している。論文Figure9によると、点火後、数十μsまでは、ボルツマン分布と異なり、特定の振動準位にピークがあり、200μs以上かけて、ボルツマン分布に近付くらしい。デトネーションでも、平衡に達するのに、同じような時間スケールが必要だとすると、例えば、速度3000(m/s)で、200μsの間に、60cm進む。これは、平均自由行程より遥かに長いから、化学反応が一瞬で終わって、短時間で平衡に達するという近似の方が問題かもしれない。


【補足】Chapman-Jouguet理論で必要な熱力学パラメータは、未知化合物でも、量子化学の力によって原理的には計算できる。0(K)に於ける生成エンタルピーは、電子エネルギー+分子内振動の零点エネルギーで計算すればいい。電子エネルギーは、量子化学計算で、ほぼ必ず計算する量で、電子系のハミルトニアンの最低固有値に相当する。分子内振動の零点エネルギーは、分子内振動を調和振動子で近似して最低エネルギーを見る。例えば、H2(g)の場合、基準振動波数ν=4160(/cm)とすれば、
0.5*6.02214086e23*2.99792458e8*6.62607004e-34*4160e2/1000 = 24.88(kJ/mol)
くらいの影響がある。実際上の問題としては、反応熱には、複数の分子の電子エネルギーの差が寄与することになるけど、通常、その大きさは電子エネルギーと比較すると非常に小さい。従って、電子エネルギー自体が1%程度の誤差で計算できても、反応熱の誤差は、割と大きくなり得る　

higher genus AGMは超幾何級数表示を持つか

a,bの算術幾何平均を、M(a,b)として、FがGauss超幾何関数の時

\dfrac{1}{M(1,x)} = F(\dfrac{1}{2} , \dfrac{1}{2},1 ; 1-x^2)

であることは、よく知られてる。


19世紀後半に、Borchardtという人が、種数2のテータ関数(テータ定数)の倍角公式に基づいて、算術幾何平均の一般化を作ったようだ。

[PDF] 算術平均の周辺：Borchardtによる４項算術幾何平均とThomaeの公式
https://www2.tsuda.ac.jp/suukeiken/math/suugakushi/sympo18/18_2shiga.pdf

Borchardtの論文は、オンラインでは見つけられなかった

Über das arithmetisch-geometrische Mittel aus vier Elementen
https://books.google.co.jp/books/about/%C3%9Cber_das_arithmetisch_geometrische_Mitt.html?id=iA0inQAACAAJ&redir_esc=y

通常の算術幾何平均では、楕円曲線y^2=4(x-a_1)(x-a_2)(x-a_3)(a_1 \lt a_2 \lt a_3)に対して、2つの周期
\omega_1 = \dfrac{\pi}{M(\sqrt{a_3-a_1} , \sqrt{a_3-a_2})}
\omega_2 = \dfrac{i \pi}{M(\sqrt{a_3-a_1} , \sqrt{a_2-a_1})}
を計算できるのと同様、Borchardt AGMの場合は、2つの2x2周期行列の行列式が計算できる。BorchardtのAGMは、テータ関数側から作ってるので、超幾何関数に相当するようなものは、調べた限り、明示的に書かれているものはない。時期的には、Borchardtの論文は1876年で、Appelの超幾何級数の論文が1880年、Lauricellaの超幾何級数の論文が1893年(Wikipedia情報)なので、多変数超幾何関数の方が新しい


最近は(?)、他にも、色んな"平均"が構成されていて、これらは、超幾何関数の変換公式を利用している

[PDF] 超幾何関数の関数等式と算術調和平均およびその拡張
https://www.riam.kyushu-u.ac.jp/fluid/meeting/18ME-S5/papers/Article_No_38.pdf

超幾何関数の変換公式と平均反復
http://fe.math.kobe-u.ac.jp/Movies/cm/2009-01-08-matsumoto.pdf


このへんを、何気なく、見てると、とある四項平均が、三変数のLauricella超幾何関数の二乗で書けるという式が書いてあったので、Borchardt AGMも、意外と簡単な形で書けたりするのだろうかと疑問に思った。


頑張れば、超楕円積分をLauricella超幾何関数で書いて解析的に計算できるかもしれないけど、別に重要な問題とも思えないので、そこまでモチベーションがない。Mathematicaなら、超楕円積分を適当な超幾何関数で書くくらい、やってくれないかと期待したけど、無理だったので、数値計算でやってみることにした。

何もなしには出来ないので、いくつか仮定を置く。W(a,b,c,d)をBorchardt AGMとして、
\dfrac{1}{W(1,x_1,x_2,x_3)} = F(1-x_1^2 , 1-x_2^2 , 1-x_3^2)
となる冪級数Fがあり、係数は、有理数で書けるものとする。(x1,x2,x3)=(1,1,1)の時、左辺は1なので、F(0,0,0)=1である。収束半径とかは分からないけど、どうせ原点近くでの振る舞いしか見ないので、気にしなくていいだろう

Borchardt AGMの計算は、平方根を使うので、有理数だけで済ますことはできず、倍精度程度の浮動小数点数では全然精度が足りない。ということで、任意精度演算をやるために、pythonのmpmathを使う。やることは簡単なので、別になんでもいいけど

最初は、3変数の内、2つは0とする。
F(z,0,0) = 1 + c_{100}z + c_{200} z^2 + \cdots
で、zが十分小さければ、高次の項の寄与は、殆ど無視できるはず

例えば、一次の項を決めるのに、以下のようなコードを実行する

from mpmath import *
mp.dps = 500


def M(x,y,z):
   a,b,c,d = 1,x,y,z
   max_iter = 1000
   for n in range(max_iter):
      a2 = (a+b+c+d)/4
      b2 = (sqrt(a*b) + sqrt(c*d))/2
      c2 = (sqrt(a*c) + sqrt(b*d))/2
      d2 = (sqrt(a*d) + sqrt(b*c))/2
      a,b,c,d=a2,b2,c2,d2
   return a


if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c = (1/M(x1,x2,x3)-1)/z1
    print(str(c)[:50])

結果は、

0.125001171883984440643884854175753723800930555241
0.125000000011718750000898437927005455667004923082
0.125000000000000117187500000000098949851548046299
0.125000000000000000001171874999999999935735474360

とかなるので、まぁ、0.125=1/8に収束しそう、と分かる。どこに収束するかは勘で判断してるだけので、合ってるかどうか証明とかはない

同様に、z_1^2の係数は、以下のようなコードで決める

if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20,1.0e-30]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c=(1/M(x1,x2,x3)- (1+z1/8))/(z1*z1)
    print(str(c)[:50])

結果は

0.058594492194493161915357892508282237042562173272
0.058593750007421875000699310573188476983608793132
0.058593750000000074218750000000075698227920533509
0.058593750000000000000742187499999999959300436760
0.058593750000000000000000000000074218750000000006

で、0.0589375 = 15/256に収束しそう

続ける。

if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20,1.0e-30]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c=(1/M(x1,x2,x3)- (1+z1/8+15*z1*z1/256))/(z1*z1*z1)
    print(str(c)[:50])

で、

0.037109910207646803092967250526967109948749607120
0.037109375005352020264233913616668406369705086927
0.037109375000000053520202636718810362650910452159
0.037109375000000000000535202026367187470651540126
0.037109375000000000000000000000053520202636718754

とかになるので、0.037109375=19/512だろう

だんだん辛くなってきたけど

if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20,1.0e-30]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c=(1/M(x1,x2,x3)- (1+z1/8+15*z1*z1/256+19*z1*z1*z1/512))/(z1*z1*z1*z1)
    print(str(c)[:50])

が

0.026760516142735359850255545761188337409031621442
0.026760101322507572174538118990437334607904260568
0.026760101318359416481971740722706058694628176862
0.026760101318359375000414819717407226539753154096
0.026760101318359375000000000000041481971740722659

は、0.026760101318359375=7015/262144らしい。などとやっていくと、2つの変数が0の場合は決まる。

0でない変数が複数ある場合は、面倒くさくなる。とりあえず、二次、三次の係数は、普通にやればいい。z_1^3z_2なんかは、例えば

def F3(z1,z2,z3):
   v1 = (z1+z2+z3)/8
   v2 = 15*(z1*z1+z2*z2+z3*z3)/256 + 3*(z1*z2+z2*z3+z3*z1)/128
   v3 = 19*(z1*z1*z1+z2*z2*z2+z3*z3*z3)/512 + 3*(z1*z1*z2+z2*z2*z1+z1*z1*z3+z2*z2*z3+z3*z3*z1+z3*z3*z2)/256 + 3*z1*z2*z3/512
   return 1+v1+v2+v3

if __name__=="__main__":
  mp.dps = 1000
  for eps in [1.0e-20,1.0e-40,1.0e-80,1.0e-120]:
    for ratio in [0.9,0.0011,1.1e-30,1.1e-60]:
       x1,x2,x3=1-mpf(eps),1-mpf(eps)*mpf(ratio),1
       z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
       diff = (1/M(x1,x2,x3)- F3(z1,z2,z3) - 7015*(z1**4+z2**4+z3**4)/262144)
       if diff < z1*z1*z1*z2 and diff > z1*z1*z2*z2:
          c = diff/(z1*z1*z1*z2)
          print(str(c)[:50])

とやると

0.007620766506805420299495398378483141958519102149
0.007620766506805419922275006967238197476573541568
0.007614135745958588340065696042324398959715200301
0.007620766506805419922275006967238197472801337654
0.007614135742187500000000000000006621551513671874
0.007614135742187500000377108834006569608349931553
0.007620766506805419922275006967238197472801337654
0.007614135742187500000000000000006621551513671874
0.007614135742187500000000000000000000000000000000

なので、0.0076141357421875=499/65536と求まる


そんな感じで計算してみた表が、以下。合ってるのかどうか分からないけど、分母が、2のべきになったので、多分合ってる。以下が計算値のまとめ

i j k c(i,j,k)
1 0 0 1/8
2 0 0 15/256
1 1 0 3/128
3 0 0 19/512
2 1 0 3/256
1 1 1 3/512
4 0 0 7015/262144
3 1 0 499/65536
2 2 0 789/131072
2 1 1 201/65536
5 0 0 43487/2097152
4 1 0 11685/(2^21)
3 2 0 4161/(2^20)
3 1 1 1059/(2^19)
2 2 1 15123/(2^21)
6 0 0 282387/(2^24)
5 1 0 36591/(2^23)
4 2 0 49185/(2^24)
7 0 0 946293/(2^26)
8 0 0 3323953575/(2^38)

意外と単純な係数が出たけど、超幾何関数そのものではないし、超幾何関数の二乗とかいう可能性も、なさそう。3変数のままやると辛いので、2つの変数は0にすると
F(z,0,0) = 1 + z/8 +\dfrac{15}{256} z^2 + \dfrac{19}{512} z^3 + \dfrac{7015}{262144} z^4 + \cdots
なので、形式的冪級数として、平方根を計算ことができるが、
\sqrt{F(z,0,0)} = 1 + \dfrac{1}{16} z + \dfrac{7}{16^2} z^2 + \dfrac{69}{16^3} z^3 + \dfrac{6267}{2^19} z^4 + \cdots
で、だめっぽい。

なんか、それらしいLauricella超幾何級数のいくつかの積になったりしないだろうかと思って、脳死状態で、コンピュータに計算してもらったけど、候補を見つけることはできなかった。まぁ、あんまり単純な式はないのかもしれない

とりあえず、数値計算で、冪級数の係数を予想するというのは初めてやったけど、意外と、うまくいくらしい


種数3以上の超楕円積分に対して、Borchardt AGMの一般化が、以下で与えられている。
Hyperelliptic integrals and generalizedarithmetic–geometric mean
https://link.springer.com/article/10.1007/s11139-011-9353-7

防御力の物理〜cavity expansion model

弓矢の貫通力
https://m-a-o.hatenablog.com/entry/2019/09/22/185608

に書いたcavity expansion modelは、由来を、ちゃんと調べなかったけど、起源は、1945年の論文

[1] The theory of indentation and hardness tests
https://doi.org/10.1088/0959-5309/57/3/301

にあるらしい。オリジナルのcavity expansion modelは、延性破壊を起こす材料を対象とした、比較的単純なモデルだけど、その後、色々な魔改造の末、土壌や脆性材料への適用や、時間依存性を考慮した計算(dynamic cavity expansion model)なども現れたという経緯っぽい

著者は、R.F.Bishop , R.Hill , N.F.Mottとなっていて、Mottは、Mott散乱公式やMott転移、Mott絶縁体に名前が残っている物理学者Mottと同一人物っぽい。Mottは、自伝を書いていて(絶版になってるっぽいけど、日本語訳も出ている)、戦時中は、軍事研究に携わったことが書かれているけど、この論文について、明示的な言及は別にない

cf)科学に生きる―ネビル・モット自伝
https://www.amazon.co.jp/gp/product/4532062748

R.Hillは、おそらくRodney Hillで、塑性力学の分野では、そこそこ知られた名前のようである。Rodney Hillについては、以下に詳しく書かれている

Rodney Hill. 11 June 1921 — 2 February 2011
https://royalsocietypublishing.org/doi/abs/10.1098/rsbm.2014.0024

"In 1943 he moved to the Armament Research Department at Fort Halstead in Kent, for three years. Here he was involved in, for example, the modelling of armour penetration by projectiles."だそうな。cavity expansion modelも、装甲の弾丸による貫通を調べる過程で生まれたのかもしれない。

R.Hillは、1950年に以下の本を出版している。

The Mathematical Theory Of Plasticity
https://www.amazon.co.jp/dp/0198503679

日本語訳もあったが、絶版らしい

塑性学 (1954年)
https://www.amazon.co.jp/dp/B000JB7J1I

R.F.Bishopという人のことは、よく分からない。


論文[1]のタイトルにある"indentation and hardness tests"は、押し込み硬さ試験(Indentation hardness tests)のことではあるんだろうけど、通常の押し込み硬さ試験より荷重を大きくしていって、押し込み深さが非常に深くなるような、穿孔に近い状況を調べているっぽい。

いわゆる押し込み硬さは、Brinell硬度とかVickers硬度とかで、くぼみができるまで圧子を押し付けて測定する。圧子の大きさと比べて、そんなに深くは押し込まないと思う。押し込み硬さの歴史はそんなに古くなくて、1900年に、スウェーデンのBrinellが導入したBrinell硬さが、現在使われてる中で、一番古い基準らしい。何故か有名なモース硬度とかは、もっと古くからあるけど、あんまり定量的ではない。


論文[1]で扱われてる対象は、穿孔機(先端形状は、円錐で、先端の角度は複数ある)を、非常に深く押し込んだ時の(単位面積当たりの)抗力で、論文[1]のFigure3とFigure4で与えられた実測結果から、深さが十分であれば、それはある一定値に収束すると見なせそうではある。この抗力の大きさは、論文で使っている材質(hardened copperとannealed copperと書いてある)では、降伏応力の4〜5倍に達する。

Figure3とFigure4の結果では、極限抗力は、多少幅があるけど、おおまかには、100tons/sq.inch.程度っぽい。単位がわかりづらいけど、100[tons/sq.inch]=1.0e3*100*9.8/(25.4*25.4)=1520[N/mm^2] = 1520[MJ/m^3]程度。これは、銅なので、前回の銃弾によるスチールの貫通抗力よりは、小さくなっている。

前回と同じく、比較として、単位体積の銅を融解するのに必要な熱エネルギーを見ておくと、比熱(※)は、常温以上の金属では、大体、3R≒25(J/mol/K)で、融点が1000度強らしいので、融点まで熱するのに、ざっくり25(kJ/mol)で、融解熱が13.26(kJ/mol)とある。モル質量が63.55(g/mol)で、密度が8960(kg/m^3)から、8960*(25+13.26)*1.0e-3/63.55e-3≒5400(MJ/m^3)程度だろう。1520/5400=0.28…で、割合としては、球状弾でスチールを貫いた時より、多少小さいけど、大きな差はない

※)面倒なので、Dulong-Petit則が使える定積比熱を用いてるけど、定圧比熱を使うほうが適切な気はする。常温の固体では、定積比熱と定圧比熱は大差なくても、高温では、それなりの差が出る可能性がある。けど、オーダーを見るだけなので気にしない

穿孔の手段として、現代だとレーザーなんかもあると思うけど、エネルギー的には、レーザーより、ずっと効率がよさそう。


この抗力の大きさを見積もるための論文[1]の方針は、以下のように書かれている。

Denoting ths maximum load by p_0 A, where A is the area of the cross-section of the punch, it is found that p_0, for a lubricated punch is about twice the hardness, or five times the yield stress, of the work-hardened material. A theoretical method is given for calculating p_0, as follows : the pressures p_c, and p_s, required to enlarge a cylindrical and a spherical hole in a material showing any kind of strain hardening can be calculated. It is plausible to assume that p_0, should be between p_c, and p_s, and since p_s, is only slightly greater than p_c, an approximate theoretical estimate of p_0, is obtained. This is in good agreement with experiment.

抗力の大きさは、hardnessの二倍とか書いてある。このhardnessが何を指してるのか明言されてないけど、何らかの押し込み硬さだろうと思う

cavity expansion modelで計算する力は、物体内に作られた円筒状あるいは球状の穴の径を拡大する時に受ける力であって、穿孔機がプレートを貫く時の抗力とかとは(少なくとも、見かけ上は)違うものである。plausibleとか言ってるので、Bishop-Hill-Mottの主張も、それほど強い根拠があってのものでもなさそう

なので、安易に同一視しないため、論文では使われてない用語であるけど、p_cをcylindrical cavity expansion pressure、p_sをspherical cavity expansion pressureと呼んでおく。

また、摩擦については本文内に言及はあるものの見積もれないと書いてある。


論文[1]の計算では、線形硬化弾塑性体(降伏応力以下では、線形弾性体で、降伏応力以上でも、ひずみと応力に一次の関係式が成立する。降伏応力以上で、ひずみの応力に対する比例係数を塑性係数と呼ぶ)を仮定した計算を実行しているけど、他のモデルを利用することもできる。ただ、塑性係数を0にしても比較的影響は小さいようなので、真面目に扱ってみても、大した意味はないかもしれない。

計算の主要な部分、特に、cylindrical cavityに関する部分は、1931年出版の

Plasticity: A Mechanics Of The Plastic State Of Matter
https://archive.org/details/in.ernet.dli.2015.166225/page/n213

という本に既に書かれている。この本は、塑性力学に関する、(少なくとも、英語で書かれたものとしては)世界最初の教科書らしい。The thick-walled tubed under internal pressureというタイトルの章がある。論文[1]独自の部分には、妥当性の疑わしい議論が含まれると思う。

論文は、古い時代のもあって、よく分からない計算が含まれているので、現代の適当な教科書を見たほうがいい。塑性力学の教科書は沢山はないけれど、日本語で読めるものから探すと、必要な計算の(全てではないが)殆どは、例えば、

基礎塑性力学 (実用理工学入門講座)
https://www.amazon.co.jp/dp/481730149X
http://www.nissinpb.co.jp/jitsuyorikogaku-kikai.htm#%E5%9F%BA%E7%A4%8E%E5%A1%91%E6%80%A7%E5%8A%9B%E5%AD%A6

の6.4節、6.5節とかに書いてある。


spherical cavityの場合、球の中心からの距離をsとして、cavity半径がaであるとし、 a \leq s \leq cを塑性領域、c \leq sを弾性領域とする。

Poisson比を\mu < 0.5として、弾性領域c \leq sでは、問題なく、consistentな解
u(s) = \dfrac{(1+\mu)Y}{3E}\left( \dfrac{c^3}{s^2} \right)
\sigma_r(s) = -\dfrac{2Y}{3}\left( \dfrac{c}{s} \right)^3
\sigma_{\theta}(s) = \dfrac{Y}{3}\left( \dfrac{c}{s} \right)^3
\epsilon_{r}(s) = -\dfrac{2(1+\mu)Y}{3E}\left( \dfrac{c}{s} \right)^3
\epsilon_{\theta}(s) = \dfrac{(1+\mu)Y}{3E}\left( \dfrac{c}{s} \right)^3
が見つかる。uはradial displacementで、Eはヤング率、Yが降伏応力etc.。\sigma_{r},\sigma_{\theta},\epsilon_r,\epsilon_{\theta}は、radial stress, tangentical stress, radial strain, tangentical strain。

現在の大抵の教科書は、球全体の半径bを有限としているが、式が複雑になるので、単に、b→∞の極限を取った。弾性領域の計算は、古典的な弾性論通りなので、特に問題はない

これらの関数が満たすべき条件は、いくつかある(Hookeの法則とか平衡方程式とか)けど、全部書くと面倒なので省略。降伏条件は、Misesの降伏条件でもTrescaの降伏条件でも同じで
\sigma_{\theta}(c) - \sigma_{r}(c) = Y
を採用している。

弾性領域と塑性領域の境界に於いて
u(c) = \dfrac{(1+\mu)Y}{3E}c
であるけど、金属などでは、ヤング率は降伏応力よりずっと大きいので、u(c)/cは小さい。

塑性領域の体積が空洞形成前後で不変だと仮定すると
\dfrac{4\pi}{3}(c^3 - a^3) = \dfrac{4\pi}{3} (c-u(c))^3
で、近似的に
u(c) \approx \dfrac{a^3}{3c^2}
となって、
\left( \dfrac{c}{a} \right)^3 = \dfrac{E}{(1+\mu)Y}
が導ける。

塑性領域a \leq s \leq cで平衡方程式
\dfrac{d \sigma_r}{ds} = \dfrac{2(\sigma_{\theta} - \sigma_r)}{s}
が成立して、塑性係数が0なら、塑性領域に於いて
\sigma_{\theta}(s) - \sigma_{r}(s) = Y
となって、平衡方程式を積分して、応力が弾塑性領域境界で連続だと仮定すると
\sigma_r(s) = \sigma_{r}(c) - 2Y \log(\dfrac{c}{s}) = -\dfrac{2Y}{3} - 2Y \log(\dfrac{c}{s})
と計算できる。

知りたいcavity expansion pressureは、cavityと塑性領域の境界位置s=aに於ける-\sigma_{r}(a)の値なので、塑性領域に於ける歪みの詳細が分からなくても、欲しい答え
P = \dfrac{2Y}{3} + 2Y \log(\dfrac{c}{a})
が得られる(\dfrac{c}{a}は既に決定した)。

論文[1]では、塑性係数が0とは仮定していないので、かなり乱暴に、塑性領域での歪みを決め打ちしているけど、弾塑性領域境界での連続性も成立してないし、空洞・塑性領域境界では発散しているので、あまり鵜呑みにできない。しかしまぁ、計算された結果を信じるなら、
P = \dfrac{2Y}{3} + 2Y \log(\dfrac{c}{a}) +\dfrac{2\pi^2}{27}A
となる。但し、Aは塑性係数。



前回計算した"破断エネルギー(密度)"と比較されるべき量は、(cavity expansion pressureではなく)弾性変形領域に蓄えられている歪みエネルギーと塑性仕事の和を、cavityの体積で割った量になるはず。cavity形成に投入されたエネルギー(銃弾による貫通の場合は、着弾時運動エネルギー)は、弾性域のひずみエネルギー+塑性仕事+その他散逸に等しい。その他散逸というのは、表面摩擦熱などだけど、今は無視する

塑性仕事の計算には、塑性歪み増分を知る必要があるけれど、塑性係数が0の時は、適当な仮定を置くことで、塑性域の歪み(とradial displacement)をモデル化する方法がある(論文[1]には書いてないけど、教科書"基礎塑性力学"には書いてある)。

塑性領域では
\epsilon_r = \dfrac{1}{E} \left(\sigma_r - 2\nu \sigma_{\theta}\right) - \dfrac{2}{3}(\sigma_{\theta} - \sigma_{r})\phi
\epsilon_{\theta} = \dfrac{1}{E} \left( (1-\nu)\sigma_{\theta} -\nu\sigma_r \right) + \dfrac{1}{3}(\sigma_{\theta} - \sigma_{r})\phi
の形の歪みを仮定しても矛盾は、ない。\phiは関数で、それぞれ、歪みの第一項は、Hookeの法則と同じ形なので、弾性歪み成分と呼ばれ、第二項は塑性歪み成分と呼ばれる。

特に、この形の歪みは、
\epsilon_r + \epsilon_{\theta} = \dfrac{1-2\nu}{E} (\sigma_r + 2\sigma_{\theta})
を満たす。radial displacementは
\epsilon_r(s) = \dfrac{du}{ds} , \epsilon_{\theta}(s) = \dfrac{u}{s}
を満たすように決める。既に、応力は計算したので、これらの条件から、歪みを決められる。具体的には
\phi(s) = \dfrac{3(1-\nu)}{E} \left( \dfrac{c^3}{s^3} - 1 \right)
となり、歪みの塑性成分は、
\epsilon_r^{(p)}(s) = \dfrac{-2Y}{3} \phi(s)
\epsilon_{\theta}^{(p)}(s) = \dfrac{Y}{3} \phi(s)
のように決まる。

計算したい量は、半径aのcavityを準静的に形成する塑性仕事W(a)で、塑性仕事増分は、
\dfrac{dW}{da} = \displaystyle \int_{a}^{c} 4\pi s^2 \left( \sigma_r(s) \dfrac{d\epsilon_{r}^{(p)}}{da} + 2\sigma_{\theta}(s)\dfrac{d\epsilon_{\theta}^{(p)}}{da} \right) ds
と求まるだろう。
c^3 = \dfrac{E}{(1+\nu)Y} a^3
に注意すると、(当然、W(0)=0なので)私の計算では
W(a) = \dfrac{4 \pi a^3}{3} \dfrac{6(1-\nu)Y}{1+\nu} \log(\dfrac{c}{a})
となった。 ポアソン比が、0.5なら、
W(a) = \dfrac{4 \pi a^3}{3} (2 Y \log(\dfrac{c}{a}))
となって、弾性変形による歪みエネルギー(この計算は、弾性論の範疇なので省略)をcavity体積で割った量が、丁度、2Y/3なので、足すとspherical cavity expansion pressureに等しくなる。

金属とかだと、ポアソン比は、0.3〜0.4くらいだったりするので、塑性仕事+弾性歪みエネルギー/空洞体積の値は、spherical cavity expansion pressureより、もう少し大きい値になる。文献では見当たらない結果なので、見なかったことにしよう(まぁオーダーは合ってるし...

実際のとこ、c/aは定数としているけど、弾丸がプレートを貫通する前とかは、a=c=0で、これらの仮定の妥当性が疑われる。また、c/aが上で決めた値と違っても、定数であれば、W(a)は空洞体積に比例するので、c/aが素早く定数値に近付くのであれば、結果は少し変わるけど、計算自体に変更は要らない。例えば
\left( \dfrac{c}{a} \right)^3 = \dfrac{E}{3(1-\nu)Y}
であれば、
W(a) = \dfrac{4 \pi a^3}{3} (2 Y \log(\dfrac{c}{a}))
が(ポアソン比に依らず)常に成り立つ。



加えられた塑性仕事の大部分は、熱として散逸する(らしい)。もし、空洞形成が(熱伝導より)素早く起きるなら、塑性仕事を塑性域の体積で割って、更に比熱で割ると、塑性域の温度が何度くらい上昇するか見積もれそうに思える。金属の場合、塑性域の体積は空洞体積よりずっと大きいので、簡単のため、塑性域体積+空洞体積で割る。どうせ概算なので、ポアソン比は0.5としておく。すると
W(a)/(\dfrac{4 \pi c^3}{3}) = \dfrac{3Y^2}{E} \log(\dfrac{c}{a})
で、E=200(GPa)とY=500(MPa)を採用すると、19(MJ/m^3)くらいの塑性発熱があるだろう。鉄の場合は、比熱が3.5[MJ/m^3/K]程度なので、ざっくり5度くらいの温度上昇がある計算になる。

ヤング率も降伏応力も、温度によって変化する(なので、鉄を溶かせなかった時代でも、熱した鉄を叩いて加工できた)が、5度くらいなら影響を無視できるだろう

現在では、高速な塑性変形で材料加工する際の温度変化を計測して、塑性変形量を計測する技術というものが存在するっぽい

高速変形試験（4）～高速変形時の加工発熱分布の計測技術～
https://www.jfe-tec.co.jp/jfetec-news/27/5p.html

cylindrical cavityでもspherical cavityでも、cavity周りに塑性領域が存在すると考える点では同じだが、c/a(塑性領域半径/空洞半径)の値は、ヤング率、降伏応力、ポアソン比が同じ場合でも結構違う。例えば、上のE=200(GPa)、Y=500(MPa)、ν=0.5だと、cylindrical cavityでは、c/a≒15.2で、spherical cavityでは、c/a≒6.4なので、弾丸が装甲を貫いた時に、穴の周囲に、どれくらいの大きさの塑性領域が形成されるか分かれば、どちらのcavity modelに近いか判定する基準になるかもしれない



論文[1]では、(cold-workedとあるので、冷間加工した)銅に対して、
降伏応力Y=285(MPa)
塑性係数A=100(MPa)
Young率E=130(GPa)
Poisson比ν=0.34
という値を使って(応力などの単位は、tons/sq.inchなので、換算した)いて、spherical cavity expansion pressureはP=1375(N/mm^2)程度となる。穿孔機での実測値と比較するとオーダーは合ってる

論文[1]のFigure7で、軟鋼に対する応力歪み曲線が掲載されていて、鋼のヤング率は、大体200GPaくらい、ポアソン比は0,3くらいで、それ以外の値は、私が目視で読み取った限り、
降伏応力Y=275(MPa)
塑性係数A=350(MPa)
Young率E=200(GPa)
Poisson比ν=0.3
で、spherical cavity expansion pressure=1600(MPa)=1600(N/mm^2)となる。論文著者たちの計算だと、120[tons/sq.inch] ≒1960[N/mm^2]程度。計算に使った値が書かれてないけど、まぁ、こんなものだろう



前回のマスケット銃に対するGrazテストでは、ブリネル硬さ290のcold-worked mild steelを使ったとか書いてあった

出典:Material Culture and Military History: Test-Firing Early Modern Small Arms
https://journals.lib.unb.ca/index.php/MCR/article/view/17669/22312

引張強度(近似値)換算で、965(MPa)程度と思われるが、降伏応力の値は分からない。引張強度が400(N/mm^2)のSS400という鋼材だと、降伏応力が235~245(N/mm^2)と規格で決まってるので、ざっくり、引張強度の60%程度である。965(MPa)の60%で、Y=580(MPa)として、A,E,νは上記の値を流用すると、spherical cavity expansion pressure P=2800(MPa)と計算され、Grazテストの破断エネルギーとオーダーは合ってる(勿論、偶然かもしれない)

同じパラメータを使って、cylindrical cavity expansion pressureだと、P=2350(MPa)程度。

表を再掲。前回書き忘れた(元論文には書くまでもないのか何も言及されてない)けど、銃弾の素材は、(密度を計算すると)鉛っぽい。

銃器 試験 運動エネルギー "破断エネルギー" 着弾速度(m/s)
"Doppelhaken" G 284 PS100 1779.6J 3138(MJ/m^3) 305
"Doppelhaken" G 358 PS100 2992.7J 2335(MJ/m^3) 349
Matchlock LG 1514 PS30 1241.7J 3467(MJ/m^3) 378
Matchlock LG 1514 PS100 605.7J 3382(MJ/m^3) 264
Wheellock RG 33 PS30 2333.2J 3347(MJ/m^3) 394
Wheellock RG 33 PS100 1238.0J 2664(MJ/m^3) 287
Wheellock RG 117 PS30 660.2J 2779(MJ/m^3) 349
Wheellock RG 117 PS100 307.0J 2539(MJ/m^3) 238
Wheellock RG 272 PS30 1874.9J 3898(MJ/m^3) 342
Wheellock pistol RP 2895 PS30 602.4J 2754(MJ/m^3) 355
Flintlock STG 1287 PS30 2269.8J 2560(MJ/m^3) 406
Flintlock STG 1287 PS100 1166.1J 2630(MJ/m^3) 291
Flintlock STG 1288 PS30 2032.8J 3131(MJ/m^3) 390
Flintlock STG 1288 PS100 1055.3J 2438(MJ/m^3) 281
Flintlock STG 1316 PS30 2458.3J 3368(MJ/m^3) 391
Flintlock STG 1316 PS100 1324.5J 2722(MJ/m^3) 287
Flintlock STG 1318 PS30 2806.5J 2917(MJ/m^3) 426
Flintlock STG 1318 PS100 1438.6J 2991(MJ/m^3) 305
Flintlock pistol STP 1128 PS30 753.8J 2633(MJ/m^3) 323
アサルト・ライフル58 PS100 2801.5J 4824(MJ/m^3) 770
アサルト・ライフル77 PS100 1375.0J 5987(MJ/m^3) 874
Glock80ピストル PS30 467.9J 3444(MJ/m^3) 342

Grazテストでは、スプルース(木材の一つ)をターゲットにしたデータもあるので、同様の計算をしてみる

銃器 試験 運動エネルギー "破断エネルギー" 着弾速度(m/s)
"Doppelhaken" G 284 PW100 1779.6J 41.0(MJ/m^3) 305
"Doppelhaken" G 358 PW100 2992.7J 49.4(MJ/m^3) 349
Matchlock LG 1514 PW30 1241.7J 53.0(MJ/m^3) 378
Matchlock LG 1514 PW100 605.7J 40.6(MJ/m^3) 264
Wheellock RG 33 PW30 2333.2J 52.9(MJ/m^3) 394
Wheellock RG 33 PW100 1238.0J 66.6(MJ/m^3) 287
Wheellock RG 117 PW30 660.2J 42.1(MJ/m^3) 349
Wheellock RG 117 PW100 307.0J 30.8(MJ/m^3) 238
Wheellock RG 272 PW30 1874.9J 46.4(MJ/m^3) 342
Wheellock RG 272 PW100 1083.6J 43.7(MJ/m^3) 260
Wheellock pistol RP 2895 PW30 602.4J 45.5(MJ/m^3) 355
Flintlock STG 1287 PW30 2269.8J 89.0(MJ/m^3) 406
Flintlock STG 1287 PW100 1166.1J 63.4(MJ/m^3) 291
Flintlock STG 1288 PW100 1055.3J 61.0(MJ/m^3) 281
Flintlock STG 1316 PW30 2458.3J 51.8(MJ/m^3) 391
Flintlock STG 1316 PW100 1324.5J 37.0(MJ/m^3) 287
Flintlock STG 1318 PW30 2806.5J 63.8(MJ/m^3) 426
Flintlock STG 1318 PW100 1438.6J 52.5(MJ/m^3) 305
Flintlock pistol STP 1128 PW30 753.8J 46.2(MJ/m^3) 323
アサルト・ライフル58 PW100 2801.5J 120.0(MJ/m^3) 770
アサルト・ライフル77 PW100 1375.0J 187.8(MJ/m^3) 874
Glock80ピストル PW30 467.9J 54.7(MJ/m^3) 342

材料力学的パラメータは
Wood, Panel and Structural Timber Products - Mechanical Properties
https://www.engineeringtoolbox.com/timber-mechanical-properties-d_1789.html
に、
ヤング率E〜10(GPa)
降伏応力Y〜70(MPa)
とあるけど、降伏応力が破断エネルギーより大きい。そんなわけ無いと思うのだけど、木材の強度は、方向依存性があって、繊維に直交する方の降伏応力は、もっとずっと弱いようだ。cavity expansion modelは等方性を仮定したモデルなので、そのままでは適用できない


論文[1]の実験だと、抗力がサチるには、(穿孔機の直径と比較して)ある程度の深さまで、押しこむ必要があるっぽい。一方、Grazテストでは、直径1cm程度の弾丸で、数mmのプレートを貫通するという状況がメインなので、最大抗力が発揮されるのか謎である

論文[1]より少し前に、BetheやG.I.Taylor、プレートの厚みが比較的薄い状況を念頭に置いてモデル化を試みている。

まず、Betheは1941年、

Attempt of a Theory of Armor Penetration
https://books.google.co.jp/books/about/Attempt_of_a_Theory_of_Armor_Penetration.html?id=AYzkwgEACAAJ&redir_esc=y

という報告を書いてるけど、オンラインでは読めないっぽい。これに対して、G.I.Taylorは、Notes on H.A. Bethe's "Theory of armor penetration"という文書を書いていて、ネット上で読むことができる。G.I.Taylorの結果は、Betheのものとは少し違うようだけど、極端な差はないらしい

Notes on H.A. Bethe's "Theory of armor penetration", I. Static penetration,
https://apps.dtic.mil/dtic/tr/fulltext/u2/1009963.pdf

Notes on H.A. Bethe's "Theory of armor penetration", II. Enlargement of a hole in a flat sheet at high speed
https://apps.dtic.mil/dtic/tr/fulltext/u2/1009965.pdf

この2つの報告より、戦後に書かれた以下の論文を読む方がいいかもしれない。

The formation and enlargement of a circular hole in a thin plastic shee
https://doi.org/10.1093/qjmam/1.1.103

cavity expansion modelでは、空洞の周囲に、塑性変形領域と弾性変形領域の二層ができると考えていた。G.I.Taylorは、プレートの厚みが薄いことにより、空洞の周囲では、厚みが増大するとし、厚みの増大が無視できない塑性変形領域、厚みの増大が無視できる塑性変形領域、弾性変形領域の三層で考えている(Betheも、同様の仮定を採用したようだ)。

一番内側の厚みの増大が無視できない塑性変形領域では、軸中心からの距離sに於ける厚さh(s)を、新しい変数として導入している。プレートは薄く、完全に貫通している状況を考えているので、計算するべきは、抗力ではなく、穴を形成するのに必要なエネルギーということになる。

G.I.Taylorの議論の詳細は見てないけど、結論は、"the work done in expanding to a given radius"が式(81)に書いてあるけど、単位体積あたりの破断エネルギー/降伏応力は、1.33になるらしい。これだと、Grazテストの結果とは合わなそうに思える。モデル化自体に問題はない気がするので、塑性仕事の計算を間違えてるんじゃないかと思うけど、面倒くさいので、詳細は確認してない



ちなみに、G.I.Taylorの1946年の講義録(?)

THE TESTING OF MATERIALS AT HIGH RATES OF LOADING
https://www.icevirtuallibrary.com/doi/pdf/10.1680/ijoti.1946.13699

には、以下のように書いてある。

Lead bullets weighing from 1.8 gram to 2.7 grams were fired from B 0.22-inch rifle at the mid-point of the anvil. We found that the lead behaved like a fluid, spreading out in a thin film over the plane surface of the anvil. me measured the force exerted by the bullet and found it to be, in fact, just what a fluid jet of diameter 0.22 inches and the density of lead would exert.

I have already mentioned that a lead bullet fired at a speed of about 1,OOO feet per second produces a thin spray of lead on striking a steel target, just as a jet of water would. To understand in a quantitative way what happens when impact occurs at such speeds would, in general, be beyond the power of our present analytical methods, even if we knew beforehand the stress-strain characteristics of the material.



もうひとつ、上の話とは直接関係ないけど、800m/s〜2500m/sの高速度領域ではあるが、タングステン・カーバイド製の非常に硬い球状の弾丸による貫通深度の実測値が以下の論文に記載されているのを見つけた。

[2]Penetration of HSLA-100 steel with tungsten carbide spheres at striking velocities between 0.8 and 2.5 km/s
https://www.sciencedirect.com/science/article/pii/S0734743X03000800

ターゲットのHSLAプレートは、Table3によれば、
ヤング率E=197(GPa)
降伏応力Y=730(MPa)
Poisson比ν=0.29
らしい。流石に塑性係数は記載されてないので、0としておくと、spherical cavity expansion pressureはP〜3100(MPa)である。

この論文で使用している弾丸は、直径6.4mmの球状弾丸で、タングステン・カーバイドの密度は、14.9(g/cm^3)とTABLE3に記載されているので、弾丸質量は、2.045(g)程度

実際の実測値は、Table2を見ると、まず、速度が1000m/s程度までは、形成されたクレーター直径は、弾丸直径と殆ど変わらないようであるが、それ以上では、クレーター直径の増大が起きるらしい。高速になるにつれて、貫通深度は伸び悩むのに、クレーター直径は、ますます大きくなる。

データに含まれるクレーター体積のデータを使って、衝突時の運動エネルギー/クレーター体積によって"破断エネルギー密度"を計算してみる。データは、12個の異なる速度のものがあるけど、そのうち、8個を気分で選んで計算した結果が以下

Shot 速度(m/s)　 体積(ml) 破断エネルギー 深さ(mm)
1 830 0.11 6404(MJ/m^3) 4.57
3 980 0.16 6138(MJ/m^3) 5.46
4 1270 0.27 6872(MJ/m^3) 7.09
6 1500 0.35 6574(MJ/m^3) 8.51
8 1910 0.5 7461(MJ/m^3) 8.53
9 2150 0.73 6475(MJ/m^3) 8.41
11 2460 0.93 6654(MJ/m^2) 9.50
12 2550 1.0 6649(MJ/m^3) 9.60

破断エネルギー密度のばらつきは、Grazテストの計算結果よりも小さく、ほぼ一定っぽい。spherical cavity expansion pressureと比較すると、まぁ、二倍くらい？

これくらいずれるとなると、Grazテストの時とは、大きく異なることが起きてる可能性があるが、弾丸材質が違うせいということもありえる(プレートの材質も微妙に違うけど、こっちの差は小さそうに思える)

[3]Impact Phenomena at High Speeds
https://aip.scitation.org/doi/10.1063/1.1722215

[4]Crater Formation in Metallic Targets
https://aip.scitation.org/doi/abs/10.1063/1.1723437

などの古い論文を読むと、様々な材質の球を1000m/sから5000m/s(もしくは750m/sから2250m/s)まで加速して、同じターゲットに衝突させると、形成されたクレーター体積と弾丸の運動エネルギーの比は、弾丸の材質が同じなら一定値を取るということが書いてある。どちらの論文も、単位が、m^3/jouleになっているが、逆数を取ると

弾丸材質 [4]TABLE I [3]TABLE II
亜鉛 1700MJ/m^3 5000MJ/m^3
錫 960MJ/m^3
銅 1700MJ/m^3
鉄 3600MJ/m^3
Steel 9000MJ/m^3
鉛 200MJ/m^3 610MJ/m^3
Brass 1000MJ/m^3
マグネシウム 1390MJ/m^3
銀 1250MJ/m^3
アルミニウム 1300MJ/m^3
Wax

という感じ。かぶってる材質が殆どないが、亜鉛、鉛では、比率が大体、1:3になっていて、この違いは、プレート側の材質に依存するものだろう。鉄とSteelは一応、区別した

鉛と鉄、あるいは、鉛とsteelでは、10倍以上違いが出てるので、論文[2]の結果とGrazテストの結果の差は、弾丸速度の差によるものかなぁという気がする


高速衝突で見られる破壊モードは、いくつかあるとか書いてあったりする(プラグ破壊とかスポール破壊とか...)けど、意味のある分類なのか、よく分からない

cf)スポール破壊条件の研究
https://www.jstage.jst.go.jp/article/jsms1963/41/468/41_468_1396/_article/-char/ja/

The interrelation of failure modes observed in the penetration of metallic targets
https://www.sciencedirect.com/science/article/pii/0734743X84900010

以下には、アルミニウム合金プレートに、900m/sから1500m/sの弾丸を打ち込んだ時、衝突速度の増大に伴って、破壊モードが、プラグ破壊からdiscingに変化したとか書いてある

Thick AA7020-T651 plates under ballistic impact of fragment-simulating projectiles
https://doi.org/10.1016/j.ijimpeng.2015.08.001

(論文[2][3]を信じるなら)ある程度の高速度領域に於いて、装甲と弾丸材質が同じなら、着弾速度に依存しない破断エネルギー密度が出るのは、不思議なことではあるが、破壊モードの違いにも関わらず、破断エネルギー密度が一定なのだろうかとか、色々疑問は沸くけど、細かいことは、よく分からない。

これらの値を(解析的に)計算できるモデルは、ないようである。また、運動エネルギーが同じなら、柔らかい弾丸の方が、大きなクレーターを形成できるようだけど、これもなんだか不思議な気がする

低次視覚野の計算モデル

今まで、(低次)視覚神経系のモデルについて、見聞きした雑多な話題として、

(1)カブトガニの眼で側方抑制がどうたら...(1950年代後半〜1960年代初頭、Hartline & Ratliff)
(2)ネコの一次視覚野(V1野)に、特定方向の線分に反応する細胞があるとか何とか(1959,Hubel & Wiesel)→1981年のノーベル医学・生理学賞
(3)V1野の方位選択性の機能は、Gaborフィルターでモデル化できる(Jones & Palmer 1987)
(4)infomax原理がどうとか

くらいがあったけど、それぞれの話の繋がりがどうなってるのか、よく分かってなくて、理解しようと思った。


(1)は、生理学とか生物学の結果で、画像処理的には、エッジ強調/画像鮮鋭化に対応する処理が行われているとされる。論文としては
Inhibitory interaction of receptor units in the eye of Limulus.
https://www.ncbi.nlm.nih.gov/pubmed/13398569

など。元々の話は、カブトガニを使ったものだったけど、マッハバンドとかハーマン格子のような古くから知られる錯視を説明でき、人間の網膜にも存在するよう。いずれにせよ、基本的には、網膜で処理されるレベルの話のよう。

マッハバンドは、Wikipediaには「微妙に濃淡の異なるグレーの領域が接触している場合に、暗い方の領域の境界付近はより暗く、明るい方の領域の境界付近はより明るく強調されて見える、錯視の一種である」と説明されている。実際の処理は、一種のhigh-passフィルタでモデル化できると考えられ、このようなhigh-passフィルタは、比較的単純な神経回路で作ることができるということが、基本的には(1)で解明されたらしい。



網膜にある視細胞で受けた情報は、次に、LGN(外側膝状体)という領域に送られる。視細胞には、rod cellとcone cellの二種類があるけど、基本的には、視細胞→双極細胞→網膜神経節細胞(RGC)→(視神経)→LGN→一次視覚野(V1)と伝わっていくらしい。rod cellは、明暗にのみ反応し、cone cellは色情報を担当するとされ、前者は、片目で一億個くらいあるのに対して、後者は500〜700万個くらいらしい。

cone cellも一様に分布してるわけじゃなく、外縁部は低密度らしいけど、平均的には10000x10000くらいの解像度があるということで、8Kモニターより上。cone cellは、大多数の人類が、三色認識するので、一色あたり200万くらいと考えると、色情報についてはフルHD程度の解像度(?)。

視細胞〜網膜神経節細胞には、水平細胞とアマクリン細胞というのもいて、沢山種類があるけど、機能の全容が解明されているわけではないらしい

視覚系の神経細胞は、それぞれが、視野の中の一部の領域からの入力にのみ反応するようになっていて、この領域を受容野と呼んでいる。画像処理で、Gaussian blurとかのフィルター処理をかける際には、サイズが3x3とか5x5の領域の情報から、一ピクセルの色情報を計算するのと同じような話。一般に、受容野は、低次視覚野の神経細胞では狭く、高次視覚野では、広くなっていくらしい。

RGC細胞数が約100万個とか書かれてたりするので、rod cellの数との比は1:100くらい。単純に考えると、10x10くらいのrod cellからの入力が一つにまとめられている計算になる(?)。この比率は、網膜位置によって、かなり異なるらしいので、目安程度にしかならない。

低次視覚野は、局所的な幾何学的特徴に反応するのに対して、高次視覚野では、顔や、特定の物体に反応する細胞とかがあるとされている。


で、神経細胞は、受容野に入ってきた入力をまとめて、ある大きさの応答を起こす。RGCやLGNの応答は、DoGフィルタとしてモデル化でき、一次視覚野の単純細胞の応答は、Gaborフィルタでモデル化できるとされている。というのが、(2)や(3)の話。DoGフィルタによるモデル化は、誰が言い出したのか分からないけど、よく書いてある

と言われても、フ〜ンという感想しかなかったので、実際に、DoGフィルタとかGaborフィルタを通した結果が、どんな感じになるか見てみた。

DoGフィルタやGaborフィルタには、いくつかのパラメータが存在し、動物や人間の神経系のモデルとして、どのようなパラメータを選ぶべきか、通常、何の説明もない。
Analysis of multidimensional difference-of-Gaussians filters in terms of directly observable parameters
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3789628/

Modeling lateral geniculate nucleus response with contrast gain control. Part 1: Formulation
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3918962/

などでは、DoGカーネルを3つのパラメータで規定される
f(x,y) = \dfrac{1}{2 \pi \sigma_C^2} \exp(-\dfrac{x^2+y^2}{2\sigma_C^2}) - \dfrac{\beta}{2 \pi \sigma_S^2} \exp(-\dfrac{x^2+y^2}{2\sigma_S^2})
の形に取っている。βをbalance parameterと呼び、無次元量である。他の2つのパラメータは、長さの次元を持つので、面倒だけど、例えば、2つの比は無次元量である。積分すると、1-βになる。

CV業界では、β=1として、残り2つのパラメータ比は、1.6にすることが多いっぽい。上の論文は、観測可能な量から各パラメータを見積もろうというノリらしいけど、実際の値は与えられていない。後者の論文の§2に於ける記述によれば、LGNにおけるβの値は、生物種や細胞種ごとに、0.95から0.52まで、色んな見積もりがあるらしい。

DoGフィルタは、フィルタを通しただけでは、あまり意味はなくて、ゼロ交差の検出をしなければならない。このへんは、工夫の余地があるのかもしれないけど、何も考えない実装になっている。

あと、適当に検索すると、gray画像のピクセル値をuint8のまま、2つのGaussianフィルタに通して(cv2.filter2Dを使うと、適当に近い整数値を選んで、出力もuint8のままにしてくれる)、そのまま引くという実装をしている人がいた。大体の場合は、引き算した値は、0に近い値になるけど、uint8同士の引き算なので、オーバーフローを起こして、たまに大きい値になる。これが、ゼロ交差判定の代替として、うまく働くのか、こっちの方が、真面目に計算したものより、結果が綺麗に見えるっぽい。あまり意味はないけど、参考までに実装した

また、(1)の話からすると、網膜では、エッジ強調が行われるらしいので、簡単なエッジ強調フィルタも通してみることにした。で、以下がコード

import cv2
import numpy as np
import math


def sharpen(img):
    kernel = np.array([[-1,-1,-1], [-1,9,-1], [-1,-1,-1]], np.float32)
    return cv2.filter2D(img, -1, kernel)


def DoG(gray , W, sigma1, sigma2 ,beta = 1.0):
    assert(W > 0)
    assert(gray.dtype==np.uint8)
    k1 = cv2.getGaussianKernel(W , sigma1)
    k2 = cv2.getGaussianKernel(W , sigma2)
    dog_kernel = np.dot(k1 , k1.T) - beta*np.dot(k2 , k2.T)
    res = cv2.filter2D(np.float64(img), -1, dog_kernel)
    minNN = cv2.morphologyEx(res, cv2.MORPH_ERODE, np.ones((3,3)))
    maxNN = cv2.morphologyEx(res, cv2.MORPH_DILATE, np.ones((3,3)))
    zeroCross = np.logical_or(np.logical_and(minNN >= 0 , maxNN <= 0), np.logical_and(minNN<=0 , maxNN>=0))
    return np.uint8(zeroCross)*255


def DoG_int(gray , W, sigma1, sigma2 ,beta = 1.0):
    assert(W > 0)
    assert(gray.dtype==np.uint8)
    k1 = cv2.getGaussianKernel(W , sigma1)
    k2 = cv2.getGaussianKernel(W , sigma2)
    G1,G2 = np.dot(k1,k1.T) , np.dot(k2,k2.T)
    res = cv2.filter2D(img, -1, G1) - cv2.filter2D(img , -1, G2)
    return res


def gabor(gray,lam=10.0):
   assert(gray.dtype==np.uint8)
   sigma = lam*0.4
   gamma = 1.0
   imgs = [cv2.filter2D(gray , -1 ,cv2.getGaborKernel((int(lam), int(lam)), sigma, np.radians(5*n), lam, gamma, 0)) for n in range(36)]
   ret = 0*gray
   for img in imgs:
      ret += img
   return np.clip(ret , 0, 255)


if __name__=="__main__":
   img = cv2.cvtColor(cv2.imread("Lenna.png"),cv2.COLOR_BGR2GRAY)
   contours = []
   contours.append( DoG(img, 5, 2.2 , 2.2*1.6) )
   contours.append( DoG(sharpen(img), 5, 2.2 , 2.2*1.6) )
   contours.append( DoG_int(img, 5, 2.2 , 2.2*1.6) )
   contours.append( DoG_int(sharpen(img), 5, 2.2 , 2.2*1.6) )
   contours.append( DoG(img, 7, 7*math.sqrt(0.2459)/math.pi , 7*math.sqrt(0.9837)/math.pi , 0.891) )
   contours.append( DoG(sharpen(img), 7, 7*math.sqrt(0.2459)/math.pi , 7*math.sqrt(0.9837)/math.pi , 0.891) )
   contours.append( DoG_int(img, 7, 7*math.sqrt(0.2459)/math.pi , 7*math.sqrt(0.9837)/math.pi , 0.891) )
   contours.append( DoG_int(sharpen(img), 7, 7*math.sqrt(0.2459)/math.pi , 7*math.sqrt(0.9837)/math.pi , 0.891) )
   contours.append( cv2.Canny(img,100,200) )
   contours.append( cv2.Canny(sharpen(img),100,200) )
   contours.append( gabor(img) )
   contours.append( gabor(sharpen(img)) )
   contour = np.concatenate([np.concatenate(contours[0:4], axis=1),np.concatenate(contours[4:8], axis=1), np.concatenate(contours[8:12], axis=1)] , axis=0)
   cv2.imwrite("result.png",contour)

結果発表。合計12パターンの計算をやってみた
f:id:m-a-o:20190922185224p:plain

1行目:左2つが、真面目に計算したDoGで、エッジ強調なしとあり。右2つは、オーバーフローするDoGで、エッジ強調なし、あり
2行目:パラメータを変えて、一行目と同じ組み合わせの計算をしたもの
3行目:1列目と2列目は、Cannyエッジフィルタ(エッジ強調なし、あり)。3列目と4列目は、Gaborフィルタを5度ずつ回転して足したもの(エッジ強調なし、あり)


まぁ、エッジ強調なしで、Canny検出器が一番いいように思える。Cannyは、エッジ強調すると、(当然だが)全然ダメになる。Gaborも割といい結果になってる気がする。こっちは、エッジ強調なしでもありでも、それなりの成績を出す。

肝心のDoGフィルタは、まぁ、どう解釈したものか、よく分からない。とりあえず、エッジ強調の影響が全然見られない。2行目の1,2列目とか見ると、全然取れてないってこともないので、パラメータ次第では、よい結果が出ることもあるのかもなーという感じではある。

以上は、生理学側からのアプローチで、既に形成された神経回路の機能を解析した話。



一方で、そもそも、こういう都合の良い神経回路が、どうやって自然に形成されるのかという問題もある。

脳内情報表現への情報理論的アプローチ
http://www.jnns.org/previous/niss/1999/text/sakakaba.pdf

によれば、「1970年代には，視覚野に見られる特徴抽出細胞が生後の感覚体験を通じて形成されることを示唆する実験報告が相次いだ」らしい。

これに対して、色々なモデルを作って、シミュレーションで再現するという研究が色々出たっぽい。上の総説は2000年頃のもののようだけど、現在でも、この問題が解決してるのかは、よく分からない。

大体、基本になってるのはHebb則で、加えて、いくつかの補助的な機構が仮定されるという感じで、網膜〜一次視覚野でHebb則による学習が存在するのかは知らないけど、Hebb則自体は、海馬なんかでのLTPの発見によって、脳神経の学習の適切なモデル化だと考えられている。


1980年代後半に、Linskerという人も、この手のシミュレーションをして、Linskerは、シミュレーション結果から一歩進んで、Hebb則による学習が何を最適化しているのかという視点から、infomax原理を提唱した。今では、"入出力間の相互情報量を最大化する"と説明されてることが多い

From basic network princples to neural architecture:Emergence of spatial-opponent cells
https://www.pnas.org/content/83/19/7508

From basic network princples to neural architecture:Emergence of orientation-selective cells
https://www.pnas.org/content/83/21/8390

From basic network princples to neural architecture:Emergence of orientation columns
https://www.pnas.org/content/83/22/8779

Self-organization in a perceptual network
https://ieeexplore.ieee.org/document/36
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.183.1750


Linskerの論文は、結局、主成分分析が出るみたいな話に見えるけど、主成分分析とEdge detectorやGabor filterが、どう繋がってるのかは定かでない

Linskerが、infomax原理からPCAを導くに当たって、二次の統計量までを考慮するという近似と、noiseはgaussianであるという性質を仮定したが、1996年の以下の論文では、その点を批判している。

The "independent components" of natural scenes are edge filters
https://www.ncbi.nlm.nih.gov/pubmed/9425547

論文の主な主張は、PCAではなく、ICAをやると、エッジフィルターが出るという話。PCAとICAの比較がされていて、PCAだと、エッジフィルターとは思えないものが得られている

論文でも、Gabor-likeという表現が使われているけど、本当にGaborなのか？とか、ICAだと何故Gaborが出るのか？とかは説明がない。この論文の結果を追試しようかと思ったけど、論文を読んでから長時間放置してたら、やる気がなくなったので、まぁ、そのうち

【補足】使用した画像データなんかは、既に入手できなくなってるっぽいけど、2013年の以下の論文でも同じようなことをやっていて、こっちはデータが現時点でまだ入手できるので、こっちのデータを使えば、まぁいいかという感じ

Double-Gabor Filters Are Independent Components of Small Translation-Invariant Image Patches
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3693455/


1990年代とかに、どうだったのか知らないけど、PCAを適用してみた事例は、今では沢山あって、まぁ概して、特にうまくいってるという感じもない。ICAだったら、どうかっていうのは、よく分からない。

同じ著者による論文で、自然音声に対して、同様のテクニックを適用したという論文が以下。

Learning the higher-order structure of a natural sound.
https://www.ncbi.nlm.nih.gov/pubmed/16754385

以下の論文も同じようなことをやってるっぽいけど、Gammatoneフィルター(?)っぽいものが得られるとか書いてある
Efficient coding of natural sounds.
https://www.ncbi.nlm.nih.gov/pubmed/11896400

Efficient auditory coding.
https://www.ncbi.nlm.nih.gov/pubmed/16495999

Gammatone filter
https://en.wikipedia.org/wiki/Gammatone_filter

以下は、二次視覚野もICAでイケソウというような話っぽい
Statistical model of natural stimuli predicts edge-like pooling of spatial frequency channels in V2.
https://www.ncbi.nlm.nih.gov/pubmed/15715907