沖縄返還50周年

今日は，1972年5月15日の沖縄の施政権返還から50周年，沖縄旅行の話は以前書いている。

1972年は，3月の3日から5日まで大学入試があり，4月からは大阪での生活がはじまった年だ。阪大物理学科の同窓会でも入学50年を迎えてという話があったがどうなることやら。

というわけで，1972年にあった様々なイベントは全て50年目を迎えることになる。そこで，Wikipediaをざっと調べて，自分の印象に残っているものをリストアップしてみた。千日前デパートビル火災と日中共同声明の記事が日本語版Wikipediaの1972年にはなかったが，英語版の1972年にはあった。

入試前の2月にいろいろな事件が集中していてたいへんだった。大学入試で金沢と大阪を往復した雷鳥の中で，パイオニア10号の記事をみて感動したのをおぼえている。カール・セーガンの発案で，地球外知的生命に向けてのメッセージを記録した金属板が搭載されていたからだ。

2月02日　横井庄一帰国
2月03日　札幌オリンピック開幕

2月21日　ニクソン訪中
2月28日　あさま山荘事件終結
3月01日　ローマクラブ「成長の限界」
3月02日　パイオニア10号打ち上げ
3月07日　連合赤軍リンチ殺人事件発覚（山岳ベース事件）
4月16日　川端康成ガス自殺
5月13日　千日前デパートビル火災
5月15日　沖縄返還・沖縄県発足
5月30日　テルアビブ空港乱射事件
9月05日　ミュンヘンオリンピック事件
9月29日　日中共同声明　日中国交正常化
12月7日　アポロ17号打ち上げ（アポロ計画終了）
12月28日　金日成国家主席就任

写真：パイオニア10号に搭載された知的生命体へのメッセージ金属板

［１］日本復帰への道（沖縄公文書館）

仏滅の周期

 トイレのカレンダーを見ていると，周期性のある仏滅の位置が今月末あたりでズレている。

六曜は，暦に記載される日時・方位などの吉凶、その日の運勢などの事項である暦注のひとつである。先勝（せんしょう）・友引（ともびき）・先負（せんぶ）・仏滅（ぶつめつ）・大安（たいあん）・赤口（しゃっこう）の六種がこの順に毎日繰り返す。

調べてみると，六曜の周期性の乱れは次のルールによる。旧暦の月の朔日が一月と七月は先勝，二月と八月は友引，三月と九月が先負，四月と十月が仏滅，五月と十一月が大安，六月と十二月が赤口と，あらかじめ定まった六曜が朔日に割り当てられており，これから順に並べていくことになる。そこで旧暦の月Mと日Dがわかれば，(M+D-2) mod 6 を計算しよう。この剰余の集合 { 0, 1, 2, 3, 4, 5 } に対して，先ほどの六曜の集合 {先勝, 友引, 先負, 仏滅, 大安, 赤口} を順番に対応させれば，各月各日の六曜が決まる。

今年の5月30日は旧暦の皐月朔日（5月1日）になるので，大安にセットされる。このため仏滅が，新暦の5/7−5/13−5/19/−5/25−5/31 という系列からずれて，6/4になったというわけだ。

あまり役に立たない豆知識シリーズ終了。

LDLコレステロール（１）

人間ドック（２）からの続き

以前から，中性脂肪やLDLコレステロール値にしばしば警告が出ていたが，あまり気にすることはなかった。ところが，今回はいよいよ，要治療領域まであと一歩という説明をドクターから直接受けたので早速以前の結果を見直してみた。 

図：人間ドックにおける2001年以降の血清脂質値

この22年の人間ドックの検査値の平均は中性脂肪（TG=トリグリセリド）が TG=113±24 mg/dL（35-149），LDLコレステロール（悪玉）が LDL-C=139±11 mg/dL（70-139），HDLコレステロール（善玉）が HDL-C=61±4 mg/dL （40-）となった。総コレステロールは，TC=LDL-C+HDL-C+TG/5で求まり，TC=225±12 mg/dL（130-219）である。括弧内が正常値だ。なお 善玉を除いたnon-HDLコレステロールは，non-HDL-C = TC − HDL-C となっている。

中性脂肪とHDLコレステロールはほぼ問題ないのだが，これまでのLDLコレステロールの平均値が，許容値の70-139 mg/dL の上限に達し，上回ることの方が多かった。LDLコレステロールが140 mg/dL以上は高コレステロール血症に該当し，180 mg/dL以上が要治療対象だ。加齢や運動量の不足からか，定年後の3-4年でLDL-CやTCが上昇傾向にあるのがよろしくない。

高コレステロール血症の 2-3割は食物や生活習慣が原因で，7-8割がコレステロールをつくる肝臓の機能の問題らしい。肝臓の機能のうち遺伝的な原因が主ならばなかなか対応が難しいわけだ。ネットで検索するとちょっと怪しいサプリページがどんどんと引っかかってくる。

とりあえず，次の検査の3ヶ月後まで，動物性脂肪（飽和脂肪酸）を少しでも減らして，有酸素運動を始めろということか。

［１］脂質異常症（eヘルスネット）

［２］人間ドック学会判定区分（2021年度版）

人間ドック（２）

人間ドック（１）からの続き

昨日は，天理市立メディカルセンターの人間ドック健診日だった。

天理市の国民健康保険加入者に対する人間ドック補助制度がある，今年から条件に該当することになったので，早速申し込んでいた。メディカルセンターの受付時間は8:30から8:45までで，10分前につくと6番目だ。大手前病院では，家を6:00前に出て，7:15に病院に到着してやっと6番目前後だったので，えらい違いである。

天理市立メディカルセンターは，高井病院の福祉医療法人である高清会が運営していて，設備も比較的新しく普通の病院とあまりかわらない。健診を受けている人数がそもそも少ないので対応はよい。タニタの最新式体重計では，体脂肪率だけでなく身体各部の脂肪量や筋肉量まで算出してくれる。脂肪量は平均値なのだが，筋肉量が平均より2段階下である。どおりで懸垂ができなくなるわけだ。

気になった点は，(1) 眼底検査がなかなかうまくいかず，３回取り直した。これは私のせいなのか？(2) 聴力検査はかなりアバウトなのであった。(3) 胃Ｘ線（バリウム）は検査技師の方はとても丁寧で親切だったけれど，発泡剤を補助液なしにそのまま口に入れ，大量のバリウムで直接流し込む方式には閉口した。診断台で身体を支えるのもなかなか困難になってきたので，来年は胃カメラかな。

問診は，超音波検査を担当していたお医者さんに血液検査のデータを見ながら丁寧に説明してもらった。LDLコレステロールが170mg/dLと後一歩で要治療の180mg/dLに達するということで脅された。さっそく看護師の保健指導をうけつつ3ヶ月後の再検査を申し込むことに。

［１］バックス発泡顆粒（「透視開始に際して、造影剤投与開始直前あるいは投与開始後、年齢、胃内容積の個人差、造影の体位に応じて、約100~400mLの炭酸ガスの 発生量に相当する量を、少量の水または、造影剤と共に経口投与する」なので，補助液なしでもかまわないらしい。）

理想気体のエントロピー

 （１）熱力学的なエントロピー$S(U,V)$を，$dS=\frac{d'Q}{T}\ $で定義する。これを熱力学第一法則の中に埋め込むと$\ dU=TdS-pdV$，すなわち，$dS=\frac{dU}{T}+\frac{p dV}{T}$ となる。

ここに，理想気体の状態方程式$\ pV = nRT = N k_B T\ $と，単原子分子気体の内部エネルギーの表式，$U = n C_v T = n \frac{3}{2} R T = N \frac{3}{2} k_B T\ $ を代入する。

$dS= N k_B \Bigl( \frac{3}{2} \frac{dT}{T} + \frac{dV}{V} \Bigr )\quad \therefore \int dS = N k_B \Bigl( \frac{3}{2} \int \frac{dT}{T} + \int \frac{dV}{V} \Bigr )$ 

$S = S_0 + N k_B \Bigl( \frac{3}{2} \log \frac{T}{T_0} + \log \frac{V}{V_0} \Bigr) = S_0 +  N k_B \log \frac{T^{3/2}V}{T_0^{3/2}V_0}$

（２）一方，統計力学において，自由粒子の小正準集団のエントロピーは，ボルツマンの原理から，$S = k_B \log W(E) = k_B \log \frac{\partial}{\partial E} \Omega_0(E)\  \delta E$となる。

質量$m$の$N$粒子系のエネルギー$E$までの状態数$\ \Omega_0(E)\ $は，$N$粒子系の$6N$次元の位相空間の体積を$6N$次元細胞の体積$\ h^{3N}\ $と同一粒子が区別できないことによる因子$N!$で割ったものになり，運動量空間での半径$\ \sqrt{2mE}\ $の$3N$次元超球の体積の式を使うと，

$\displaystyle \Omega_0(E) = \frac{1}{h^{3N} N!} \int d{\bm q} \int _{\Sigma p_i^2/2m < E} d{\bm p} = \frac{V^N (2\pi m E)^{3/2}}{h^{3N}\Gamma(3N/2+1)}$

$W(E) = \frac{3N}{2} \frac{V^N}{N! (3N/2)!} \Bigl( \frac{2\pi m E}{h^2} \Bigr)^{3/2} \frac{\delta E}{E} =  \frac{3N}{2} \Bigl(\frac{V}{N}\Bigr) ^N\Bigl( \frac{2\pi m E}{3N/2} \Bigr)^{3N/2} e^{5N/2} \frac{\delta E}{E}$

$\therefore S(E) = k_B \log W(E) = N k_b \Bigl\{\frac{3}{2}\log \Bigl( \frac{4\pi m E}{3 h^2 N}\Bigr)  + \log \frac{V}{N} +\frac{5}{2} \Bigr\}$

$=N k_B  \log \frac{(2\pi m k_B T)^{3/2} V e^{5/2}}{N h^3}$

ここで，$\frac{1}{T} = \frac{\partial S}{\partial U}= \frac{\partial S}{\partial E} = \frac{3 N k_B }{2  E}$より，$\frac{E}{N}=\frac{3}{2}k_B T$を用いた。$S(E)$の式の$\log$の中身は無次元であり，$V/N$があるので，全体として示量変数ではないことが保証されている。

理想気体の熱力学的エントロピーと統計力学的エントロピーは，ともに $N k_B \log T^{3/2}V + const.$の形をしているが，定数部分まで含めて同じかどうかがよく理解できていない。

摩擦のあるカルノーサイクル（３）

 摩擦のあるカルノーサイクル（２）からの続き

摩擦のあるカルノーサイクルでクラウジウスの不等式を説明するためには，前回のように$W'$のなかの摩擦力による仕事$\delta w$として表現するかわりに摩擦力で生じた熱$\delta_{\rm q}=-\delta w < 0$として扱うこともできる。

仕事として表現すると：
等温過程 A$\rightarrow$B：（$Q_{\rm H}>0, \quad \delta w_{\rm AB}>0$）
$W'_{\rm AB}=\int_{V_{\rm A}}^{V_{\rm B}} pdV -\delta w_{\rm AB}=W_{\rm AB} -\delta w_{\rm AB}=Q'_{\rm H}$
等温過程 C$\rightarrow$D：（$Q_{\rm L}<0, \quad \delta w_{\rm DC}>0$）
$W'_{\rm CD}=\int_{V_{\rm C}}^{V_{\rm D}} pdV -\delta w_{\rm DC}=W_{\rm CD} -\delta w_{\rm DC}= Q'_{\rm L}$
熱として表現すると：
等温過程 A$\rightarrow$B：（$Q_{\rm H}>0, \quad \delta q_{\rm H}<0$）
$W'_{\rm AB}=\int_{V_{\rm A}}^{V_{\rm B}} pdV + \delta q_{\rm H} = Q_{\rm H} +\delta q_{\rm H} = Q'_{\rm H}$
等温過程 C$\rightarrow$D：（$Q_{\rm L}<0, \quad \delta q_{\rm L} < 0$）
$W'_{\rm CD}=\int_{V_{\rm C}}^{V_{\rm D}} pdV + \delta q_{\rm L} = Q_{\rm L} + \delta q_{\rm L}= Q'_{\rm L}$

カルノーサイクルにおいては，系に入る熱量を温度でわった，エントロピーに対応する状態量$\frac{Q}{T}$の和が保存していた。すなわち，$\dfrac{Q_{\rm H}}{T_{\rm H}}+\dfrac{Q_{\rm L}}{T_{\rm L}} = n R \log \dfrac{V_{\rm B}}{V_{\rm A}} + n R \log \dfrac{V_{\rm D}}{V_{\rm C}} = 0$

一方，摩擦のあるカルノーサイクルで，出入りする熱量に対して，温度で割ったものの和を考えると，$\dfrac{Q'_{\rm H}}{T_{\rm H}}+\dfrac{Q'_{\rm L}}{T_{\rm L}}=\dfrac{Q_{\rm H}+\delta q_{\rm H}}{T_{\rm H}}+\dfrac{Q_{\rm L}+\delta q_{\rm L}}{T_{\rm L}} = \dfrac{\delta q_{\rm H}}{T_{\rm H}}+\dfrac{\delta q_{\rm L}}{T_{\rm L}} \le 0$

これを一般化すると，可逆過程だけのサイクルについては $\displaystyle{ \oint \dfrac{d'Q}{T_{\rm ex}} = 0}$，不可逆過程を含むサイクルについては，$\displaystyle{\oint \dfrac{d'Q}{T_{\rm ex}} \le 0}$，ここで，$d'Q$は系が受け取る熱量で，$T_{\rm ex}$はその熱量を与えた熱源の温度である。これがクラウジウスの不等式。

図：不可逆過程におけるクラウジウスの不等式

準静的過程がわからない

エントロピーがわからないからの続き

熱力学の入門的教科書を手元に並べて呻吟している。

そういえば，教養課程で物理学科の専門科目として大学１年のときにクラス担任の国富信彦先生が担当したのが「物理学要論？」だった（科目名も忘れてしまった）。そこで，最初につまづいたのが応力テンソルと準静的過程だった。熱力学の初歩のところでは，覆水盆に返らずの話をしながら準静的過程の説明があったので，これは可逆過程なのか不可逆過程なのかどうなっているの？と混乱したのだった。

さて，並べているやさしい教科書は以下のとおり

１．フェルミ熱力学（エンリコ フェルミ，三省堂 1973）

２．熱・統計力学（戸田盛和，岩波書店 1983）

３．熱・統計力学の考え方（砂川重信，岩波書店 1993）

４．熱学入門（藤原邦男・兵藤俊夫，東京大学出版会 1995）

５．ゼロからの熱力学と統計力学（和逹三樹・十河清・出口哲生，岩波書店 2005）

フェルミには準静的過程というワードは出てこない。それに相当するものは熱平衡状態をつないでいく可逆過程で考えるという立場だ。戸田さんは，この本（熱力学）で扱う可逆過程は準静的過程に限定すると注意している。砂川さんは，力学的な過程も可逆過程に含め，可逆過程は必ずしも時間反転してたどる必要がないとしている。準静的過程は可逆過程の集合に含まれる。藤原さんや和逹さんは，可逆過程＝準静的過程としている。

広島大学の戸田昭彦さんは準静的過程と可逆過程に対してもっと細かな議論を展開していた。

摩擦のあるカルノーサイクル（２）

 摩擦のあるカルノーサイクル（１）からの続き

「エントロピーについての理解を図るため，不可逆過程の具体的な例を構成したい。そのためカルノーサイクルの等温過程においてのみピストンに散逸のある抵抗力=摩擦力が働くモデルを考える。この摩擦力は，ピス トンの運動方向と逆向きに作用し，その仕事はピストンに熱として放出され，カルノーサイクル の作業物質である理想気体には影響を及ぼさないものとする。この考察において作業物質の系=理想気体がする仕事は，ピストンを用いて測定されることに留意する。すなわち，ピストンに働く力の総和とピストンの変位の積によって作業物質系が「する」仕事や「される」仕事（＝負の「する」仕事）が定義される」

ということで，前書きをかいて計算をはじめてみたもののなかなか難渋するのであった。

図：摩擦のあるカルノーサイクルの散逸過程

等温過程 A $\rightarrow$ B：（$\delta_{\rm AB}>0$）

$W'_{\rm AB}=\int_{V_{\rm A}}^{V_{\rm B}}p dV - \int_{d_{\rm A}}^{d_{\rm B}} f dx = nRT_{\rm H} \log \frac{V_{\rm B}}{V_{\rm A}} - \delta w_{\rm AB} = Q_{\rm H} - \delta w_{\rm AB} = Q'_{\rm H}$

断熱過程 B $\rightarrow$ C：

$W'_{\rm BC}=\int_{V_{\rm B}}^{V_{\rm C}}p dV = \int_{\rm B}^{\rm C} -dU = U_{\rm B}-U_{\rm C} = n C_V (T_{\rm H}-T_{\rm L})$

等温過程 C $\rightarrow$ D：（$\delta_{\rm DC}>0$）

$W'_{\rm CD}=\int_{V_{\rm C}}^{V_{\rm D}}p dV - \int_{d_{\rm C}}^{d_{\rm D}} f dx = nRT_{\rm L} \log \frac{V_{\rm D}}{V_{\rm C}} - \delta_{\rm DC} = Q_{\rm L} - \delta w_{\rm DC}= Q'_{\rm L}$

断熱過程 D $\rightarrow$ A：

$W'_{\rm DA}=\int_{V_{\rm D}}^{V_{\rm A}}p dV = \int_{\rm D}^{\rm A} -dU = U_{\rm D}-U_{\rm A} = n C_V (T_{\rm L}-T_{\rm H})$

したがって，この摩擦のあるカルノーサイクルの効率は次の式で与えられる。

$W'_c = W'_{\rm AB}+W'_{\rm BC}+W'_{\rm CA}+W'_{\rm AD}$

$= nRT_{\rm H}\log \frac{V_{\rm B}}{V_{\rm A}} -\delta_{\rm AB}  + nRT_{\rm L} \log \frac{V_{\rm D}}{V_{\rm C}} - \delta_{\rm DC} = Q_{\rm H} -\delta_{\rm AB} + Q_{\rm L} - \delta_{\rm DC}$

$\eta' = W'_c/Q'_{\rm H}=1-\frac{-Q'_{\rm L}}{Q'_{\rm H}} \approx 1-\frac{-Q_{\rm L}}{Q_{\rm H}} \bigl( 1+\frac{\delta w_{\rm DC}}{-Q_{\rm L}} \bigr) \bigl(1+\frac{\delta w_{\rm AB}}{Q_{\rm H}} \bigr)$

$\le 1-\bigl( \frac{-Q_{\rm L}}{Q_{\rm H}}\bigr) =  1-\frac{T_{\rm L}}{T_{\rm H}} = \eta_c$

こうして，摩擦のあるカルノーサイクルの効率$\eta'$は，カルノーサイクルの効率$\eta_c$より小さくなる。しかし，このままでは，$\int \frac{d'Q}{T} \le 0$の説明にうまくつながらない。

