phys-aのブログ

物理に関することちょこちょこまとめられたらな、と

ファインマンパラメータ積分

ファインマンパラメータ積分は (ファインマン・パラメータ積分 - Wikipedia, Feynman parametrization - Wikipedia) 何度か導出した気がするけど毎回忘れて時間を浪費してしまう、というわけでしっくり来たものを一旦見えるところに書き記しておこうというのが今回の内容です。例のごとく上にWikipediaのリンクを貼っているので、そこで納得できればここに書かれてることを読む必要ありません。英語版Wikipediaの方は導出をちゃんと書いており、ここではその行間を埋めて書いているだけです。

 

ここで考える  n 個のパラメータ  a_1, \cdots , a_n に対するファンマンパラメータ積分

\begin{align} \frac{1}{a_1 \cdots a_n} = \int ^1_0 \mathrm{d} x_1 \cdots \mathrm{d} x_n \frac{ \Gamma(n) \delta(1- x_1 - \cdots - x_n)}{ [ a_1 x_1 + \cdots + a_n x_n]^n} \end{align}

である。

 

まずは基本的な話として冪函数はガンマ函数を使うことで

\begin{align} \frac{1}{a^s} = \frac{1}{\Gamma(s)} \int^\infty _0 \mathrm{d} t \, t^{s-1} e^{-a t} \end{align}

のような積分の形で表わすことが出来る(若干の参考: 演算子のゼータ函数 - phys-aのブログ)。これは右辺の式で  t \to t/a のように変数変換を行えば元に戻ることから確認できる。この式に  s=1 を代入し、  n 個の積を考えれば

\begin{align} \frac{1}{a_1 a_2\cdots a_n} = \int^\infty _0 \mathrm{d} t_1 \, e^{-a_1 t_1} \, \int^\infty _0 \mathrm{d} t_2 \, e^{-a_2 t_2} \cdots \int^\infty _0 \mathrm{d} t_n \, e^{-a_n t_n}\end{align} 

となる。これが出発点である。

 ここで独立だった  n 個の無限積分 n 次元空間でのひとまとまりの積分と見なし*1、以下のような変数変換を施す。

\begin{align} \begin{cases} t_i = t x_i & ( i = 1, 2, \dots, n-1) \\ t_n = t(1 - x_1 + x_2 + \cdots + x_{n-1}) \end{cases}. \end{align}

これで積分変数が  (t_1, \dots, t_n) \to (x_1, \dots, x_{n-1}, t=x_n) の組に置き換えられることになる。またそれに伴い積分区間 x_i \in [0,1] \, \text{for} \, i=1, 2, \dots, n-1 かつ  t \in [0, \infty) となる。

変数変換を施した積分

\begin{align} \frac{1}{a_1 a_2\cdots a_n} = \int^1 _0 \mathrm{d} x_1 \cdots \mathrm{d} x_{n-1} \int ^\infty _0 \mathrm{d} t \, \mathrm{det} \! \left( \frac{ \partial t_j }{\partial x_k} \right) e^{-t \left( a_1 x_1 + \cdots + a_n (1- x_1 + \cdots + x_{n-1}) \right)} \end{align}

と表わされる。 x_i \, (i=1,2, \ldots, n-1) はそれぞれ同じ積分範囲なのでひとつのみ書いてあとは省略した。積分測度は

\begin{align} \mathrm{det} \begin{pmatrix} \frac{ \partial t_1 }{ \partial x_1} & \frac{ \partial t_1 }{\partial x_2} & \cdots & \frac{ \partial t_1 }{\partial x_n} \\ \frac{ \partial t_2 }{\partial x_1}& \frac{ \partial t_2 }{\partial x_2} & \cdots & \frac{ \partial t_2 }{\partial x_n} \\ \vdots& \vdots& \ddots & \vdots \\ \frac{ \partial t_n }{\partial x_1} & \frac{ \partial t_n }{\partial x_2} & \cdots & \frac{ \partial t_n }{\partial x_n} \end{pmatrix} = & \mathrm{det}  \begin{pmatrix} t & 0 & \cdots & x_1 \\ 0 & t & \cdots & x_2 \\ \vdots & \vdots & \ddots & \vdots \\ -t & -t & \cdots & 1 - \sum_{j=1}^{n-1} x_j \end{pmatrix} \\ = & \mathrm{det}  \begin{pmatrix} t & 0 & \cdots & x_1 \\ 0 & t & \cdots & x_2 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{pmatrix} \\ = & \mathrm{det}  \begin{pmatrix} t & 0 & \cdots & 0 \\ 0 & t & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{pmatrix} \\ = & t^{n-1} \end{align}

