ツイスパラメーター統計

加速器のビームは，多数の粒子が位相空間に分布した集合です．個々の粒子の軌道だけを見ても，ビーム全体の広がりや傾き，輸送中の変形は把握しにくくなります．そのため，加速器物理では粒子分布を二次モーメントで要約し，ツイスパラメーターとエミッタンスで表すことがよくあります．

このページでは，一つの横方向位相空間 \((x,\,x^\prime)\) を例に，ツイスパラメーター，二次元正規分布，共分散行列，RMS エミッタンスの関係をまとめます．目的は，「粒子分布から何を計算すればツイスパラメーターが得られるのか」を明確にすることです．

粒子分布の二次モーメントから，共分散行列，RMS エミッタンス，ツイスパラメーターを求める流れです．

粒子の運動は，本来は六次元位相空間 \((x,\,y,\,z,\,p_x,\,p_y,\,p_z)\) で記述されます．しかし，加速器の横方向線形光学では，運動量の代わりにビーム軸からの軌道の傾き \((x^\prime=p_x/p_z,\,y^\prime=p_y/p_z)\) をよく使います．縦方向は，位相とエネルギー偏差 \((\phi,\,\Delta W)\) で表すことがよくあります．したがって実務では，\((x,\,x^\prime)\)，\((y,\,y^\prime)\)，\((\phi,\,\Delta W)\) という三組の二次元位相空間プロットが使われます．

ツイスパラメーター \((\alpha,\,\beta,\,\gamma)\) とエミッタンス \(\varepsilon\) を用いると，位相空間楕円は

\begin{align} \begin{bmatrix} x & x^\prime \end{bmatrix} \begin{bmatrix} \gamma & \alpha \\ \alpha & \beta \end{bmatrix} \begin{bmatrix} x \\ x^\prime \end{bmatrix} = \varepsilon \label{eq:html_twiss} \end{align}

と書けます．これが楕円になるためには，

\begin{align} 0 < \beta, \qquad 0 < \gamma, \qquad 0 < \beta\gamma-\alpha^2 \label{eq:html_ellipse_cond} \end{align}

が必要です．さらに加速器では，通常

\begin{align} \beta\gamma-\alpha^2 = 1 \label{eq:html_twiss_constraint} \end{align}

という規約を採用します．こうすると，式\eqref{eq:html_twiss} の右辺 \(\varepsilon\) は楕円の面積を \(\pi\) で割った量になり，エミッタンスの定義と一致します．粒子分布の RMS 楕円を扱う場合は，この \(\varepsilon\) が後で定義する \(\varepsilon_{rms}\) に対応します．

正規化座標

\begin{align} X=\frac{x}{\sqrt{\beta}}, \qquad P=\frac{\alpha x+\beta x^\prime}{\sqrt{\beta}} \label{eq:html_normalized_coord} \end{align}

を導入すると，式\eqref{eq:html_twiss} は \(X^2+P^2=\varepsilon\) になります．つまり，物理座標では傾いた楕円でも，正規化座標では円になります．このことから，\(\beta\) は位置方向の広がり，\(\gamma\) は傾き方向の広がり，\(\alpha\) は位置と傾きの相関を表す量と解釈できます．

位相空間分布の統計的なモデルとして，二次元正規分布は最も基本的な分布です．説明のため，ここでは変数を \((x,\,y)\) と書きます．加速器の位相空間へ戻すときは，\((x,\,x^\prime)\)，\((y,\,y^\prime)\)，\((\phi,\,\Delta W)\) などと読み替えます．

平均ベクトル \(\boldsymbol{\mu}\)，共分散行列 \(\boldsymbol{\Sigma}\) をもつ \(k\) 次元正規分布の確率密度関数は，

\begin{align} f_X(x_1,\cdots,x_k)= \frac{\exp\left[-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right]} {\sqrt{(2\pi)^k\lvert\boldsymbol{\Sigma}\rvert}} \label{eq:html_multi_normal} \end{align}

です．二次元の場合の共分散行列は，

\begin{align} \boldsymbol{\Sigma}= \begin{bmatrix} \sigma_x^2 & \rho\sigma_x\sigma_y \\ \rho\sigma_x\sigma_y & \sigma_y^2 \end{bmatrix} \label{eq:html_covariance} \end{align}