（注）ここではサイクルに入る熱量をすべて正，サイクルが外部にする仕事を全て正にとる。

摩擦のあるカルノーサイクル（１）

エントロピーがわからないからの続き

熱力学の第二法則と仮に導入したエントロピーの違いを明らかにしたい。普通の教科書ではクラウジウスの原理やトムソンの原理によるわけだが，そのためには，不可逆過程の考察が必要となる。その一番簡単な例は，力学でもなじみのある摩擦現象だと思う。

摩擦現象は力学的に細かく詰めて考えると何だか複雑で面倒なことになるが，とりあえず，摩擦現象は運動や仕事が熱に変換されて，力学的エネルギーが熱エネルギーとして環境中に散逸する過程だと考えることにする。環境中の熱エネルギーが直接集まってきて力学的エネルギーになるような現象は，巨視的には観察されないので，摩擦現象は不可逆過程である。

そこで，これまで練習してきたカルノーサイクルに摩擦を導入すれば，理解がしやすいのではないかと考えた。どの教科書をみてもそんな具体的な議論はされていない。クラウジウスやトムソンにしたがって，より抽象的な熱機関（可逆機関，不可逆機関）の組み合わせでの議論が進んでいくわけだ。

というわけで，より具体的な摩擦のあるカルノーサイクルを構成してみることに（続く）。

エントロピーがわからない

カルノーサイクルからの続き

エントロピーについての熱力学的な導入の論理がすっきりしないと，統計力学の授業が進めにくい。もちろん天下りでボルツマンの原理を導入してしまえばあとは計算だけになる。でも，それでは熱力学との関係もうやむやになりそうだ。

熱力学の第一法則で，$dU=d'Q+d'W=d'Q-pdV$のd'Qの部分も状態量の組み合わせで書けるとありがたい。$p$は示強変数，$V$は示量変数であり，その積がエネルギーの次元を持つ示量変数になっている。使える状態量として示強変数である温度$T$があるので，これに相補的な示量変数で温度との積がエネルギーの次元を持つ状態量をエントロピー$S$として導入して，$d'Q=TdS$とおくことにする。