と求まるので

\begin{align} \frac{1}{a_1 a_2\cdots a_n} &= \int^1 _0 \mathrm{d} x_1 \cdots \mathrm{d} x_{n-1} \int ^\infty _0 \mathrm{d} t \, t^{n-1} e^{-t \left( a_1 x_1 + \cdots + a_n (1- x_1 + \cdots + x_{n-1}) \right)} \\ =& \int^1 _0 \mathrm{d} x_1 \cdots \mathrm{d} x_{n-1} \frac{\Gamma(n)}{ \left[ a_1 x_1 + \cdots + a_n (1- x_1 + \cdots + x_{n-1}) \right]^n} \end{align}

が得られる。途中で変数変換を行って  t 積分を行ってガンマ函数にしている。積分は一つ実行できてしまったので  n-1 個になっている。この形でも充分だが、式の対称性を気にするならデルタ函数を使って積分を一つ増やすことで

\begin{align} \frac{1}{a_1 a_2\cdots a_n} = \int^1 _0 \mathrm{d} x_1 \cdots \mathrm{d} x_{n-1} \int_0^1 \mathrm{d} x_n \frac{\Gamma(n) \delta(1 - x_1 - \cdots - x_n) }{ \left[ a_1 x_1 + \cdots + a_n x_n \right]^n} \end{align}

 が得られる。これが冒頭で記していたファンマンパラメータ積分の式である。

 

もうすこし一般的な形にしたければ最初の冪函数をガンマ函数で表わす部分で  s=1 としないでそのまま計算を進めればよいことが分かる。

*1:ガウス積分の計算と同じお気持ち。

周期的な強制振動(減衰項無し)の補足の補足

前回の周期的な強制振動の計算で適当な時刻で発散しなかったというのがあったが、改めてやったら一応求まったので載せておく。時刻  s 、パラメータ  \alpha において

\begin{equation} \lim_{\alpha\to1} \frac{\sin \alpha s}{1-\alpha^2} \end{equation}

を考える。

この極限は一般には発散するが  s = l\pi \quad (l=0,1,2,3,\dots) *1のときは

\begin{align} \frac{\sin \alpha s}{1-\alpha^2} &= -\frac{-1}{1+\alpha} \frac{\sin (1 - \alpha l)\pi}{1-\alpha} \\ &= -\frac{(-1)^l  }{1+\alpha} \frac{\sin (1 - \alpha)\pi l}{1-\alpha} \\ &= -\frac{(-1)^l \pi l }{1+\alpha} \frac{\sin (1 - \alpha)\pi l}{(1-\alpha)\pi l } \\ &= -\frac{(-1)^l \pi l  }{2} \qquad (\alpha \to 1 ) \end{align}

が得られる。前回の記法に合わせるには  (-1)^l = \cos \pi l とすれば良い。

 

意味があるかわからないが、とりあえず正弦函数を無限積でバラして適当な時刻で発散がないことは確認できた一方でそれが解析的に分かるか謎だったので、ちゃんと関係が与えられて一安心。高校数学の練習問題だった。

 

グラフで見た時に面白いなと思ったのはうえで得られるのは、元々 \alpha=1 で計算して得られる(振幅が単調増加する)グラフの極大値というわけでも別になかったという点。

*1:普通の時刻として見てるから0以上で考えている。

周期的な強制振動(減衰項無し)の補足

前回(周期的な強制振動(減衰項無し) - phys-aのブログ)は減衰項なしで強制振動を考え、固有周期と強制振動を与える周期によって構成されるパラメータ  \alpha の値で場合分けをして解く、ということを述べた。その時  \alpha \neq 1 において  \alpha \to 1 を考えると発散するとも述べたが、それは一部誤っていることに気づいたので、そしてそれが何なのか自分でもよく分からなくて面白いと思ったのでここに書いておこうと思った。

 

まずは  \alpha \neq 1 において得られる運動方程式の解を述べておくと

\begin{equation} x = A \sin(s+\delta) + \frac{X}{1-\alpha^2} \sin\alpha s \end{equation}

であった。1項目の斉次での解は単なる振動を表していて、もちろんここで注目したいのは第2項。 そもそも  \alpha \neq 1 が前提として解いているのだが、一応その極限も見てみたいというのが人というもので、ぱっと見発散しそうである、というか実際基本的には発散する。ただし一部発散しないときがあるというのに気づいた。というのもサイン函数は周期関数で  \alpha s = \pi n, \quad (n \in \mathbb{Z}) でゼロになるからで、その時は分母がゼロになるからと言って一概に発散するとは限らない。それを調べるために思い出されるはサイン函数の無限乗積表示