と書けます．ここで，\(\sigma_x\)，\(\sigma_y\) は標準偏差，\(\rho\) は相関係数です．通常の二次元正規分布として密度関数を書くには，共分散行列が正則でなければならないため，\(-1<\rho<1\) とします．\(\rho=\pm1\) では分布が直線上につぶれ，共分散行列の逆行列を使えません．

平均を \(\boldsymbol{\mu}=0\) とすると，確率密度関数は

\begin{align} f_X(x,y) &=\frac{1}{2\pi\sigma_x\sigma_y\sqrt{1-\rho^2}} \exp\left[ -\frac{1}{2(1-\rho^2)} \left( \frac{x^2}{\sigma_x^2} -\frac{2\rho xy}{\sigma_x\sigma_y} +\frac{y^2}{\sigma_y^2} \right) \right] \label{eq:html_normal2d} \end{align}

となります．したがって，等密度曲線は

\begin{align} \begin{bmatrix} x & y \end{bmatrix} \boldsymbol{\Sigma}^{-1} \begin{bmatrix} x \\ y \end{bmatrix} = \lambda \label{eq:html_density_ellipse} \end{align}

という楕円族になります．共分散行列の行列式を使うと，この楕円の面積は

\begin{align} S(\lambda)=\pi\lambda\sqrt{\det(\boldsymbol{\Sigma})} =\pi\lambda\sigma_x\sigma_y\sqrt{1-\rho^2} \label{eq:html_area} \end{align}

で与えられます．加速器の横方向位相空間で \(y\) を \(x^\prime\) に読み替えると，\(\sqrt{\det(\boldsymbol{\Sigma})}\) が後で定義する RMS エミッタンスになります．

粒子分布の統計量

粒子分布の統計量は，共分散行列で表すのが基本です．連続分布 \(f(x,\,x^\prime)\) が規格化されているとき，任意の関数 \(g(x,\,x^\prime)\) の期待値は

\begin{align} \langle g \rangle = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} g(x,\,x^\prime)f(x,\,x^\prime)\,\mathrm{d}x\,\mathrm{d}x^\prime \label{eq:html_expectation} \end{align}

で定義します．平均を引いた中心化座標を使えば，共分散行列は

\begin{align} \boldsymbol{\Sigma} = \begin{bmatrix} \langle x^2\rangle & \langle xx^\prime\rangle \\ \langle xx^\prime\rangle & \langle x^{\prime 2}\rangle \end{bmatrix} \label{eq:html_covariance_cont} \end{align}

となります．ここでは，\(x\) と \(x^\prime\) は平均を引いた後の量として扱います．

シミュレーションでは，通常 \(N\) 個のマクロ粒子 \((x_i,\,x_i^\prime)\) が与えられます．重みなしの場合の平均は

\begin{align} \bar{x}=\frac{1}{N}\sum_{i=1}^N x_i, \qquad \bar{x}^\prime=\frac{1}{N}\sum_{i=1}^N x_i^\prime \label{eq:html_mean_discrete} \end{align}

です．\(\tilde{x}_i=x_i-\bar{x}\)，\(\tilde{x}_i^\prime=x_i^\prime-\bar{x}^\prime\) を使うと，共分散行列は

\begin{align} \boldsymbol{\Sigma} &= \frac{1}{N} \sum_{i=1}^N \begin{bmatrix} \tilde{x}_i \\ \tilde{x}_i^\prime \end{bmatrix} \begin{bmatrix} \tilde{x}_i & \tilde{x}_i^\prime \end{bmatrix} \nonumber\\ &= \begin{bmatrix} \langle x^2\rangle & \langle xx^\prime\rangle \\ \langle xx^\prime\rangle & \langle x^{\prime 2}\rangle \end{bmatrix} \label{eq:html_covariance_discrete} \end{align}

で与えられます．ここでの共分散は，粒子分布そのものの二次モーメントとして扱います．したがって，統計的な不偏推定量で使う \(1/(N-1)\) ではなく，\(1/N\) を使います．

粒子に重み \(w_i\) がある場合は，総重み \(W=\sum_i w_i\) を用いて