もしこれができれば，状態量空間中の点を$A$，基準点を$O$として，$S(A)=\int_O^A \frac{d'Q}{T}$は状態量になる。この積分が状態量であるということは，平衡状態Aのみに依存して積分の経路にはよらないはずである。

そこで，カルノーサイクルの断熱過程で実際にこの量を計算してみれば，断熱過程ではエントロピー$S$が一定になる。つまり，カルノーサイクルというのは，エントロピーと温度を2軸とする状態図において，等温線と等エントロピー線に囲まれた長方形領域になる。

ここまでの議論は，準静的過程＝可逆過程について成り立つ話である。不可逆過程だとどうなるのか。肝腎の熱力学の第二法則との関係がついていないわけなのでさらなる検討が必要だ。

クローニッヒ・ペニーモデル（４）

クローニッヒ・ペニーモデル（３）からの続き

これまでの結果を用いて，具体的な数値を代入したグラフを作成する。与えられた$k_n$の値に依存して，$-1 \le \cos(k_n a) \le 1$の条件から，1電子のエネルギー$E$がとりうる範囲についての条件が定まる。このとき，電子が取り得るエネルギー領域を許容帯（allowed band），取りえないエネルギー領域を禁制帯（forbidden band）とよぶ。相互作用のない電子の多粒子系では，これらの離散化された許容帯に電子が充填されることになる。

Mathematicaによって，この様子を確認してみれば次のようになる。

a = 2; b = 0.04; hc = 1973; m = 1.022*10^6; p = 2*1.05017;
v = 2*hc^2*p/(1.022*10^6*a*b)
(hc)^2/(m a^2)
200.001
0.952233

\[Alpha][e_] := Sqrt[m*e]*(a - b)/hc;
\[Beta][e_] := Sqrt[m*(v - e)]*b/hc;
\[Gamma][e_] := Sqrt[m*e]*a/hc;
X[e_] := Cos[\[Alpha][e]] Cosh[\[Beta][ e]] + ((\[Beta][e]/b)^2 - (\[Alpha][e]/(a - b))^2)/(2*\[Alpha][ e]/(a - b)*\[Beta][e]/b) Sin[\[Alpha][e]] Sinh[\[Beta][e]]

Plot[{ X[e], 1, (Cos[\[Gamma][e]] + p*Sin[\[Gamma][e]]/\[Gamma][e]), -1, Cos[\[Gamma][e]]}, {e, 0, 200}, PlotRange -> {-1.5, 3.5}, PlotStyle -> {Red, Gray, Blue, Gray, Orange}]
Table[FindRoot[X[e] == 1, {e, (hc a n )^2/(2 m) }], {n, 1, 5}]
{{e -> 2.95059}, {e -> 37.6031}, {e -> 45.1073}, {e -> 150.412}, {e -> 158.267}}
Table[ FindRoot[X[e] == -1, {e, 0.9*(hc a n )^2/(2 m) }], {n, 1, 6}]
{{e -> 9.40077}, {e -> 16.0036}, {e -> 84.6069}, {e -> 92.3771}, {e -> 242.886}, {e -> 242.886 + 0. I}}

d0 = Plot[{10^6 (e - 2.95059), 10^6 (e - 9.40077), 10^6 (e - 16.0036), 10^6 (e - 37.6031), 10^6 (e - 45.1073), 10^6 (e - 84.6069), 10^6 (e - 92.3771), 10^6 (e - 150.421), 10^6 (e - 158.267)}, {e, 0, (4 Pi)^2}, PlotStyle -> Table[{Gray, Dotted}, 9], PlotRange -> {0, 13}];
f0 = Plot[{0, Pi, 2 Pi, 3 Pi, 4 Pi}, {e, 0, (4 Pi)^2}, PlotStyle -> Table[{Green, Dotted}, 5]];
g0 = Plot[\[Gamma][e], {e, 0, (4 Pi)^2}, PlotStyle -> {Orange, Dashed}];
g1 = Plot[ArcCos[X[e]], {e, 0, Pi^2}];
g2 = Plot[2 Pi - ArcCos[X[e]], {e, Pi^2, (2 Pi)^2}];
g3 = Plot[2 Pi + ArcCos[X[e]], {e, (2 Pi)^2, .98 (3 Pi)^2}];
g4 = Plot[4 Pi - ArcCos[X[e]], {e, .98 (3 Pi)^2, .95 (4 Pi)^2}];

Show[{d0, f0, g0, g1, g2, g3, g4}, PlotRange -> {-1, 15}]

図１：境界条件から決まる（energy-cos(ka)）のグラフ

図２：図１から決まる分散関係（energy-ka）のグラフ

ポテンシャルの幅$b$を格子周期$a$の2%に設定したため，$b \rightarrow 0$の近似がよく当てはまる。そこで，$X(pa)$と$\cos \gamma(e) + P \sin \gamma(e)/\gamma(e)\ $の誤差は1%以下となり，グラフ上では区別できなくなっている。

なお，$e=\frac{(\hbar k)^2}{2m}=\frac{(\hbar c)^2}{2mc^2 a^2}\cdot (ka)^2 = 0.952233 (ka)^2\ $であることに注意する。

クローニッヒ・ペニーモデル（３）

クローニッヒ・ペニーモデル（２）からの続き

1次元ポテンシャルに周期性があるときに，ブロッホの定理から$\psi(x)=e^{ikx}\varphi(x)$と表わせて，$\psi(x+a)=e^{ik(x+a)}\varphi(x+a)=e^{ika} e^{ikx}\varphi(x)=e^{ika}\psi(x)$が成り立つ。このときの波動関数は運動量演算子の固有状態なのだろうか？違います。前回やったように，このハミルトニアンは有限の並進操作に対して不変だけれど，運動量に対応する無限小並進操作については不変ではないから。

ところで，この長さ$L=N a$の1次元周期ポテンシャルモデルの両端を同一視する周期境界条件をつけると（$N$はポテンシャルステップの数＝原子数，$a$はポテンシャルの周期＝原子間隔），$\psi(L)=\psi(0) \quad \psi(L)=e^{i k a \cdot N}\psi(0) \quad \therefore e^{i k a N}=1$

これから$k$に対する条件，$k_n = \frac{2\pi n}{a N}\quad (n=0,\pm 1, \pm 2 \cdots)\ $が得られる。$k_n$は量子数 $n$ で特徴づけられるこの状態の波数という意味をもつ。

前回得られた境界条件は，系のエネルギーを$E$，ポテンシャルの深さと幅を$V_0, b$，ポテンシャル周期を$a$として，$p=\frac{\sqrt{2mE}}{\hbar}$，$q=\frac{\sqrt{2m(V_0-E)}}{\hbar}$とおくと，$\cos k a  =  \cos p(a-b) \cosh qb + \frac{q^2-p^2}{2 p q} \sin p(a-b) \sinh qb$ である。これは，与えられた$k = k_n$に対して，系のエネルギーを決定する式になる。

(1) $b \rightarrow 0$ の極限では$p_n=k_n$となり，$E_n=\frac{\hbar k_n^2}{2m} = \frac{2 \hbar^2 \pi^2 n^2}{m a^2 N^2}$となる。

(2) 次に，$V_0 b$を一定に保ちながら，$b \rightarrow 0,\ V_0 \rightarrow \infty$とするδ関数型極限を考える。このとき，$\sinh qb \rightarrow qb$であり，右辺第2項は，$\frac{(q^2-p^2)ba}{2} \frac{\sin p(a-b)}{p a}$となる。最終的に，$\cos k a  =  \cos p a + \frac{m c^2 V_0 b a}{(\hbar c)^2} \frac{\sin pa}{pa}$ という近似式が得られる。

例えば，$a=2$ Å，$b=0.04$ Å，$mc^2 = 0.511 \times 10^6$ eV，$\hbar c =1973$ eV Å， $V_0=100$ eVとすると， 無次元のポテンシャル強度パラメータは，$\frac{m c^2 V_0 b a}{(\hbar c)^2}=1.05$となる。

日本國憲法前文

憲法記念日だが，ウクライナへのロシアの侵略戦争が続いている。先の見えない円安で疲弊した我々の心の隙に，右翼デマゴーグ達の好戦的なイデオロギーが陽に陰に染み込んでいく。こんなときは，日本國憲法の前文を写経して心を鎮めるしかない。若者は「最高法規の意志～ 憲法の本質と改正の動向 ～」をチラ見する方が役に立つかも。

　　　日本國憲法 

　日本國民は、正當に選擧された國會における代表者を通じて行動し、われらとわれらの子孫のために、諸國民との協和による成果と、わが國全土にわたつて自由のもたらす惠澤を確保し、政府の行爲によつて再び戰箏の惨禍が起ることのないやうにすることを決意し、ここに主權が國民に存することを宣言し、この憲法を確定する。そもそも國政は、國民の嚴肅な信託によるものであつて、その權威は國民に由來し、その權力は國民の代表者がこれを行使し、その福利は國民がこれを享受する。これは人類普遍の原理であり、この憲法は、かかる原理に基くものである。われらは、これに反する一切の憲法、法令及び詔勅を排除する。

　日本國民は、恒久の平和を念願し、人間相互の關係を支配する崇高な理想を深く自覺するのであつて、平和を愛する諸國民の公正と信義に信頼して、われらの安全と生存を保持しようと決意した。われらは、平和を維持し、專制と隷從、壓迫と偏狹を地上から永遠に除去しようと努めてゐる國際社會において、名譽ある地位を占めたいと思ふ。われらは、全世界の國民が、ひとしく恐怖と缺乏から免かれ、平和のうちに生存する權利を有することを確認する。

　われらは、いづれの國家も、自國のことのみに專念して他國を無視してはならないのであつて、政治道徳の法則は、普遍的なものであり、この法則に從ふことは、自國の主權を維持し、他國と對等關係に立たうとする各國の責務であると信ずる。

　日本國民は、國家の名譽にかけ、全力をあげてこの崇高な理想と目的を達成することを誓ふ。

クローニッヒ・ペニーモデル（２）

クローニッヒ・ペニーモデル（１）からの続き

前回，周期ポテンシャル中で正のエネルギーを持った電子の運動について考えた。自由電子とはいえ，金属中に束縛されているのだからポテンシャルの上端に対して負のエネルギーを持った電子を考えなければならなかった。

上智大学名誉教授の清水清孝さん（元有馬研）が，電子情報通信学会の知識ベース知識の森の12群5編量子力学・電子物理・相対論を執筆していて，そこにバンド理論入門についての話もあった。

そこで，前回のポテンシャルの符号を $-V_0 \rightarrow V_0$としたモデルで$0<E<V_0$の場合を考える。すなわち，ポテンシャルは，つぎの形を繰り返したものになる。
領域 Ⅰ：$V(x) =  \ 0 \ \cdots \ (0 < x < a-b)\ $
領域 Ⅱ：$V(x) = V_0 \ \cdots \ (-b < x < 0)$

自由電子のエネルギーは正だが，ポテンシャルの高さより小さいので，領域Iでは解は平面波，領域IIでは指数関数になる。$p=\frac{\sqrt{2mE}}{\hbar}, q=\frac{\sqrt{2m(V_0-E)}}{\hbar}$とおいて一般解は，

領域 Ⅰ：$\psi_{\rm I}(x) = A e^{i p x} + A' e^{-i p x} \ \cdots \ (0 < x < a-b)\ $

領域 Ⅱ：$\psi_{\rm II}(x) =  B e^{q x} + B' e^{- q x} \ \cdots \ (-b < x < 0)\ $

前回のように波動関数とその導関数が，領域Iと領域IIの境界で連続である条件から，

$\psi_{\rm I}(0) = \psi_{\rm II}(0);\ \psi_{\rm I}'(0) = \psi_{\rm II}'(0)$，$\psi_{\rm I}(a-b) = e^{i k a }\psi_{\rm II}(-b);\ \psi_{\rm I}'(a-b) = e^{i k a }\psi_{\rm II}'(-b)$

$A+A'=B+B'$
$i p(A-A')=q(B-B')$
$A\ e^{ip(a-b)}+A'e^{-ip(a-b)}=(B\ e^{-qb}+B'e^{qb})e^{ika}$
$i pA\ e^{ip(a-b)}-i pA'e^{-ip(a-b)}=(qB\ e^{-qb}-qB'e^{qb})e^{ika}$

ここで，$\alpha = p(a-b)$，$\beta = q b$と置くと，

$\left( \begin{array}{cccc} 1 & 1 & -1 & -1 \\ i p & -i p & -q & q \\ e^{i\alpha} & e^{-i\alpha} & -e^{-\beta}e^{ika} & -e^{\beta}e^{ika} \\ i p e^{i\alpha} & - i p e^{-i\alpha} & -q e^{-\beta}e^{ika} & q e^{\beta}e^{ika} \end{array} \right)  \left( \begin{array}{c} A\\ A' \\ B \\ B' \end{array} \right) = 0$

4元連立方程式が自明でない解を持つためには，行列式が0であり，列の加減により等価な行列式を求める（3行目から$e^{ika}$を，2列目から$i$を外にくくり出した）。

$i e^{ika} \cdot \left| \begin{array}{cccc} 1 & 0 & -1 & 0 \\ 0 & p & 0 & q \\ \cos \alpha \ e^{-ika} & \sin \alpha \  e^{-ika}& -\cosh \beta & - \sinh \beta \\ -p\ \sin \alpha & p\ \cos \alpha & q\ \sinh \beta \ e^{ika} & q\ \cosh \beta \ e^{ika} \end{array} \right|  = 0$

これから4x4行列式の値$|M|$を計算する。

$|M|=  p\ q\ e^{i k a} \left | \begin{array}{cc} -\cosh \beta  & - \sinh \beta  \\ \sinh \beta & \cosh \beta \end{array} \right | + q \left | \begin{array}{cc}  \sin \alpha e^{-i k a } & - \cosh \beta \\  p \cos \alpha & q \sinh \beta e^{ika} \end{array} \right |$

$\qquad \quad - q\ p\ e^{-i k a}\ \left | \begin{array}{cc} \cos \alpha  &  \sin \alpha \\ - \sin \alpha & \cos \alpha  \end{array} \right | - p \left | \begin{array}{cc} - \sinh \beta  & \cos \alpha e^{-ika} \\ q \cosh \beta e^{ika} & -p \sin \alpha \end{array} \right |$

$\quad = -2 p q \cos ka + q( q \sin \alpha \sinh \beta + p \cos \alpha \cosh \beta ) -p (p \sin \alpha \sinh \beta -q \cos \alpha \cosh \beta )$

$\quad = -2 p q \cos ka +2 p q \cos \alpha \cosh \beta +(q^2-p^2) \sin \beta \sinh \beta = 0$

最終的に得られる関係式は，次のとおりである。

$\cos ka  = \cos \alpha \cosh \beta -\frac{p^2-q^2}{2 p q} \sin \alpha \sinh \beta$

$\qquad = \cos p(a-b) \cosh qb -\frac{p^2-q^2}{2 p q} \sin p(a-b) \sinh qb$

（付）この式は，$E>V_0$の場合は次のように変形できるので，このまま使うことができる。$q=\frac{\sqrt{2m(V_0-E)}}{\hbar} \rightarrow i \frac{\sqrt{2m(E-V_0)}}{\hbar} \equiv i\bar{q}\ $として，$\sinh ix = i \sin x,\ \cosh ix = \cos x$を用いると，

$\cos ka  =  \cos p(a-b) \cos \bar{q}b -\frac{p^2 + \bar{q}^2}{2 p \bar{q}} \sin p(a-b) \sin \bar{q} b\ $となる。これは（１）で得られた式と同等なもの。

図：クローニッヒ・ペニーモデル

クローニッヒ・ペニーモデル（１）

クローニッヒ・ペニーモデルは周期性を持った１次元井戸型ポテンシャルのモデルであり，結晶のバンド構造の定性的な特徴を説明することができる。

図：クローニッヒ・ペニーモデルの設定（Wikipediaより引用）

電子の質量を$m$として，1粒子の1次元シュレーディンガー方程式は，$-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\psi(x) +V(x) \psi(x) = E \psi(x)$である。ポテンシャルは周期性を持ち，$n$を整数として$V(x+n a)=V(x)\ $であり，ブロッホの定理から波動関数は，$\psi(x) = e^{i k x} \varphi(x)$ かつ $\varphi(x + n a)=\varphi(x)$を満たす。

ポテンシャルは，つぎの形を繰り返したものになる。

領域 Ⅰ：$V(x) =  \ 0 \ \cdots \ (0 < x < a-b)\ $

領域 Ⅱ：$V(x) = -V_0 \ \cdots \ (-b < x < 0)$

波動関数の周期性からは，$\psi(x+a) = e^{i k (x+a)}\varphi(x+a) = e^{i k (x+a)}\varphi(x) = e^{i k a} \psi(x) \ $が成り立つ。その導関数は$\psi'(x+a) = e^{i k a}\psi '(x)$となる。なお，$\psi'(x) = i k e^{i k x}\varphi(x) + e^{i k x}\varphi'(x) = i k \psi(x) + e^{i k x}\varphi'(x)$である。

自由電子のエネルギーは正であり，ポテンシャルが定数であるため，どちらの領域でも解は平面波になる。$p=\frac{\sqrt{2mE}}{\hbar}, q=\frac{\sqrt{2m(E+V_0)}}{\hbar}$として一般解は，

領域 Ⅰ：$\psi_{\rm I}(x) = A e^{i p x} + A' e^{-i p x} \ \cdots \ (0 < x < a-b)\ $