\begin{equation} \sin x = x \prod_{n=1}^{\infty} \left( 1 - \left(\frac{x}{n \pi}\right)^2 \right) \end{equation}

である。運動方程式の第2項をこの無限乗積で書いてやると

\begin{equation} \frac{\sin \alpha s}{1-\alpha^2} = \frac{\alpha s}{1-\alpha^2} \prod_{n=1}^{\infty} \left( 1 - \left(\frac{\alpha s}{n \pi}\right)^2 \right) \end{equation}

である。これを見れば時刻が  s = l \pi, ( l = 0,1,2,3,\dots) の時は

\begin{equation} \frac{\sin \alpha  l \pi }{1-\alpha^2} = \frac{\alpha s}{1-\alpha^2} \prod_{n=1}^{\infty} \left( 1 - \left(\frac{ l }{n}\right)^2 \alpha^2 \right) \end{equation}

となり  \alpha \to 1 で有限な値になる……だけです。これに何か意味があるのかよくわかっていません。MathematicaでPlotとしてたらなんか楽しかった程度の内容です。値も数値計算でしかわからないんですかね……。一応ちゃんと  \alpha =1 で得られる  -\frac{l\pi}{2} \cos \alpha l \pi \alpha \to1 で一致してるようですが。

 

gnuplotで適当にプロットしたのを以下に貼っておきます。ひとつ目は発散するパターン。ふたつ目が時刻が適当な値の時の発散しないパターン。比較対象として  \alpha =1 での定数の直線も併せてプロットしてる。  \alpha=1 で交わっている(一致している)ように見る。

f:id:phys-a:20181117151501p:plain

 

f:id:phys-a:20181117151041p:plain

補足しておくと、この第2項の導出の仕方によっては  \alpha \to 1 の極限で微分の形になって  \alpha = 1 での形、時刻と共に振幅が増加する解が得られる、というのもあります。適当に力学の本やネットを漁れば出てくると思います。

 

周期的な強制振動(減衰項無し)

強制振動は単振動の方程式  m \ddot{x} = - k x に外力に相当する項を加えた方程式で表わされる。ここでは周期的な外力として  f\, \sin \Omega t を加えることにする。いつも通り周波数を  \omega= \sqrt{\frac{k}{m}} と置いて、運動方程式を簡略化するために  s=\omega t のように時間を無次元化し、それに合わせて外力の係数を  X=\frac{f}{m\omega^2} 、周波数を  \alpha =\frac{\Omega}{\omega}, (\omega,\Omega \gt 0) と置く。考える運動方程式

\begin{equation} \ddot{x} + x = X \sin \alpha s \end{equation}

と書ける。ここでのドットは  s での微分を表わすことにする。

 

ここでは愚直に  \alpha が充分に小さい(  \alpha \ll 1 )と仮定して順次次数を上げていく戦略を取ってみる。これは強制振動の周期  \Omega が固有周 \omega より充分に短くするということ、つまり小刻みに外力を与えることを意味する。ブランコが揺れてる時にジタバタしてるイメージすれば良いと思う。

 

まずは  \alpha の1次までを見る。つまり  \ddot{x} + x = X \alpha s を解けばよく、これはさらに2階微分すると  \ddddot{x} + \ddot{x} = 0 となることからカーネルを考えれば定数  A, \delta に対して  \ddot{x} = A \sin(s+\delta) が得られる。単純に積分して元の微分方程式に代入してやると

\begin{equation} x = A \sin(s+\delta) + X\alpha s \end{equation}

が解となる。

次に  \alpha の3次まで見る*1 \ddot{x} + x = X \left( \alpha s - \frac{1}{3!}(\alpha s)^3\right) を解く事になる。これをさらに4階微分すると  x^{(6)} + x^{(4)} = 0 *2が得られるので4階微分 x^{(4)} = A \sin(s+\delta) となり、再び積分して元の微分方程式に代入してやると

\begin{equation} x = A \sin(s+\delta) + X\left( \left( 1+\alpha^2\right) \alpha s - \frac{1}{3!}(\alpha s)^3 \right)  \end{equation}

が得られる。

同様に5次、7次と見ていけば次数を上げるほどに良い近似式が得られる。一方で  \alpha が充分小さいと仮定したままでは強制振動の周期性が見えないままである。

 