\begin{align} \bar{x}_w=\frac{1}{W}\sum_i w_i x_i, \qquad \langle x^2\rangle_w=\frac{1}{W}\sum_i w_i(x_i-\bar{x}_w)^2 \label{eq:html_weighted_moment} \end{align}

のように，重み付き平均と重み付き二次モーメントを定義します．\(x^\prime\) と \(xx^\prime\) についても同様に計算します．

ツイスパラメーターと統計量の関係

まず，RMS エミッタンスを

\begin{align} \varepsilon_{rms} = \sqrt{\det(\boldsymbol{\Sigma})} = \sqrt{\langle x^2\rangle\langle x^{\prime 2}\rangle-\langle xx^\prime\rangle^2} \label{eq:html_rms_emittance} \end{align}

と定義します．この量は，共分散楕円の面積を \(\pi\) で割った値です．

ツイス楕円と共分散楕円を同じ楕円として比較すると，両者は同じ二次形式を別の記号で書いていることが分かります．具体的には，

\begin{align} \boldsymbol{\Sigma} = \varepsilon_{rms} \begin{bmatrix} \beta & -\alpha \\ -\alpha & \gamma \end{bmatrix} \label{eq:html_sigma_twiss} \end{align}

が成り立ちます．右辺の行列は，\(\beta\gamma-\alpha^2=1\) により行列式が 1 です．そのため，式\eqref{eq:html_sigma_twiss}の両辺の行列式を取ると，\(\det(\boldsymbol{\Sigma})=\varepsilon_{rms}^2\) となり，式\eqref{eq:html_rms_emittance}の定義と一致します．

この関係から，統計量を使ってツイスパラメーターを求めると，

\begin{align} \beta=\frac{\langle x^2\rangle}{\varepsilon_{rms}}, \qquad \alpha=-\frac{\langle xx^\prime\rangle}{\varepsilon_{rms}}, \qquad \gamma=\frac{\langle x^{\prime 2}\rangle}{\varepsilon_{rms}} \label{eq:html_twiss_from_statistics} \end{align}

で求められます．逆に，

\begin{align} \langle x^2\rangle=\beta\varepsilon_{rms}, \qquad \langle xx^\prime\rangle=-\alpha\varepsilon_{rms}, \qquad \langle x^{\prime 2}\rangle=\gamma\varepsilon_{rms} \label{eq:html_statistics_from_twiss} \end{align}

と書けます．この式から，\(\beta\) は位置方向の広がり，\(\gamma\) は傾き方向の広がり，\(\alpha\) は位置と傾きの相関を表す量と分かります．ただし，\(\beta\)，\(\alpha\)，\(\gamma\) は一般に寸法をもちます．そのため，「共分散行列を無次元化した量」ではなく，「共分散行列を RMS エミッタンスで規格化して得られる量」と考えるのが正確です．

統計学の記号 \((\sigma_x,\,\sigma_{x^\prime},\,\rho)\) を使うと，

\begin{align} \varepsilon_{rms}&=\sigma_x\sigma_{x^\prime}\sqrt{1-\rho^2}, \nonumber\\ \beta&=\frac{\sigma_x}{\sigma_{x^\prime}\sqrt{1-\rho^2}}, \nonumber\\ \alpha&=-\frac{\rho}{\sqrt{1-\rho^2}}, \nonumber\\ \gamma&=\frac{\sigma_{x^\prime}}{\sigma_x\sqrt{1-\rho^2}} \label{eq:html_eps_rho} \end{align}

となります．相関係数 \(\rho\) が正なら，この符号規約では \(\alpha\) は負になります．

実際の計算手順

マクロ粒子分布からツイスパラメーターを求めるときは，次の順序で計算すると間違いにくいです．

1. 粒子座標 \((x_i,\,x_i^\prime)\) の平均値 \((\bar{x},\,\bar{x}^\prime)\) を求めます．
2. 各粒子から平均値を引き，中心化座標 \((\tilde{x}_i,\,\tilde{x}_i^\prime)\) を作ります．
3. 中心化座標から \(\langle x^2\rangle\)，\(\langle xx^\prime\rangle\)，\(\langle x^{\prime 2}\rangle\) を計算します．
4. 式\eqref{eq:html_rms_emittance}で \(\varepsilon_{rms}\) を求めます．
5. 式\eqref{eq:html_twiss_from_statistics}で \(\alpha\)，\(\beta\)，\(\gamma\) を求めます．