$\qquad \quad = e^{i k x} (A e^{i (p-k) x} + A' e^{-i (p+k) x} )$

領域 Ⅱ：$\psi_{\rm II}(x) =  B e^{i q x} + B' e^{-i q x} \ \cdots \ (-b < x < 0)\ $

$\qquad \quad = e^{i k x} (B e^{i (q-k) x} + B' e^{-i (q+k) x})$

波動関数とその導関数が，領域Iと領域IIの境界で連続であるという条件を書く。
$\psi_{\rm I}(0) = \psi_{\rm II}(0);\ \psi_{\rm I}'(0) = \psi_{\rm II}'(0)$，$\psi_{\rm I}(a-b) = e^{i k a }\psi_{\rm II}(-b);\ \psi_{\rm I}'(a-b) = e^{i k a }\psi_{\rm II}'(-b)$

$A+A'=B+B'$
$p(A-A')=q(B-B')$
$A\ e^{ip(a-b)}+A'e^{-ip(a-b)}=B\ e^{-iqb+ika}+B'e^{iqb+ika}$
$pA\ e^{ip(a-b)}-pA'e^{-ip(a-b)}=qB\ e^{-iqb+ika}-qB'e^{iqb+ika}$

ここで，$\alpha = p(a-b)$，$\beta = q b$と置く。

$\left( \begin{array}{cccc} 1 & 1 & -1 & -1 \\ p & -p & -q & q \\ e^{i\alpha} & e^{-i\alpha} & -e^{-i\beta}e^{ika} & -e^{i\beta}e^{ika} \\ p e^{i\alpha} & - p e^{-i\alpha} & -q e^{-i\beta}e^{ika} & q e^{i\beta}e^{ika} \end{array} \right)  \left( \begin{array}{c} A\\ A' \\ B \\ B' \end{array} \right) = 0$

波動関数の係数に対するこの4元連立方程式が自明でない解を持つためには，行列式が0でなければならない。この条件をつかって，2列と1列の和と差，4列と3列の和と差から等価な行列式を求めれば，次のようになる（3行目から$e^{ika}$を外にくくり出した）。

$e^{ika} \cdot \left| \begin{array}{cccc} 1 & 0 & -1 & 0 \\ 0 & p & 0 & q \\ \cos \alpha \ e^{-ika} & i \sin \alpha \  e^{-ika}& -\cos \beta & -i \sin \beta \\ i p\ \sin \alpha & p\ \cos \alpha & i q\ \sin \beta \ e^{ika} & q\ \cos \beta \ e^{ika} \end{array} \right|  = 0$

これから4x4行列式の値$|M|$を計算する。

$|M|= p\ q\ e^{i k a} \left | \begin{array}{cc} -\cos \beta  & - i \sin \beta  \\ i \sin \beta & \cos \beta \end{array} \right | + q \left | \begin{array}{cc} i \sin \alpha e^{-i k a } & - \cos \beta \\ p \cos \alpha & i q \sin \beta e^{ika} \end{array} \right |$

$\qquad \quad - q\ p\ e^{-i k a}\ \left | \begin{array}{cc} \cos \alpha  & i \sin \alpha \\ i \sin \alpha & \cos \alpha  \end{array} \right | - p \left | \begin{array}{cc} -i \sin \beta  & \cos \alpha e^{-ika} \\ q \cos \beta e^{ika} & i p \sin \alpha \end{array} \right |$

$\quad = -2 p q \cos ka + q( -q \sin \alpha \sin \beta + p \cos \alpha \cos \beta ) -p (p \sin \alpha \sin \beta -q \cos \alpha \cos \beta )$

$\quad = -2 p q \cos ka +2 p q \cos \alpha \cos \beta -(p^2+q^2) \sin \beta \sin \beta = 0$

最終的に得られる関係式は，次のとおりである。

$\cos ka  = \cos \alpha \cos \beta -\frac{p^2+q^2}{2 p q} \sin \alpha \sin \beta$

$\qquad = \cos p(a-b) \cos qb -\frac{p^2+q^2}{2 p q} \sin p(a-b) \sin qb$

4行4列の行列式は，Mathematicaを使えば手軽に計算できるのだけれど，手計算でもなんとかなる場合があるということを学ぶ。

１次元周期ポテンシャル

 月曜日の予習シリーズ。

１次元の周期ポテンシャル中を運動する粒子の問題を考える。

(1) 並進演算子：$\psi(x+\delta x) \approx \psi(x) + \delta x \cdot \frac{d}{dx} \psi(x)  = (1 + i \ \delta x \cdot p_x / \hbar ) \psi(x)$から，運動量演算子は微小並進操作と関係している。そこで，ユニタリ演算子，$U(a) = \exp( i\ a \cdot p_x / \hbar )$が，有限の並進操作を行う演算子となる。つまり，$U(a) \psi (x) = \sum_{k=0}^\infty \frac{1}{k!} (\frac{i\ a \cdot p_x}{\hbar})^k \psi(x) =  \sum_{k=0}^\infty \frac{a^k}{k!} (\frac{d}{dx})^k \psi(x) = \psi(x+a)$

(2) １次元周期ポテンシャル：1次元のポテンシャル$V(x)$中を運動する質量$m$の粒子に対する定常状態のシュレーディンガー方程式は，$H \psi(x) = \{ -\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+V(x) \}\psi(x)  = E \psi(x)$ である。このポテンシャルが周期$a$を持つとき，すなわち，$V(x+a)=V(x)$のとき，$U(a)V(x)\psi(x)=V(x+a)\psi(x+a)=V(x) U(a)\psi(x)$なので，$[U(a), V(x)]=0$，また，$U(a)$は演算子$p$から構成されるので，$[U(a),\frac{p^2}{2m}]=0$である。

(3) 固有関数の並進対称性：したがって，$[U(a),H]=$であり，$H$の固有関数は，$U(a)$との同時固有関数（絶対値1の複素固有値）になるから，$U(a) \psi(x)  = \exp(ika)\psi(x)$とかける。つまり，$\psi(x+a) = \exp(ika) \psi(x)$であり，$\psi(x)=\exp(ikx) \phi(x)$とすると，$\phi(x+a)=\phi(x)$を満足することになる。すなわち，ブロッホの定理「周期ポテンシャルの固有関数は同じ周期性を持つ関数と平面波の積となる」が成り立つ。

うーん，ここからクローニッヒ＝ペニーモデルに持ち込むにはちょっと覚悟が必要だということがわかったので，宿題にする。

フェルミ分布

月曜の授業の予習シリーズ。

位相空間（$\mu$空間）の細胞に含まれる状態数は，プランク定数を$h$として，$dn=\frac{1}{h^3} dx dy dz dp_x dp_y dp_z$である。電子のようなスピン1/2のフェルミ粒子を考えると，位相空間の各状態にスピンアップとダウンの2状態がともなう。そこで，単位体積をとって，運動量空間細胞に含まれる状態数は $dn' = \frac{2}{V_0} \int_V dn = \frac{2}{h^3}dp_x dp_y dp_z  = \frac{8 \pi}{h^3} p^2 dp$となる。ただし，$p^2=p_x^2+p_y^2+p_z^2$。粒子の質量を$m$，エネルギーを$w=\frac{p^2}{2m},\ p=\sqrt{2mw}$とすると，$dw=\frac{p}{m}dp\ $より，$dn'=\frac{8 \pi}{h^3} \sqrt{2m^3 w}\ dw\ $となる。

この系のフェルミ分布関数は，$f(w_i)=[\exp(\frac{w_i - w_F}{kT}) + 1]^{-1}$である。ただし，$k$はボルツマン定数，$T$は系の絶対温度，$w_F$はフェルミ準位を表わす。これから，エネルギーの分布関数は，$N(w)dw= \frac{8 \pi}{h^3} \sqrt{2m^3 w} [\exp(\frac{w - w_F}{kT}) + 1]^{-1} dw$となる。これを速度空間の表式にひきもどして，$v_x,v_y$で積分すると$\ v_z$の分布関数が得られる。

そこで，次の関係式に留意する。$\int_{-\infty}^{\infty}dv_x\int_{-\infty}^{\infty}dv_y\int_{-\infty}^{\infty}dv_z f(v_x,v_y,v_z) = 4\pi \int_{0}^{\infty} f(v^2)\ v^2 dv$

$w=\frac{m}{2}v^2$，$dw = m v\ dv$なので

$4\pi \int_{0}^{\infty}f(v^2)  v^2 dv = 4\pi \int_{0}^{\infty} f(v^2)  \frac{v}{m} dw =  \int_{0}^{\infty} 4\pi f(v^2)  \sqrt{\frac{2w}{m^3}} dw$

そこで，$N(w)dw =  4\pi f(v^2)  \sqrt{\frac{2w}{m^3}} dw\ $とおけば，$f(v^2)=N(w) \frac{1}{4\pi} \sqrt{\frac{m^3}{2w}}\ $となるので，$\ v_z$の分布関数 $N(v_z) dv_z$は次式で与えられる。

$N(v_z) dv_z = \int_{-\infty}^{\infty}dv_x\int_{-\infty}^{\infty}dv_y f(v_x,v_y,v_z) dv_z =  \int_{-\infty}^{\infty}dv_x\int_{-\infty}^{\infty}dv_y N(w) \frac{1}{4\pi} \sqrt{\frac{m^3}{2w}} dv_z$

$= \int_{-\infty}^{\infty}dv_x\int_{-\infty}^{\infty}dv_y \frac{2m^3}{h^3} [\exp(\frac{w - w_F}{kT}) + 1]^{-1} dv_z$

$v_x, v_y$平面での積分を2次元の極座標によって実行するため，$u^2=v_x^2+v_y^2$とおいて，$dv_x dv_y = 2\pi u du$となる。

$\therefore \quad N(v_z) dv_z = \frac{4\pi m^3}{h^3} \int_0^\infty du u  [\exp(\frac{\frac{m}{2}(u^2+v_z^2) - w_F}{kT}) + 1]^{-1} dv_z$

$= \frac{4\pi m^3}{h^3} \int_0^\infty \frac{1}{2} dt  [\exp(\frac{\frac{m}{2}(t+v_z^2) - w_F}{kT}) + 1]^{-1} dv_z$

ここで，$a=\exp(\frac{\frac{m}{2}v_z^2 - w_F}{kT})$，$b=\frac{m}{2kT}$とおけば，必要な積分は$\int_0^\infty \frac{1}{a \exp(bt) + 1}dt$となり，その値は$\ \frac{1}{b}\log(1+1/a)\ $である。これより，$N(v_z) dv_z = \frac{4\pi m^2 kT}{h^3} \log (1 + \exp(\frac{w_F - \frac{m}{2}v_z^2}{kT})) dv_z$

これらの分布関数をMathematicaでプロットすると次のようになる。

f[w_, kT_] := Sqrt[w]/(Exp[(w - 1)/kT] + 1)
Plot[Table[f[w, 0.01*k], {k, 1, 10, 2}], {w, 0, 2},PlotRange -> {0, 1}]
g[v_, kT_] := kT/2 Log[(Exp[(1 - .5*v^2)/kT] + 1)]
Plot[Table[g[v, 0.01*k], {k, 1, 10, 2}], {v, 0, 2},PlotRange -> {0, 0.5}]

図１：エネルギー分布関数　N(w)

図２：速度分布関数　N(v_z)

モル比熱

かつて 中学校で熱について学んだとき，もっとも重要な基本法則は熱量と温度と比熱の関係だった。これが重要であることは，大学の熱力学でもそうなのだけれど，あくまでも熱力学第一法則と第二法則の脇役であって，電磁気学のオームの法則のようなものだ。

物質量が$\ n\ $モルの体系に熱量$\ d'Q\ $を与えたときに，温度が$\ dT\ $だけ増えたとする。系の温度を1K上げるために必要な熱量である熱容量$\ {\rm [J/K]}\ $は，$\frac{d'Q}{dT}$で与えられる。このとき，系の体積を一定にするならば定積熱容量 $C_V$，系の圧力を一定にするならば定圧熱容量 $C_p$ とよぶ。これらは物質量に比例する示量変数である。

熱力学の第一法則より $\ d'Q = dU + pdV = dU + d(pV)-V dp\ $が成り立つ。したがって，$C_V = \frac{dU}{dT}$，$C_p=\frac{dU}{dT} + \frac{d(pV)}{dT}$となる。ここで，理想気体を考えると，状態方程式 $\ pV = n R T\ $が成り立ち，$C_p=C_V + n R$と表わされる。

単位質量あるいは単位物質量あたりの熱容量が比熱容量＝比熱となる。定積モル比熱は$c_V=\frac{1}{n}C_V$，定圧モル比熱は$c_p=\frac{1}{n} C_p$と小文字の$c$で表わすことになるが，教科書を眺めると，そのあたりの定義や記号の使い方は必ずしもそろっているわけではなかった。

カルノーサイクル

熱力学の復習シリーズ，カルノーサイクルの練習をする。

　熱力学第一法則： $dU = d'Q + d'W = d'Q - p dV$

　理想気体の状態方程式： $pV=nRT$

　理想気体のポアソンの法則： $pV^\gamma = const,\quad T V^{\gamma-1} = const'$

　エントロピー： $dU = TdS -p dV,\quad dS = \frac{dU + pdV}{T}$

　内部エネルギー： $U = nC_{V}T$

■過程 A $\rightarrow$ B（$p_{\rm A}V_{\rm A}=p_{\rm B}V_{\rm B}$）

理想気体が高温熱源$T_{\rm H}$と接触を保ちつつ，一定の温度$T_{\rm H}$の状態を保ちつつ，熱量$Q_{\rm H}$をもらって膨張し，外へ仕事$W_{\rm AB}$をする。理想気体の温度は一定なので，内部エネルギーは$U_{\rm B}=U_{\rm A}$であり，熱力学第一法則より$W_{\rm AB}=Q_{\rm H}$である。

外部にした仕事は，$W_{\rm AB}=\int_{V_{\rm A}}^{V_{\rm B}}p dV = \int_{V_{\rm A}}^{V_{\rm B}}\frac{nRT_{\rm H}}{V} dV=nRT_{\rm H}\log \frac{V_{\rm B}}{V_{\rm A}} = Q_{\rm H}$

エントロピー変化は，$S_{\rm AB}= \int_{\rm A}^{\rm B} dS = \int_{V_{\rm A}}^{V_{\rm B}} \frac{p}{T} dV = \int_{V_{\rm A}}^{V_{\rm B}} \frac{nR}{V} dV = nR \log \frac{V_{\rm B}}{V_{\rm A}} = \frac{Q_{\rm H}}{T_{\rm H}}$

■過程 B $\rightarrow$ C（$p_{\rm B}V_{\rm B}^\gamma=p_{\rm C}V_{\rm C}^\gamma \quad T_{\rm H} V_{\rm B}^{\gamma-1} = T_{\rm L} V_{\rm C}^{\gamma-1}$）

断熱壁と接触する理想気体が，熱の流入なしに断熱的に膨張して外に仕事$W_{\rm BC}$をする。熱力学第一法則によって，理想気体の内部エネルギーは$U_{\rm B}$から$U_{\rm C}$まで減少し，温度は$T_{\rm L}$まで下がる。熱の出入りがないのでエントロピーは変化しない。

外部にした仕事は，$W_{\rm BC}=\int_{V_{\rm B}}^{V_{\rm C}}p dV = \int_{\rm B}^{\rm C} -dU = U_{\rm B}-U_{\rm C} = U(T_{\rm H}) - U(T_{\rm L})= n C_V (T_{\rm H}-T_{\rm L})$

■過程 C $\rightarrow$ D（$p_{\rm C}V_{\rm C}=p_{\rm D}V_{\rm D}$）

理想気体が低温熱源$T_{\rm L}$と接触を保ちつつ，一定の温度$T_{\rm L}$の状態を保ちつつ，熱量$Q_{\rm L}$を放出して収縮し，外から仕事$W_{\rm CD}$がなされる。理想気体の温度は一定なので，内部エネルギーは$U_{\rm D}=U_{\rm C}$であり，熱力学第一法則より$W_{\rm CD}=Q_{\rm L}$である。

外部からされる仕事は，$W_{\rm CD}=\int_{V_{\rm C}}^{V_{\rm D}}-p dV = \int_{V_{\rm C}}^{V_{\rm D}}-\frac{nRT_{\rm L}}{V} dV=nRT_{\rm L}\log \frac{V_{\rm C}}{V_{\rm D}} = Q_{\rm L}$

エントロピー変化は，$S_{\rm CD}= \int_{\rm C}^{\rm D} dS = \int_{V_{\rm C}}^{V_{\rm D}} \frac{p}{T} dV = \int_{V_{\rm C}}^{V_{\rm D}} \frac{nR}{V} dV = nR \log \frac{V_{\rm D}}{V_{\rm C}} = - \frac{Q_{\rm L}}{T_{\rm L}}$

■過程 D $\rightarrow$ A（$p_{\rm D}V_{\rm D}^\gamma=p_{\rm A}V_{\rm A}^\gamma \quad T_{\rm L} V_{\rm D}^{\gamma-1} = T_{\rm H} V_{\rm A}^{\gamma-1}$）

断熱壁と接触する理想気体を，熱の流入なしに断熱的に圧縮して外から仕事$W_{\rm BC}$がされる。熱力学第一法則によって，理想気体の内部エネルギーは$U_{\rm C}$から$U_{\rm D}$まで増加し，温度は$T_{\rm H}$まで上がる。熱の出入りがないのでエントロピーは変化しない。

外部からされる仕事は，$W_{\rm DA}=\int_{V_{\rm D}}^{V_{\rm A}}- p dV = \int_{\rm D}^{\rm A} dU = U_{\rm A}-U_{\rm D} = U(T_{\rm H}) - U(T_{\rm L})= n C_V (T_{\rm H}-T_{\rm L})$

■カルノーサイクルの効率

1サイクルの過程${\rm A \rightarrow B \rightarrow C \rightarrow D \rightarrow A}$において，理想気体（作業物質）が外部にする正味の仕事は，$W = W_{\rm AB} + W_{\rm BC} - W_{\rm CD} -W_{\rm DA} = W_{\rm AB} - W_{\rm CD}$

$= nRT_{\rm H}\log \frac{V_{\rm B}}{V_{\rm A}} - nRT_{\rm L} \log \frac{V_{\rm C}}{V_{\rm D}} = Q_{\rm H}-Q_{\rm L}$

このカルノーサイクルの効率は $\eta = \frac{W}{Q_{\rm H}} = \frac{Q_{\rm H}-Q_{\rm L}}{Q_{\rm H}} = 1 - \frac{Q_{\rm L}}{Q_{\rm H}}$で与えられる。

ところで，ポアソンの法則の温度と体積の関係式を組み合わせると，

$\frac{T_{\rm H}}{T_{\rm L}}=\Bigl( \frac{V_{\rm C}}{V_{\rm B}}\Bigr)^{\gamma-1}=\Bigl(\frac{V_{\rm D}}{V_{\rm A}}\Bigr)^{\gamma-1}$，

したがって，$\frac{V_{\rm C}}{V_{\rm B}}=\frac{V_{\rm D}}{V_{\rm A}} \quad \frac{V_{\rm A}}{V_{\rm B}}=\frac{V_{\rm D}}{V_{\rm C}}$

$\therefore \frac{Q_{\rm L}}{Q_{\rm H}}=\frac{nRT_{\rm L}\log \frac{V_{\rm C}}{V_{\rm D}}}{nRT_{\rm H}\log \frac{V_{\rm B}}{V_{\rm A}}}=\frac{T_{\rm L}}{T_{\rm H}}, \quad \eta = 1 - \frac{T_{\rm L}}{T_{\rm H}}$

■カルノーサイクルのエントロピー

$S = S_{\rm AB} + S_{\rm BC} +S_{\rm CD} +S_{\rm DA} = \frac{Q_{\rm H}}{T_{\rm H}} -  \frac{Q_{\rm L}}{T_{\rm L}} = nR \log \frac{V_{\rm B}}{V_{\rm A}} + nR \log \frac{V_{\rm D}}{V_{\rm C}} = 0$

カルノーサイクルでは，${\rm A} \rightarrow {\rm B}$の等温膨張過程で熱を吸収するとともに，理想気体のエントロピーが増加し，${\rm C} \rightarrow {\rm D}$の等温圧縮過程で熱を放出するとともに，理想気体のエントロピーが減少する。その結果，1サイクルが終了後にはエントロピーの増減はなくなり，エントロピーが状態量であることが保証されている。

図：カルノーサイクルのp-V図

ミクロカノニカル分布

 ミクロカノニカル分布について。

小正準集団（ミクロカノニカルアンサンブル）とは，外界から孤立した系の熱平衡状態を記述するための統計集団。考えている孤立系と粒子数（$N$），体積（$V$）が等しく，エネルギー（$E$）がある幅（$\delta Ｅ$）の範囲で等しい系（コピー）の集団である。その数は，下記で定義される微視的な状態数$W$に等しい。

考えている$\ N\ $粒子系の位相空間（$\Gamma\ $空間）とする。1粒子の位相空間（$\mu\ $空間）の体積が$h^f$で表わされるとき，$\Gamma\ $空間の細胞$\ d\Gamma\ $における微視的な状態数を$\ dW = \frac{1}{h^{Nf}} d\Gamma = \frac{1}{h^{Nf}} dq_1 \cdots dq_N\ dp_1 \cdots dp_N$とする。これから，$W= \frac{1}{h^{Nf}} \int_{E}^{E+\delta E} dq_1 \cdots dq_N\ dp_1 \cdots dp_N$

等重率の原理は，巨視的に観測される全エネルギーが$E$である小正準集団がすべておなじ重みで確率的な平均操作に寄与するというものである。この確率分布を小正準分布（ミクロカノニカル分布）とよぶ。このとき，ある物理量$A(q_1 \cdots q_N, p_1 \cdots p_N)$において，これを小正準分布を用いてその観測される期待値を求めると次式のようになる。

$\langle A \rangle = \int_{E}^{E+\delta E} A(q_1 \cdots q_N, p_1 \cdots p_N) dq_1 \cdots dq_N\ dp_1 \cdots dp_N / \int_{E}^{E+\delta E} dq_1 \cdots dq_N\ dp_1 \cdots dp_N$

等重率の原理を用いて，$N$粒子系の全位相空間の点を同じ確率で扱う根拠として，かつては，エルゴート定理をその根拠とする教科書が多かった。ところが，田崎晴明さんの統計力学の教科書（2008）でこれを否定してからは，こうした教科書は少なくなった。もっとも，高橋康さんの統計力学入門（1984）には，そのあたりはていねいに書いてあったのだった。

愛の設計図

松竹座で5月から始まる「藤山寛美三十三回忌追善喜劇特別公演」に関連したイベントがなんばパークスシネマで開催された。応募ハガキが当選したので授業前の時間を利用して大阪まで出る。

藤山寛美による1983年の松竹新喜劇公演の「愛の設計図」がDVD化されていて，これを上映したあとに，渋谷天外と藤山扇次郎のトークショウがあった。DVDの音源を利用しているからか，劇場の音響がちょっときつすぎて閉口した。

曽我廼家文童が情報処理プログラマーの試験に受かるところからはじまって，建設会社の現場監督と設計士を巡るドラマが進んでいく。アドリブも沢山入っているらしい藤山寛美のセリフは必ずしも全編に渡って流ちょうというわけでもないのだけれど，要所要所で観客の心をつかむ技はさすがである。

ラザフォード

 月曜日の授業の課題を考えていた。

先週のオリエンテーションの後，ゴールデンウィーク明けまでは，現代物理学の歴史の概論を講義する予定なので，現代物理学の基礎を築いた物理学者について調べてもらうことにする。で，Wikipediaの日本語版と英語版の記述量がかなり違うことを知ってもらうというのを目的として，英語版にはあって日本語版にない情報を選んで簡単に紹介してもらうという趣旨にした。

そのため，物理学者の一覧を作るべく調べていると，ノーベル物理学賞の受賞者にラザフォードの名前がない。原子の構造にたどり着いた一番肝腎な人物なのに，いったいどうしたことでしょう。というわけで，トイレからでて落ち着いてノーベル化学賞の方に，ソディやデバイなどとともに無事にリストされていた。

ちなみに，宿題に出したのは60名のリストからさらに精選した以下の24名とした。

ヴィルヘルム・レントゲン　W. C. Rontogen (1845-1923)
アルバート・マイケルソン　A. A. Michelson (1852-1931)
アンリ・ベクレル　A. H. Becquerel (1853-1908)
ジョゼフ・ジョン・トムソン　J. J. Thomson (1856-1940)
マックス・プランク　M. Planck (1858-1947)
ピエール・キュリー　P. Curie (1859-1906)
ヘンリー・ブラッグ　W. H. Bragg (1862-1942)
フィリップ・レーナルト　P. Lenard (1862-1947)
ビーター・ゼーマン　P. Zeeman (1865-1943)
マリ・キュリー　M. Curie (1867-1934)
ロバート・ミリカン　R. A. Millikan (1868-1953)
ジャン・ペラン　J. B. Perrin (1870-1942)
アーネスト・ラザフォード　E. Rutherford (1871-1937)
フレデリック・ソディ　F. Soddy (1877-1956)
マックス・フォン・ラウエ　M. von Laue (1879-1960)
ジェイムス・フランク　J. Franck (1882-1964)
マックス・ボルン　M. Born (1882-1970)
ピーター・デバイ　P. Debye (1884-1996)
ニールス・ボーア　N. Bohr (1885-1962)
エルヴィン・シュレーディンガー　E. Schrodinger (1887-1961)
アーサー・コンプトン　A. Compton (1892-1962)
ルイ・ド・ブロイ　L. V. de Broglie (1892-1987)
ヴォルフガング・パウリ　W. E. Pauli (1900-1958)
ボール・ディラック　P. A. M. Dirac (1902-1984)

原子の寿命

現代物理学の導入部分で， 古典物理学（ニュートン力学＋マクスウェル電磁気学）では原子の安定性が説明できないことがひとつの鍵になる。このためには，加速度運動する電子が電磁波を放出することを示す必要がある。ところが，加速荷電粒子からの電磁波の放出は，電磁気学で学ぶ最終コーナーであり，そこまで到達しない場合が多い。

お茶の水女子大学の理学部3年次編入試験では，この部分が次元解析で説明されていた。最初に，陽子（$+e$）の周りを円運動している電子（$-e$）の位置エネルギー$V(r)$を無限遠点を基準として求める。ただし，$r$は陽子から電子までの距離。クーロン定数は$k=\frac{1}{4\pi\varepsilon_0}$とするので，$V(r) = - k \frac{e^2}{r}$。

次に，この非相対論的な運動をしている電子（速度$\ v$，質量$\ m_e$）の全エネルギー$\ E\ $を求めると$\ E=\frac{1}{2}m v^2 - k \frac{e^2}{r}$。なお，円運動の向心力＝クーロン力から，$m_e \frac{v^2}{r} = k \frac{e^2}{r^2}\ $を用いると，$E = -\frac{k}{2} \frac{e^2}{r}\ $であり，加速度は$\ a = \frac{v^2}{r} = \frac{k e^2}{m_e r^2}\ $。

加速度運動する電子から単位時間に放出される電磁波のエネルギー $S\ {\rm [kg \cdot m^2 \cdot s^{-3}]}$を，次元解析によって表わす。電子の加速度$\ a\ {\rm [m \cdot s^{-2}]}$，微細構造定数$\ \alpha\ $を使って$\ e^2\ k=\alpha \hbar c\ {\rm [kg \cdot m^3 \cdot s^{-2} ]}$ ，光速度$\ c\ {\rm [m \cdot s^{-1}]} \ $を用いると，$S \sim \alpha \hbar c \frac{a^2}{c^3}$ となる。以下では数係数を1とする。

円運動する電子が単位時間に失うエネルギーは，$-\frac{d E}{d t} = \frac{k e^2}{2} \frac{-\dot{r}}{r^2}=  \frac{\alpha \hbar c}{2} \frac{-\dot{r}}{r^2}$。これが上記の$S$と等しいことから，$\frac{\alpha \hbar c}{2} \frac{-\dot{r}}{r^2} = \alpha \hbar c \frac{a^2}{c^3}$，つまり，$\dot{r} = - 2 r^2 \frac{a^2}{c^3}$

ここで，$\frac{a}{c} = \frac{\alpha \hbar c}{m_e c r^2}\ $なので，$r^2 \dot{r} = - 2  \bigl( \frac{\alpha \hbar c}{m_e c^2} \bigr)^2 c\ $

初期状態で半径$r_0$の原子が電磁波を放出して半径0になるまでの時間を$\tau$とすると，

$\int_{r_0}^0 r^2 dr = \int_0^\tau - 2  \bigl( \frac{\alpha \hbar c}{m_e c^2} \bigr)^2 c\ dt\ $から，

$-\frac{r_0^3}{3} = - 2  \bigl( \frac{\alpha \hbar c}{m_e c^2} \bigr)^2 c \tau$

$\therefore \tau = \frac{r_0^3}{6 c} \bigl( \frac{m_e c^2}{\alpha \hbar c} \bigr)^2$

$r_0=10^{-10} {\rm [m]}$，$\alpha = \frac{1}{137}$，$\hbar c = 197 \times 10^{-15}{\rm [MeV \cdot m]}$，$m_e c^2 = 0.5 {\rm [MeV]}$，$c=3 \times 10^8 {\rm [m \cdot s^{-1}]}$ を代入すると，$\tau = 6.7 \times 10^{-11} {\rm [s]}$となる。

ファン・デル・ワールスの状態方程式

統計物理学の授業がはじまった。専門科目にしては受講者が多かったので，講義室は共通講義棟の１Ｆに割り当てられていた。授業の第１法則：受講者数は時間とともに指数関数的に減少する。授業の第２法則：受講者数の空間密度分布は教卓からの距離の逆数に比例して減少する。

最後にファン・デル・ワールスの状態方程式を紹介して授業を終えたところ，早速質問があった。理想気体では，$pV=n\ R\ T\ $だった状態方程式が，実在気体の場合$\bigl( p + \frac{a}{V^2} \bigl) (V - b) = R\ T\ $と書いてあるけれど，右辺にモル数の$\ n\ $が抜けているのかというものだ。

「あ，ごめんごめん，忘れてたわ」と返事して終ったものの，帰り道に380段の階段を下りながら考えてみると何だかおかしい。圧力の補正項$\ \frac{a}{V^2}\ $は分子間力によるものであり，気体密度の二乗に比例する項だと説明した。ところが，圧力は示強変数なのに，補正項の分母は示量変数の二乗になっている。つまり，$a\ $は定数ではなく示量変数の二乗に比例しなければならない。

演習課題の参考にしていた「熱学入門」（藤原邦男・兵頭俊夫）の式も同様の問題点を含んだままだった。そこで，いくつかの本などでファン・デル・ワールス方程式を調べてみると，気体のモル数を1モルに限定していたり，体積として$\ V\ $ではなく，モル当たり体積$\ V_m\ $を用いている。

ということで，正しくは，$n$モルの実在気体に対しては，$\bigl( p + \frac{a n^2}{V^2} \bigl) (V - n b) = n\ R\ T\ $としなければならない。$a, b\ $の値もこれまで考えたこともなかったが，お茶の水女子大学の理学部編入試験問題にも採用されているくらいであり，水（H2O）の場合，$a=5.54\times 10^{-1} {\rm [Pa \cdot m^6 \cdot mol^{-2}]}$，$b=3.05 \times 10^{-5} {\rm [m^3 \cdot mol^{-1}]}\ $だった。

これを使ってMathematicaで状態図をプロットしてみる。

図：水蒸気のファン・デル・ワールス状態図

a = 0.554; b = 3.05*10^-5; R = 8.314; 8 a/(27 R b)

647.331　・・・　（臨界温度）

p[V_, T_] := R T /(V - b) - a/(V^2)

q[V_, T_] := Abs[R T /(V - b) - a/(V^2)]

LogLogPlot[{q[V, 673], q[V, 647.33], q[V, 573], q[V, 473], q[V, 373], 

  q[V, 273]}, {V, 3*10^-5, 3*10^-4}, PlotRange -> {10^5, 10^9}]

Plot[{p[V, 673], p[V, 647.33], p[V, 573], p[V, 473], p[V, 373], 

  p[V, 273]}, {V, 3*10^-5, 3*10^-4}, PlotRange -> {-10^8, 10^8}]

渦 妹背山婦女庭訓 魂結び：大島真寿美

 最近，読書のスピードが極めて遅くなっている。本日，ようやく大島真寿美の「渦 妹背山婦女庭訓 魂結び」を読了した。2019年の第161回直木賞受賞作品だ。

もうすこし，重くて硬い作品かと予想して読みはじめると，意外に軽い文章だったので，始めのうちは少し慣れなかった。やがて，近松半二の物語についていくことができるようになった。この作品は操り浄瑠璃の竹本座の座付作者である近松半二（1725-1783）をテーマに選んだところが鍵だったと思う。

物語としては，やや物足りなかった。250年前に近松半二は，日高川入相花王（1759），奥州安達原（1762），本朝廿四孝（1766），傾城阿波の鳴門（1768），近江源氏先陣館（1769），妹背山婦女庭訓（1771），新版歌祭文（1780），伊賀越道中双六（1783）と，現在でも頻繁に上演されている，三大名作と近松の世話物以外の主要な浄瑠璃を次々と産み出した。それにも関わらず，操り浄瑠璃が歌舞伎に破れていくところが最も印象に残ったところだ。

もう一つは，妹背山婦女庭訓の主要登場人物のお三輪をふくめて，登場する女性陣がしなやかな強さをもっているところか。半二の妻はお佐久で，娘はおきみなのだが，これがどのようにして伊賀越道中双六の作者として名を連ねた近松加作につながっていくのかいかないのかは，次作の「結 妹背山婦女庭訓 波模様」を待たなければならないのか。

岡本綺堂の戯曲「近松半二の死」は青空文庫でも読める小編だが，そこでは祇園町の娘お作が加作につながることが暗示されている。

写真：大島真寿美の渦の書影（amazonnより引用）

abc予想（４）

 abc予想（３）からの続き

abc予想が正しければ，ただちにフェルマーの定理が簡単に証明されるとのことだが，事情は若干ややこしい。正確には，abc予想の(2)の表現において$\varepsilon = 1$として，このときに，$K(\varepsilon)=1$が満たされればという条件がついてくる。ウェブ上ではこれを強いabc予想と書いているものもあるが，その表現がよいのかどうかはあまりはっきりしない。

(強いabc予想？)　$c < K(\varepsilon=1)  {\rm rad}(abc)^{1+(\varepsilon=1)} = 1 \cdot {\rm rad}(abc)^2$ 

これからフェルマーの定理を導くのは次のようになる。

自然数$a,\ b,\ c\ $が互いに素であり，自然数$n > 6\ $に対して，$a^n,\ b^n,\ c^n\ $が abcトリプレットをなすとする。すなわち，$a^n + b^n = c^n\ $であって，$a^n < b^n$かつ$a^n,\ b^n\ $は互いに素。このとき，$c^n < {\rm rad}(a^n b^n c^n)^2 = {\rm rad}(a b c)^2 < (a b c)^2 < c^6\ $が成り立つ。

つまり，$n \ge 6\ $に対して，$a^n + b^n = c^n\ $は偽であることから，フェルマーの定理が $n \ge 6\ $に対して成り立つ。$n= 3, 4, 5\ $については別途成り立つことを示す必要があるが，これらはすでに証明されている。

これを少し詳しく見た様子は，黒川信重・小山信也による「ABC予想入門」に記述されている。

abc予想（３）

abc予想（２）からの続き

$c > rad(abc)\ $が成り立つabcトリプルを例外的abcトリプルとよぶことにする。そうはいっても前回は，これが無限個存在することも示している。例外的トリプルの個数を計算してみると，$c<1,000\ $で31個（全トリプル151,895個），$c<10,000\ $で120個（全トリプル個15,196,785），$c<100,000\ $で418個（全トリプル1,519,805,376）となる（最大10^15のオーダーまでの数の素因数分解を10億回しなければならないので，めちゃ時間かかった）。

これを計算したMathematicaコードは次の通り。

Rad[n_] := Times @@ (First /@ FactorInteger[n]);
f[n_] := {k = 0; Do[If[Rad[a*b*(a + b)] < a + b && GCD[a, b, a + b] == 1, k = k + 1], {a, 1, n/2 - 1}, {b, a + 1, n - a - 1}]; Print[k]};
g[n_] := {k = 0; Do[If[GCD[a, b, a + b] == 1, k = k + 1], {a, 1, n/2 - 1}, {b, a + 1, n - a - 1}]; Print[k]};
f[1000]
31
g[1000]
151 895

そこで，少しだけ条件を変更することで，例外的トリプルが有限個になるようにできないかを考えたところ，次のような予想が立てられた。(1)と(2)は同等であり，これがabc予想といわれるものである。

(1) 任意の$\varepsilon > 0$ に対して $c > {\rm rad}(abc)^{1+\varepsilon}\ $を満足するabcトリプルは高々有限個しか存在しない。

(2) 任意の$\varepsilon > 0$ に対してある $K(ε) > 0$ が存在し、全ての abcトリプルについて$c < K(\varepsilon) {\rm rad}(abc)^{1+\varepsilon}\ $が成り立つ。

abc予想（２）

 abc予想（１）からの続き

肝腎のabc予想である。準備のために，自然数 $n$ の根基（radical）を定義する。素数を$\{ p_i \}$として，$n = p_1^{k_1} \cdot p_2^{k_2} \cdot \cdots \ $と素因数分解できるとき，$n$ の根基が次のように定義される。${\rm rad}(n) \equiv p_1 \cdot p_2 \cdots \ $ 。例えば，${\rm rad}(120)={\rm rad}(2^2  \cdot 3^2 \cdot 5) = 2 \cdot 3 \cdot 5 = 30$ 。

3つの自然数の組，$(a, b, c)$ において，$a+b=c\ $であり，$a < b\ $が互いに素であるものを，abc-トリプルとよぶ。普通は，$c < {\rm rad}(a b c)$ であるが，$c > {\rm rad}(a b c)$も無限に存在しうる。そこで，ある形のabc-トリプルで後者の条件を満たすものが無限個あることを示す。

素数 $p>2$，自然数 $n>1$として，$a=1, b=2^{p(p-1)n}-1, c=2^{p(p-1)n}$というabc-トリプルを考える。このとき，${\rm rad}(a b c) = {\rm rad}(a) {\rm rad}(b) {\rm rad}(c) = 1 \cdot {\rm rad}(b)  \cdot 2 = 2\ {\rm rad}(2^{p(p-1)n}-1)$となる。

ところで，$b = 2^{p(p-1)n} - 1 = \bigl( 2^{p(p-1)} \bigr)^n - 1 = \bigl( 2^{p(p-1)}-1 \bigr) \bigl( \cdots  \bigr) = p^2 \cdot \bigl( \cdot \  \cdots  \bigr)$ とかけるので(*)，$\frac{b}{p^2}$は自然数である。そこで，${\rm rad}(a b c) = 2 \cdot {\rm rad} (p^2 \cdot \frac{b}{p^2}) \le 2 \cdot p \cdot  \frac{b}{p^2} = 2 \cdot \frac{b}{p} \lt \frac{2}{p} c$ 。ただし，${\rm rad}(m) \le m$ を用いた。

(*) フェルマーの小定理により，$2^{p-1}\ \equiv 1 \ ({\rm mod}\ p)$，すなわち，$2^{p-1} = k\ p + 1\ $なので，$2^{p(p-1)} = (k\ p + 1)^p = \Sigma_{r=0}^{p-1} C_r^p\ (k\ p)^{p-r} 1^r +1\ $。このため，$2^{p(p-1)}-1=(k\ p + 1)^p - 1\ $を展開した $k p$の最低次の項は，$C_{p-1}^p \cdot (k\ p) = k\ p^2$であり，全体は$p^2$を因子として持つ。

つまり，このabc-トリプルは任意の$n$に対して，$c > \frac{2}{p} c > {\rm rad}(abc)$ となるので，無限個存在する。

abc予想（１）

 NHKスペシャルの「数学者は宇宙をつなげるか？abc予想証明を巡る数奇な物語」の完全版（90分）が放映された。60分版にエピソードをいくつか加えて丁寧に編集していた。60分版を見た水野義之さんによれば，「素因数という言葉もさけて人間ドラマに仕立てた大失敗作だ」ということだが，そうでもないと思うけど。

Wikipediaのabc予想の記載は，現状をほぼ客観的に説明している。英語版のほうが淡白（冷淡）で，日本語版はややていねいに経緯を示している。NHKスペシャルでは，望月新一さん本人は取材を断っており，加藤文元さん，プリンストンでの望月の指導教官のファルティングス，批判派のピータ・ウォイトなどの話が中心となった。

番組でタイトルだけ紹介されたのが望月さんのブログからのキーワード「ノーと言える人間」であり，その説明はなかった。調べてみると，新一の「こころの一票」というブログの2017年11月21日に「心壁論」と，論理構造の解明・組合せ論的整理術を「心の基軸」 とすることの本質的重要性という記事があった。

abc予想を証明することができる宇宙際タイヒミュラー理論は，京大数理解析研究所のPRIMSで査読が完了して出版されたものの，世界の数学者の大勢は，その証明を認めていない。そんなこともあって，望月さんのこころの内がブログには表出されている。

望月さんは，5歳で日本を離れてからほぼ米国で生活し，16歳でプリンストンに入学し，23歳でPh. Dをとっている。そんな彼が欧米の大学のポストではなく日本の京大に職を得たことに理由に相当する，英語圏での生活における強い違和感が詳しく述べられている。

許容される表現のイメージ全体に一つの固定された「座標系」を敷き、表現のイメージ全体を通して　同一の「座標系」＝「視点」＝「声」＝「神」＝「心の基軸」しか認めないという姿勢を徹底することにはどうしても強い違和感を覚える

「ノーと言える人間」 というのは，例の盛田昭夫と石原慎太郎の1989年の著書「NO」と言える日本―新日米関係の方策」からきている。これが，欧米のメインストリームの数学者を権威主義的に敬う空気への強烈な拒否感として現れたものにほかならない。

P. S. 新一の「こころの一票」にNHKスペシャルについての評価に関する記事が掲載された。なかなかおもしろい（2022.5.3）。

にほんごであそぼ

 4月のNHKの番組改編では例年より大幅なものとなった。その中でも，Eテレのにほんごであそぼの評判が極めて悪かったので，録画を飛ばし飛ばしちょっとだけチェックしてみた。

うーん，頭にポットをつけたお姉さんが，百発百中の説明というかゲームに無駄な時間を費やすだけで時間が過ぎてゆく。部分的に以前のフレーバーが残っているものの，伝統芸能のかけらもなく，貴重だった織太夫，清介，勘十郎の文楽のコーナーも見当たらなかった。

1979年に福音館書店から出版された「にほんご」は，安野光雅，大岡信， 谷川俊太郎，松居直によって，小学校低学年の国語の教科書に代わるものとしてデザインされた労作だった。子供が小学校にあがるまで，ふとんの中でよく読んだものだ。そのフレーバーやコンセプトがふんだんに盛り込まれたのが，NHKの「にほんごであそぼ」だったが，番組改編で大変残念なことにそれが失われてしまった。

NHKはニュース系がひどくても，Eテレがあるからと我慢していた。この調子だと本当に全部がダメになりそうで怖い。まあ，あの優等生的で擬似中立主義的な総体による隠れた洗脳が問題だというのであれば，Eテレの優良番組だってどうなのよということかもしれないが。

写真：福音館のにほんごの書影（福音館書店から引用）

Wボソンの質量

4月7日付のScienceに，CERMのCDF実験によるWボソン質量の精密測定の結果が報告された。アブストラクトを訳すと次のようになる。

素粒子間の弱い力の媒介となるWボソンの質量は，素粒子物理学の標準模型の対称性によって厳しい制約を受けている。ヒッグス粒子は，この標準模型の最後の欠落部分であった。ヒッグス粒子が観測された後，Wボソンの質量を測定することで，標準模型を厳密に検証することができる。

私たちは，フェルミ国立研究所のテバトロン加速器のCDFII検出器を用いて，重心エネルギー1.96 TeVの陽子・反陽子衝突で収集した8.8 /fb の積分光度に相当するデータを用いてWボゾンの質量M_Wを測定した。約400万個のWボゾン候補のサンプルを用いて，M_W = 80,433.5 ± 6.4 stat ± 6.9 syst = 80,433.5 ± 9.4 MeV/c^2 を求め，その精度はこれまでのすべての測定値を合わせたものを超える。この測定は、標準模型の予想と大きく異なる（注：標準理論では，M_W = 80,357 ± 6 MeV/c^2，7σの食い違い）。

標準理論の綻びは，ミューオンの異常磁気能率 のところにも現れていたので，いよいよ新しい物理が見えそうな期待感が高まる。ただ，野尻さんは論文のアペンディックスを見て懐疑的な感想を述べているので，まだどうなるかわからない。Wボソンの質量とニューオンの異常磁気能率の組み合わせが，標準理論の矛盾を大きくする方向に働いているようだ。

早速，様々な理論が話題になっているけれど，レプトクォークよりヒッグス粒子のダブレットの方が無難なのかもしれない。

 

写真：CDF 検出器（FNALから引用）

うさみ亭マツバや

大阪のきつねうどんの起源には諸説あるが，うさみ亭マツバやが元祖だというのがなんとなく有力らしい。一度行こうと思っていたので，午後の文楽のついでによってみた。

心斎橋筋を北に向かい，ちょっと右に折れると間口の狭いこじんまりした店があった。11:00開店で，11:30ごろについたが，既にお客さんがかなり入っている。早速きつねうどんとミニ親子丼のセットを注文した。

感想：味は普通。麺はやや太めでもっちゃりする。うどんの断面が丸くて箸から滑りました。油揚げも普通だ。これならば自宅で食べるスーパーの30円の生麺（断面は四角）＋市販の味付け油揚げと大差ない。ただ，ミニ親子丼は昔ながらの玉子とじ状態でおいしかった。

食べているうちにも次々にお客さんが入ってきて，2階の席にも案内されている。隣の席ではおじさんが大阪弁で商売の話をしつつ「僕はしっぽくうどんにするわ」。お二人さんは，しっぽくうどんときつねとじうどんであり，なかなか味のある会話が流れてくる。

そういえば，金沢にいたころは，甘く炊いた大判の揚げのきつねうどんは存在しておらず，刻んだ揚げがのったいなりうどんだった。近所のうどん屋お多福の出前で，土曜の午後に天ぷらうどんとセットで食べていた。夏休みには泉丘高校のESSの部室に集まって勉強していたが，高橋君らと向かいのうどんやで大盛りのいなりうどんとかき氷を食べるのが日課だった。

入るときは気づかなかったが，大将がレジに座ってお会計をしていた。

写真：うさみ亭マツバや（撮影 2022.4.11）

摂州合邦辻

国立文楽劇場の四月文楽公演には，豊竹咲太夫文化功労者顕彰記念・文楽座命名150年のタイトルがついていた。第一部が義経千本桜，第二部が摂州合邦辻，第三部が嬢景清八嶋日記と契情倭荘子だ。コロナ流行が始まってからは，以前のような三部通しの観劇がしんどいので，今回は第二部を見ることに。

住吉大社，万代池や天王寺界隈が舞台である摂州合邦辻。その万代池の段は，三輪太夫，希太夫，南都太夫，津國太夫，咲寿太夫，清友が並ぶ。公演記録の録画カメラが入っていたので，皆さん気合が入っていたようだ。最後に，吉田玉也の合邦道心によって，高安次郎丸の鬼若が，万代池に投げ込まれるのだが，池から可愛らしい顔を覗かせながら幕となる。

合邦住家の段は，中：睦太夫・清馗，前：呂勢太夫・清治，切：呂太夫・清介だった。今日は睦太夫が良かったと，家人と意見が一致した。声もはっきりと大きく変な癖がない。呂勢太夫と清治のコンビはうまいのだろうが，どうも自分にはピンとこない。一番良かったのは，呂太夫と清介だった。三味線の迫力やテクニックもすごくて（三の糸が切れた？），呂太夫も清介も熱演だった。

後半の玉手御前と浅香姫の激しいバトルでのキックには思わず笑ってしまった。この段は何度か見ているはずだけれど，百万遍の大数珠のイメージはすっかり消えていた。奴入平の吉田玉勢が，玉手御前の肝を取り出すのをビビって拒む当たりも面白い。

物語では，寅年，寅の月，寅の日，寅の刻に生まれた女の生肝の血が俊徳丸の毒を消すための薬になるという設定になっている。今年は寅年である。寅の月は旧暦の正月であり，2022年では，2月4日から3月4日までだ。寅の日は，12日ごとに巡ってきて，この度では，2月6日（日），2月18日（金），3月2日（水）。寅の刻は，24時間を12で割った3番目なので，午前3時から5時となる。

写真：浪速区下寺西方寺にある合邦辻閻魔堂（Wikipediaより引用）

真木悠介

 4月10日の夜7:00のNHKニュースで，見田宗介（1937-2022）の訃報が流れた。そういえば昔よく読んだものだと本棚を探したけれど見あたらなかった。何度か本の断捨離をしたときに，寄付に回してしまったようだ。

テレビでは，「現代社会の理論 −情報化・消費化社会の現在と未来−」を見田の代表作の一つとしてあげていた。調べてみたらこれは1995年の岩波新書であり，そのころには読む気がしなくなっていたので買わなかった・・・たぶん。大学1,2年の頃に，筑摩書房の箱入りのシリーズで，真木悠介名義の「人間解放の理論のために」とか，見田宗介名義の「現代日本の心情と論理」を買って愛読していた。いずれも，1971年の出版である。

当時は，ディヴィッド・リースマンの「孤独な群衆」から始まって，社会学の本を何冊か読んでいた。大学の教養の人文科学・社会科学では，哲学，倫理学，心理学，政治学，経済学は履修したが，社会学は取らなかったので，その埋め合わせもあり，「社会学のすすめ」から始まって色々と読みあさった・・・たぶん。

社会学といえば，1972年に阪大教養部の助教授だった井上俊（1938-）の「死にがいの喪失（1973）」はまだ本棚で生き延びていた。当時は真木悠介の本の方が良いと思っていたはずなのだけれど。

Wordle（８）

Wordle（７）からの続き

Wordle支援プログラム word.py を改良して，出現しない文字を含む単語をとり除く処理を１行加えた。第3引数に出現しない文字を連結して並べたものを与える。

# usage: warp.py w.txt a.b.c def

import sys
import re

f = open(sys.argv[1], 'r')
datalist = f.readlines()
arg1 = sys.argv[2]
arg2 = sys.argv[3]

for data in datalist:
  if(re.search(r'\d* '+arg1, data) != None):
    if(re.search('['+arg2+']', data) == None):
      print(data, end="")

f.close()

Wordle（７）

Wordle（６）からの続き

wiki-5.txtは，WIkipediaデータセットのWikiText-103 word から5文字単語を取り出したものだった。これに対して，n-gramの出現頻度を，1-gram， 2-gram， 3-gramに対して求めるpythonプログラム ngrm.py を作る。1-gramはアルファベット1文字，2-gramは連続するアルファベット2文字，3-gramは連続するアルファベット3文字を表している。

ngrm.pyの出力結果は , を区切り文字として，出現回数とn-gramパターンの組なので，sort で区切り文字を , 数値として解釈した逆順をリダイレクトして適当なファイルに格納すると，1-gramから3-gramまで混ぜた出現頻度順のファイルが出来上がる。

# usage: ngrm.py wiki-5.txt | sort -n -r -t , >! alphabet 

import sys

a=[0 for i in range(26)]
b=[[0 for i in range(26)] for j in range(26)]
c=[[[0 for i in range(26)] for j in range(26)] for k in range(26)]

def n_gram(target, n):
  return [ target[idx:idx + n] for idx in range(len(target) -n + 1) ]

f = open(sys.argv[1], 'r')
datalist = f.readlines()

for data in datalist:
  for l in range(5):
    x = data[l]
    i = ord(x) - 97
    a[i] = a[i] + 1
    if l in range(4):
      y = data[l+1]
      j = ord(y) - 97
      b[i][j] = b[i][j] + 1
      if l in range(3):
        z = data[l+2]
        k = ord(z) -97
        c[i][j][k] = c[i][j][k] + 1

for i in range(26):
  if a[i] > 49:
    print(a[i],",",chr(i + 97))

for l in range(26*26):
  i = l // 26
  j = l % 26
  if b[i][j] > 49:
    print(b[i][j],",",chr(i+97)+chr(j+97))

for l in range(26*26*26):
  i = l // (26*26)
  j = (l // 26) % 26
  k = l % 26
  if c[i][j][k] > 49:
    print(c[i][j][k],",",chr(i+97)+chr(j+97)+chr(k+97))

f.close()

wiki-5.txtにおけるベスト50は次のとおりであった（数字は出現数）。

5730707 , e

4092848 , r

3965317 , a

3759319 , t

3424032 , s

3108801 , o

3068296 , i

2603632 , h

2600893 , l

2398537 , n

1734820 , d

1574948 , c

1501211 , u

1291809 , w

1195832 , m

1086248 , g

1076779 , er

964115 , f

962898 , th

882193 , b

801081 , p

746024 , te

729576 , y

674447 , he

651813 , ou

636278 , st

604734 , ar

563176 , hi

559096 , re

556926 , v

556020 , or

533535 , in

521374 , k

521197 , es

519713 , ch

518219 , ir

514744 , al

514410 , wh

493762 , the

473288 , ea

427832 , le

417701 , ic

412433 , whi

412193 , an

409812 , il

384013 , la

372098 , at

369592 , se

364056 , ed

362684 , ter

Wordle（６）

 Wordle（５）からの続き  　

Wordleを解くための支援プログラムを次のようなものとする。「5文字単語のデータベースを用意して，ワイルドカードを含む与えられた候補にマッチする単語を全て選び出して表示する」。perlでプログラムを組みたいところだったが，pythonの練習をすることにした。

5文字単語データベースとしては，前回示した WikiText-103 word Level からwords.shで処理したものを用いる。その後は，久々の python プログラミングなので，よちよち歩きながら試行錯誤する。

(1) コマンドラインの引数は，sysモジュールをimportすれば，sys.argv[n]でアクセスできる。sys.argv[0]はpythonスクリプト名，sys.argv[1] 以下で複数の引数を受け渡す。

(2) ファイルのオープンとクローズの方法は記載の通り。第1引数は先に求めてある5文字単語データファイル名である。

(3) テキストファイルを1行づつ分割して読み込むには，datalist=f.readlines() などとするとその結果のリストが 得られる。

(4) perlのようにdatalistから1行づつ取り出して処理をするのが，for data in datalist: である。

(5) 正規表現によるマッチングには，re モジュールをimportしたうえで，re.search(r'正規表現' , 原データ) という形式による。これが Noneでなければマッチしていることになる。

(6) 第2引数で入力する5文字のワイルドカードをピリオド（. ）にしておけばそのまま正規表現として使えてるので便利である。

(7) データファイルは，wc -l で生成しているので，頭に頻度を表した数桁の数字が空白を区切りに付加されているため，'\d* '+arg としている。文字列の結合には，プラス（+）演算子を用いる。

(8) print文で余分の改行を削除するには，end="" をつければよい。

# usage: word.py w.txt a.b.c
import sys
import re

f = open(sys.argv[1], 'r')
datalist = f.readlines()
arg = sys.argv[2]

for data in datalist:
# print(arg, data, re.search(r'\d* '+arg, data))
  if(re.search(r'\d* '+arg, data) != None):
    print(data, end="")

f.close()

Wordle（５）

Wordle（４）からの続き

前回の復習：英文のテキストデータセットと単語の文字数 n が与えられたとする。文字数がnの単語を1行とした一時ファイルを作る。大文字は小文字に変換している。次に，wc -l によって一時ファイルの行数＝n文字単語数を出力したものを並べたファイルを作る。また，上記の wc -l の前にsort | uniq フィルターをかけてユニークな n文字単語数についても同様のファイルを作る。これが，前回のwords.shプログラムのアルゴリズムだった。

具体例として，最初は手元にある数冊の英書の数千ページ（80MB）のpdfファイルからユニークにソートされた5文字単語（3万4千語）を抽出したものでテストした。もう少し大きなデータがないかと，オープンコーパスを色々探してみたが，なかなか適当なものがない。

最終的に，Wikipediaの記事を使ったデータセットにたどり着いた。多言語のWIki-40Bが前処理済で良さげだったが，TensorFlowを使えという指示が面倒であり，英語データで340万ページ2GBは大きすぎる。WikiText-103 Word Levelの方は，圧縮後181MBでそのままダウンロードできそうなので，こちらを使うことにした。

WikiText-103 Word Levelの解凍後のテキストデータは539MBあった。これをwords.shにかけると，1-20文字の単語（うちユニークな語）が約8300万語（約21万語），このうち5文字の単語（うちユニークな語）が約960万語（約2万4千語）であった。

WikiText-103 Word Level における，n文字単語（うちユニークな語）の出現頻度を下図に示す（オレンジが総語数の出現頻度，ブルーがユニーク語数の出現頻度）。n文字単語の出現ピークはn=3にあって約20%，ユニークなn文字単語の出現ピークはn=7にあって，約16%となった。5文字単語の場合は，ともに約11%となっている。

図：WikiText-103 Word Levelにおける単語出現頻度（横軸は文字数）

（春休み 8）

「月明や沖にかゝれるコレラ船」　（日野草城　1901-1956） 

（春休み 7）

  「花水木咲き新しき街生まる」　（小宮和子） 

（春休み 6）

「 鹿になる考えることのなくなる」　（阿部完市　1928-2009）

（春休み 5）

 「 みつまめをギリシャの神は知らざりき」　（橋本夢道　1903-1974） 

（春休み 4）

「 さきみちてさくらあをざめゐたるかな」　（野澤節子　1920-1995） 

（春休み 3）

  「満開のふれてつめたき桜の木」　（鈴木六林男　1919-2004）

（春休み 2）

 「手をつけて海のつめたき桜かな」　（岸本尚毅　1961-）

（春休み 1）

  「つなぐ手にをさなの湿り夕ざくら」　（千代田葛彦　1917-2003）

Wordle（４）

Wordle（３）からの続き

Wordleをやっていると，そもそも英語の文章で出てくる単語のうち5文字である確率はどんなものか，また，最もよく出現するのは何文字の単語なのか，などなどが気になるようになった。

これを調べるために，与えられたテキストファイルやpdfファイルから単語を切り出して，その文字数の分布を調べるためのシェルスクリプトを作ってみた。

case $3 in
"txt")
echo "txt";
;;
"pdf")
echo "pdf";
pdftotext $2.pdf;
;;
*)
echo "undefined";
;;
esac