ここでやってることはいわゆる摂動展開と言われるものである。摂動パラメータ  \alpha が充分小さいものとしてその冪級数で展開してその近傍での挙動を見る、という戦略である。線形な微分方程式なら一般的な解法が存在するので厳密解が求まるが、今回のような非線形微分方程式の場合は一般的な解法が存在しない。したがって順次次数を上げながらパラメータで展開していくことで予測できる範囲を広げていこう、ということである。

ただしパラメータが大きくなるとこの戦略は破綻してしまうので、この戦略ではパラメータが十分小さいようなローカルな情報は分かるが解の周期性のようなグローバルな情報が分からない。

 

このままでは埒が明かないのでとりあえずこれまでの計算結果に照らしあわせて解を以下のように仮定してみる。まずはカーネルの部分は3次まで考えても特に変わりがなかったのでそのままにして、後ろの項は周期は変わらないものとして  \sin\alpha s に比例すると仮定する。またその係数は s を含まないものとする。つまり  x = A \sin(s+\delta) + \beta X \sin\alpha s と仮定してみる。これを微分方程式に代入すると

\begin{equation} \beta = \frac{1}{1-\alpha^2} \end{equation}

が得られる。つまり

\begin{equation} x = A \sin(s+\delta) + \frac{X}{1-\alpha^2} \sin\alpha s \end{equation}

が解として得られる。逆に  \alpha=0 近傍で展開すれば上記の近似式と一致することが確認できる。またこの解の形から  \alpha=1特異点になっているが分かる。物理的には固有周期と同じ周期で外力を与えることを意味する。加えて摂動展開では  0 \lt \alpha \lt 1 の範囲のみ有効だったがこの式では  \alpha \gt 1 でも有効である。

 

先述の通りこの解には1が代入できないので元の微分方程式に戻って代入してやると

\begin{equation} \ddot{x}+x=X \sin s \end{equation}

である。この微分方程式でのカーネルは上で既に書いてあるように  A\sin (s+\delta) で表わすことができる。ここで注目するのは非斉次項とこの斉次式での解とは独立になってないことである。なんとなく係数の A をいじれば解けてしまいそうに思えるのでここで改めてこの解を  A \sin s + B \cos s のようにバラしてやって  A s の式と見なおして改めて微分方程式に代入してやると

\begin{equation} \ddot{x}+x = \ddot{A} \sin s + 2\dot{A} \cos s \neq X \sin s \end{equation}

となり、これは失敗。次に  B s の式と見なおして計算すると今回は成功してカーネルの係数の数を保つように書いてやると  B = B^\prime - \frac{s}{2}X が得られる。改めてカーネルに属する方を1つの項にまとめた表記だと

\begin{equation} x = A \sin(s+\delta) - \frac{s X}{2} \cos s \end{equation}

 \alpha=1 での解ということになる。この解法は定数変化法と呼ばれるものである。

 

まとめると減衰項無しの周期的な強制振動では固有周期と強制振動の周期が一致しない時はその差でうなりが生じ、一致するときは振幅が時間に比例して増加する。

*1:2次を飛ばした理由は強制振動がサイン函数だから。

*2:冪の (n) はn階微分の略記とする。

第2種変形Bessel函数: 原点近傍

前回は第2種変形Bessel函数  K_\nu(x)積分形っぽいのを見てみた。今回は原点付近でどのように展開できるのかを見てみる。参考にしたのはみんながお世話になってるだろう岩波数学公式集のIII。

特殊函数 (岩波 数学公式 3)

特殊函数 (岩波 数学公式 3)

 

 

そのp172にある関係式

\begin{equation} K_{n+\frac{1}{2}}(x) = \sqrt{\frac{\pi}{2x}} e^{-x} \sum_{r=0}^n \frac{(n+r)!}{r! (n-r)!} \frac{1}{(2x)^r}, \quad (n\in\{0,1,2,3,\dots\} )\end{equation}

を元に  x=0 近傍を調べる。両辺に  x^{n+\frac{1}{2}} を乗じてやると

\begin{equation} x^{n+\frac{1}{2}} K_{n+\frac{1}{2}}(x) = \sqrt{\pi} e^{-x} \sum_{r=0}^n \frac{\Gamma(n+r+1)}{\Gamma(r+1) \Gamma(n-r+1)} \frac{x^{n-r}}{(2)^r} \end{equation}

となる。ここで階乗はガンマ函数に書き換えておいた。ここで  x \ll 1 とおくと