この計算では，平均値を先に引くことが重要です．平均軌道がずれているまま二次モーメントを計算すると，ビームサイズではなく，軌道のオフセットまで含んだ値になります．

RMS エミッタンスと包含率

式\eqref{eq:html_rms_emittance}で定義した RMS エミッタンスは，共分散行列の行列式で決まる量です．そのため，共分散楕円

\begin{align} \begin{bmatrix} x & x^\prime \end{bmatrix} \boldsymbol{\Sigma}^{-1} \begin{bmatrix} x \\ x^\prime \end{bmatrix} =1 \label{eq:html_rms_ellipse} \end{align}

の面積は \(\pi\varepsilon_{rms}\) になります．ただし，RMS 楕円は「全粒子を含む境界」ではありません．分布の二次モーメントから決まる代表的な楕円です．

二次元正規分布に対し，楕円

\begin{align} \begin{bmatrix} x & x^\prime \end{bmatrix} \boldsymbol{\Sigma}^{-1} \begin{bmatrix} x \\ x^\prime \end{bmatrix} \leq \lambda \label{eq:html_probability_ellipse} \end{align}

の内部に含まれる粒子の割合は，

\begin{align} P(\lambda)=1-\exp\left(-\frac{\lambda}{2}\right) \label{eq:html_probability} \end{align}

で与えられます．RMS 楕円は \(\lambda=1\) に対応するので，その内部に入る粒子の割合は

\begin{align} P(1)=1-\exp\left(-\frac{1}{2}\right)\approx 0.393 \end{align}

です．つまり，二次元正規分布では，RMS 楕円の内側に入る粒子は約 39.3 % です．一次元の \(1\sigma\) の直感を，そのまま二次元へ持ち込まないことが大切です．

二次元正規分布での楕円係数 \(\lambda\) と包含率 \(P\) の関係です．RMS 楕円は \(\lambda=1\) で，包含率は約 39.3 % です．

任意の包含率 \(P\) に対応する楕円係数は，

\begin{align} \lambda=-2\log(1-P) \label{eq:html_lambda} \end{align}

で求められます．例えば 90 % の粒子を含む楕円では \(\lambda \approx 4.605\) です．この楕円の面積を \(\pi\) で割った量は，RMS エミッタンスの約 4.605 倍になります．

また，平均値を引いた粒子座標 \((x_i,\,x_i^\prime)\) が楕円の内側か外側かは，

\begin{align} \begin{bmatrix} x_i & x_i^\prime \end{bmatrix} \begin{bmatrix} \gamma & \alpha \\ \alpha & \beta \end{bmatrix} \begin{bmatrix} x_i \\ x_i^\prime \end{bmatrix} \left\{ \begin{array}{ll} < \lambda\varepsilon_{rms} & \textrm{楕円内部} \\ = \lambda\varepsilon_{rms} & \textrm{楕円上} \\ > \lambda\varepsilon_{rms} & \textrm{楕円外部} \end{array} \right. \label{eq:html_inside_outside} \end{align}

で判定できます．この判定式を使うと，粒子分布から指定した包含率の楕円に入る粒子を数えることもできます．

なお，式\eqref{eq:html_probability}と式\eqref{eq:html_lambda}は，分布が二次元正規分布の場合の結果です．実際のビーム分布が非正規分布の場合，RMS エミッタンスは計算できますが，包含率は粒子を直接数えて確認する必要があります．

このページの要点は次の通りです．

• 粒子分布の広がりと相関は，共分散行列で表せます．
• RMS エミッタンスは，共分散行列の行列式の平方根 \(\sqrt{\det(\boldsymbol{\Sigma})}\) です．
• ツイスパラメーターは，共分散行列を RMS エミッタンスで規格化して求めます．
• RMS 楕円の面積は \(\pi\varepsilon_{rms}\) ですが，二次元正規分布でその内側に入る粒子は約 39.3 % です．
• 包含率と楕円係数の関係は，正規分布を仮定したときだけ簡単な式で書けます．非正規分布では粒子を直接数える必要があります．

参考資料

更新履歴