for ((i=1 ; i<=$1 ; i++))
do
perl -nse 'while (/\b[a-z]{$num}\b/ig) {print "$&\n";}' -- -num=$i $2.txt | tr A-Z a-z > tmp.txt
cat tmp.txt | wc -l >> $2-\s1.txt;
cat tmp.txt | sort | uniq | wc -l >> $2-\s2.txt;
rm tmp.txt
done

昨日のスクリプトを少しだけ修正すればよかったが，ポイントは，シェルスクリプト中の反復の記述法である。繰り返し変数は$をつけて，perlのワンライナーに受け渡すことができた。あとはこれを使って実験してみれば良いのだが，それはまた次回のお楽しみ。

Wordle（３）

Wordle（２） からの続き

Wordleは英語の勉強になる。簡単な5文字の英単語でも知らないものがたくさんあって，自分の語彙力はやはり10歳並みだということが確かめられる。そのため，辞書の助けがないと5-6回で答えに辿り着くことは出来ない。

巷には，Wordle Word Finderなどというツールも登場しているが，これでは英語学習の役には立たない。それでも，5文字の英単語のリストを収集したいという目的のために，次のようなツールを作った。与えられたテキストファイルあるいはpdfやhtmlファイルをテキストに変換したものから，perlのワンライナーで n 文字の単語を切り出すものだ。