\begin{equation} x^{n+\frac{1}{2}} K_{n+\frac{1}{2}}(x) = \frac{\sqrt{\pi}}{(2)^n} \frac{\Gamma(2n+1)}{\Gamma(n+1)} + O(x) \end{equation}

 のように一番ドミナントに効く定数項が得られる。次に  n+\frac{1}{2}\to\nu にすることで

\begin{equation} K_{\nu}(x) = \frac{\sqrt{\pi}}{(2)^\nu} \frac{\Gamma(2\nu)}{\Gamma(\nu+\frac{1}{2})}x^{-\nu} + O\left(x^{-\nu+1}\right), \quad(x \ll 1) \end{equation}

が得られる。ガンマ函数に書き換えてるので0と正の整数とに限定されない式になっている。ということで  x=0 近傍では一番効く発散が  x^{-\nu} に比例していることが確認できた。加えてその係数も確認できた。

 

十中八九もっと自明な調べ方があると思うけど自分はこうしたってのを挙げている。

第2種変形Bessel函数

特殊函数の1つにBessel函数というのがある(ベッセル関数 - Wikipedia)。   J_\alpha (x), Y_\alpha (x), I_\alpha (x), K_\alpha (x) のようにいくつか種類があって物理の計算の中で時折顔を出す。そのいくつかある中で第2種変形Bessel函数に関係する話として以下の積分に触れてみる。

\begin{equation} A(\nu,x) = \int_0^\infty \mathrm{d}t\, \left( 2 t \right)^{-\nu-1}  \exp\left[-x \left(\frac{1}{4t} + t \right) \right] \qquad \mbox{for } \mathrm{Re} \, x > 0  \end{equation}

この積分 x について微分すると  \frac{1}{2}\left(A(\nu+1,x)+A(\nu-1,x)\right) が得られる。第2種変形Bessel函数  K_\alpha(x)微分すると同様の漸化式を満たす。また第2種Bessel函数K_\alpha(x) = K_{-\alpha} (x) を満たすのに対し上記の積分でも積分変数を  t \to \frac{1}{4 s} と変換すれば満たされることが分かる。

というわけでこの積分 x が正の実数の範囲内で第2種変形Bessel函数を与える積分らしい(ただしちゃんと確認できていない。)。

線形微分方程式についての補足

パラメータ  t を含む変数 x の2階線形微分方程式にFourier展開した解を実際に代入する。例えば

\begin{equation} \frac{\mathrm{d}^2 x}{\mathrm{d} t^2} + 2\mu \frac{\mathrm{d} x}{\mathrm{d} t} + \omega^2 x = 0 \end{equation}

\begin{equation} x = \int_{-\infty}^{\infty} \mathrm{d}\lambda\, \tilde{x}(\lambda) e^{-i \lambda t} \end{equation}

を代入してみる。ここで  \lambda は時間を無次元化するための量であり、  \tilde{x}(\lambda) 係数である。すると

\begin{align} \int_{-\infty}^{\infty} \mathrm{d}\lambda\, \tilde{x}(\lambda )\left( -\lambda^2 + i2\mu \lambda + \omega^2 \right)e^{-i \lambda t} = 0 \end{align}

が得られる。この方程式を満たす係数  \tilde{x}(\lambda) \lambda を、任意の  \mu, \omega に対して考える。すなわち

\begin{equation} \tilde{x}(\lambda)\left( -\lambda^2 + i2\mu \lambda + \omega^2 \right)=0 \end{equation}

を考えればよくて、常に  \tilde{x}(\lambda)=0 が成立するのは自明解。一方で後ろの2次代数方程式の解  \lambda_\pm = -i\mu\pm\sqrt{\omega^2-\mu^2} においては  \tilde{x}(\lambda) がゼロである必要が無いことが分かる。まとめると係数は任意の定数  A,B に対して

\begin{equation} \tilde{x}(\lambda) = A\, \delta(\lambda-\lambda_+) + B\, \delta(\lambda-\lambda_-) \end{equation}

を満たせば良いことが分かる。ここで  \delta(\lambda) はデルタ函数である。したがって見慣れた2つの解が求まることになる。ここでのポイントは適当な階数の線形微分方程式である限り係数の形は

\begin{equation} \tilde{x}(\lambda) = \sum_k A_k \, \delta(\lambda-\lambda_k) \end{equation}

の形になってしまうことである。これならわざわざFourier展開するまでもなく  e^{\lambda t} を1つ代入してやればその代数方程式の解がそのまま微分方程式の解になってしまうことが見て取れる。