pandocは汎用の文書型変換ツールだけれど，pdfやhtmlからテキストを取り出すという想定がない。あくまでも整形された文書変換ツールだからだ。そこでテキストファイルの取り出しには，pdftotextやtextutilなどのコマンドに任せることにした。また，シェルスクリプトの引数はperlのワンライナーにそのまま送れないので， -- -perl_var=$shell_varとして，perlの中で$perl_varによって引用することになる。

#! /bin/zsh
# usage: word.sh 5 sample pdf (tst, txt, pdf, html)
# output sample-5.txt

case $3 in
"tst")
pwd;
ls -al $1.txt;
ls -al $2.*;
;;
"txt")
echo "txt";
perl -nse 'while (/\b[a-z]{$num}\b/ig) {print "$&\n";}' -- -num=$1 $2.txt | tr A-Z a-z > $2-$1.txt;
;;
"pdf")
echo "pdf";
pdftotext $2.pdf;
perl -nse 'while (/\b[a-z]{$num}\b/ig) {print "$&\n";}' -- -num=$1 $2.txt | tr A-Z a-z > $2-$1.txt;
;;
"html")
echo "html";
textutil -convert txt $2.html;
perl -nse 'while (/\b[a-z]{$num}\b/ig) {print "$&\n";}' -- -num=$1 $2.txt | tr A-Z a-z > $2-$1.txt;
;;
*)
echo "undefined"
;;
esac

ロフテッド軌道（３）

 ロフテッド軌道（２）からの続き

初速度を第１宇宙速度$v=\sqrt{\frac{GM}{R}}$に固定した場合，Mathematicaによるプロットを試みる。このとき，$u(\varphi)=A(\theta) \cos \varphi + B(\theta) \sin \varphi +\frac{GM}{h^2}$，$A(\theta)= - \frac{1}{R \tan^2 \theta }$，$B(\theta) = -\frac{1}{R \tan \theta}$，$\frac{GM}{h^2}=\frac{1}{R \sin^2 \theta}$ であり，Mathematicaのコードは以下のとおり。

In[1]:= {G, R, M} = {6.67*10^-11, 6.37*10^6, 5.97*10^24}
Out[1]= {6.67*10^-11, 6.37*10^6, 5.97*10^24}
In[2]:= {g, v} = {G M /R^2, Sqrt[G M/R]}
Out[2]= {9.81344, 7906.43}

In[3]:= A[t_] := 1/R*(1 - 1/Sin[t]^2)
In[4]:= B[t_] := -1/(R*Tan[t])
In[5]:= GM[t_] := 1/(R Sin[t]^2)
In[6]:= r[s_, t_] := 1/(A[t] Cos[s] + B[t] Sin[s] + GM[t])
In[7]:=
Table[f1[i] = Plot[{r[s, i]/R, 1}, {s, 0, Pi},
GridLines -> Automatic, PlotRange -> {0, 2.0}];, {i, 0.1, 1.5, 0.2}];
In[8]:= Show[Table[f1[i], {i, 0.1, 1.5, 0.2}]]
In[9]:=
Table[f2[i] = PolarPlot[{r[s, i]/R, 1}, {s, 0, Pi},
GridLines -> Automatic, PlotRange -> {0, 2.0}];, {i, 0.05, 1.5, 0.2}];
In[10]:= Show[Table[f2[i], {i, 0.05, 1.5, 0.2}]]

In[11]:= f[t_] :=
NIntegrate[ 1/Sqrt[G M R Sin[t]^2] *1/(A[t] Cos[s] + B[t] Sin[s] + GM[t])^2, {s, 0.001, Pi/30}]

In[12]:= Plot[f[t], {t, 0, Pi/8}]

図１：ロフテッド軌道（横軸$\varphi$ラジアン，縦軸$r/R$）

図２：ロフテッド軌道（上記と同じものを極座標表示）

この近似のもとでは，6500kmの高度に達する弾道ミサイルはほぼ地球を半周するところまで到達できるので，北朝鮮から米国本土全体が射程に入ることになる。鉛直上方から60度の角度で射出すれば，地球を1/4周するので1万kmまで到達できる。

実際には，空気抵抗があることや，初速度はロケットエンジンによる一定時間の加速によって獲得されることなどから，もう少し真面目な計算をする必要がある。

［１］ロケットの運動と人工衛星の打ち上げ（冨田信之）

ロフテッド軌道（２）

ロフテッド軌道（１）からの続き

重力場中の質点に対する２次元極座標系$(r(t), \varphi(t))$でのニュートンの運動方程式は次のとおり。

$\quad \quad m \bigl( \ddot{r} -r \dot{\varphi}^2 \bigr) = F_r = \frac{G M m}{r^2}$

$\quad \quad  m \frac{d}{dt} \bigl( r^2 \dot{\varphi} \bigr)  = F_\varphi = 0$

第2式は角運動量保存則に対応しており，$r^2 \dot{\varphi} = h$ (一定) が成り立つ。

これを第1式に代入して，$\dot{\varphi}$ を消去すると，$\ddot{r} -\frac{h^2}{r^3}  =\frac{G M}{r^2}$となる。そこで，$u(\varphi) =\frac{1}{r(\varphi)}$と置き，$\dot{\varphi} =h u^2$に注意して，運動方程式を，軌道の形を定める$\varphi$に関する微分方程式に書き換える。

$\dot{r} = \frac{d}{d\varphi} \bigl( \frac{1}{u} \bigr ) \dot{\varphi}= -\frac{1}{u^2} \frac{du}{d\varphi} \ h u^2 =  -h \frac{du}{d \varphi}$

$\ddot{r} = \frac{d}{d\varphi} \bigl( -h \frac{du}{d \varphi} \bigr ) \dot{\varphi}= - h \frac{d^2 u}{d\varphi^2}  \ h u^2 =  -h^2 u^2 \frac{d^2 u}{d \varphi ^2}$

$\dot{\varphi}$を消去した運動方程式に代入すると， $-h^2 u^2 \frac{d^2 u}{d \varphi ^2} -h^2 u^3 = -GM u^2$，つまり，$\dfrac{d^2 u}{d \varphi^2} + u = \dfrac{GM}{h^2}$であり，一般解は，$u = A \cos \varphi + B \sin \varphi + \frac{GM}{h^2}$ である。

初期条件から$A,\ B$を定めるには，$\dot{r} = -h \frac{du}{d \varphi}$の表式が必要となる。一般解を右辺に代入すると，$\dot{r} = h (A \sin \varphi - B \cos \varphi )$

半径$R$の地球の中心Oを原点にとった座標系において，$(x,y)=(R,0); (r,\varphi)=(R,0)$から鉛直上方（$x$軸正方向）となす角度$\theta$の方向に，初速度$v$で質量$m$の質点を投射する。このとき $r(0)=R$，$\dot{r}(0)=v_r=v \cos \theta$，$h = r(0)^2 \dot{\varphi}(0) = R v_\theta = R v \sin \theta$である。

これらから，$\frac{1}{R} = A+ \frac{g}{(v \sin \theta ) ^2}$，$v \cos \theta = - B R v \sin \theta$であり，$A=\frac{1}{R} - \frac{g}{(v \sin \theta ) ^2}$，$B = -\frac{1}{R \tan \theta}$となる。

なお，到達時間のオーダーは，$\frac{R \Delta \varphi}{v \sin \theta}= 800 \frac{ \Delta \varphi}{\sin \theta}$[s]となる。

図：ロフテッド軌道の例