ツイスパラメーター分布生成

ツイスパラメーターから粒子分布を作る作業は，ビームダイナミクス計算の入口です．入力として欲しいのは，一つの横方向平面のツイスパラメーター \((\alpha,\,\beta)\) と RMS エミッタンス \(\varepsilon_{rms}\) です．出力として欲しいのは，その統計量をもつ多数のマクロ粒子 \((x_i,\,x_i^\prime)\) です．

楕円の意味，二次元正規分布の基本式，統計量とツイスパラメーターの関係は，楕円，正規分布，統計 のページで説明しています．このページではそれらを必要な範囲だけ使い，分布をどのように作るかに話を絞ります．

結論は単純です．まず平均 0，共分散行列が単位行列の標準二次元正規分布を作ります．次に，目標の共分散行列をコレスキー分解して得た下三角行列を掛けます．これにより，入力したツイスパラメーターと RMS エミッタンスをもつ分布が得られます．平均がゼロでないビームを作りたい場合は，最後に平均値を加えます．

標準正規乱数から，コレスキー変換を使って目標のツイスパラメーターをもつ分布を作る流れです．

基本

考え方の中心は，線形変換と共分散行列の関係です．平均 0，共分散行列が単位行列の二次元確率変数を

\begin{align} \boldsymbol{z}_i= \begin{bmatrix} u_i \\ v_i \end{bmatrix}, \qquad \mathrm{Cov}(\boldsymbol{z}_i)= \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} =\boldsymbol{I} \label{eq:standard_normal_vector} \end{align}

とします．目標の共分散行列を \(\boldsymbol{\Sigma}_{target}\) とし，ある行列 \(\boldsymbol{L}\) が

\begin{align} \boldsymbol{\Sigma}_{target}=\boldsymbol{L}\boldsymbol{L}^T \label{eq:cholesky_basic} \end{align}

を満たしているとします．このとき，

\begin{align} \boldsymbol{w}_i=\boldsymbol{L}\boldsymbol{z}_i \label{eq:linear_transform} \end{align}

で作った新しい確率変数の共分散行列は

\begin{align} \mathrm{Cov}(\boldsymbol{w}_i) &= \boldsymbol{L}\,\mathrm{Cov}(\boldsymbol{z}_i)\,\boldsymbol{L}^T \nonumber\\ &= \boldsymbol{L}\boldsymbol{I}\boldsymbol{L}^T = \boldsymbol{L}\boldsymbol{L}^T = \boldsymbol{\Sigma}_{target} \label{eq:cov_linear_transform} \end{align}

になります．つまり，標準正規分布の円形の雲を，行列 \(\boldsymbol{L}\) で伸ばして傾ければ，欲しい共分散行列をもつ分布が得られます．この \(\boldsymbol{L}\) を正定値対称行列に対して系統的に求める方法が，コレスキー分解です．

コレスキー行列 \(\boldsymbol{L}_{tw}\) は，標準正規分布を目標の共分散行列をもつ分布へ変換します．

ツイスパラメーターの共分散行列

ツイスパラメーターと RMS エミッタンスが表す目標の共分散行列は，

\begin{align} \boldsymbol{\Sigma}_{tw} = \varepsilon_{rms} \begin{bmatrix} \beta & -\alpha \\ -\alpha & \gamma \end{bmatrix}, \qquad \beta\gamma-\alpha^2=1 \label{eq:twiss_covariance} \end{align}

です．通常は \(\beta>0\)，\(\varepsilon_{rms}>0\) です．また，\(\beta\gamma-\alpha^2=1\) なので，\(\gamma=(1+\alpha^2)/\beta\) です．この条件の下で \(\boldsymbol{\Sigma}_{tw}\) は正定値対称行列になり，コレスキー分解できます．

下三角行列を \(\boldsymbol{L}_{tw}\) とすると，一般には

\begin{align} \boldsymbol{L}_{tw} = \sqrt{\varepsilon_{rms}} \begin{bmatrix} \sqrt{\beta} & 0 \\ -\cfrac{\alpha}{\sqrt{\beta}} & \sqrt{\gamma-\cfrac{\alpha^2}{\beta}} \end{bmatrix} \label{eq:twiss_cholesky_general} \end{align}

です．ツイスパラメーターの拘束条件 \(\beta\gamma-\alpha^2=1\) を使うと，

\begin{align} \boldsymbol{L}_{tw} = \sqrt{\varepsilon_{rms}} \begin{bmatrix} \sqrt{\beta} & 0 \\ -\cfrac{\alpha}{\sqrt{\beta}} & \cfrac{1}{\sqrt{\beta}} \end{bmatrix} \label{eq:twiss_cholesky_simple} \end{align}

と簡単になります．実装ではこの形が便利です．\(\gamma\) を独立に入力するより，\(\gamma=(1+\alpha^2)/\beta\) として計算した方が，パラメーターの不整合を避けやすくなります．

ここでは，一つの横方向平面 \((x,\,x^\prime)\) の分布を作ります．\(y\) 平面や \((\phi,\,\Delta E)\) に対しても，同じ考え方を使えます．ただし，平面間の相関を入れたい場合は，二次元ではなく多次元の共分散行列を扱う必要があります．

手続き 1: 標準二次元正規分布

最初に，互いに独立な標準正規乱数 \(u_i\) と \(v_i\) を \(N\) 個ずつ生成します．これを

\begin{align} \boldsymbol{z}_i= \begin{bmatrix} u_i \\ v_i \end{bmatrix}, \qquad u_i \sim \mathcal{N}(0,\,1), \qquad v_i \sim \mathcal{N}(0,\,1) \label{eq:standard_pair} \end{align}

と書きます．この分布は原点中心の円形のガウス分布であり，共分散行列は単位行列です．言い換えると，\(u\) 方向と \(v\) 方向に同じ広がりをもち，相関がない状態から出発します．

ここで重要なのは，作りたいものが楕円の境界線ではなく，楕円状に広がるガウス分布だという点です．そのため，半径と角度から楕円周上の点を作るのではなく，まず正規乱数を使って点群を作ります．

手続き 2: コレスキー分解

次に，入力された \(\alpha\)，\(\beta\)，\(\varepsilon_{rms}\) から式\eqref{eq:twiss_cholesky_simple}の下三角行列 \(\boldsymbol{L}_{tw}\) を作ります．この行列が，標準正規分布の円形の雲を，目的の広がりと傾きをもつ雲へ変形します．

行列の形をもう一度書くと，

\begin{align} \boldsymbol{L}_{tw} = \sqrt{\varepsilon_{rms}} \begin{bmatrix} \sqrt{\beta} & 0 \\ -\cfrac{\alpha}{\sqrt{\beta}} & \cfrac{1}{\sqrt{\beta}} \end{bmatrix} \label{eq:twiss_cholesky_for_use} \end{align}

です．1 行目は位置 \(x\) の広がりを決めます．2 行目は \(x\) と \(x^\prime\) の相関，および傾き方向の広がりを決めます．この符号規約では，\(\alpha=0\) なら相関はなく，\(\alpha<0\) なら \(x\) と \(x^\prime\) は正に相関し，\(\alpha>0\) なら負に相関します．

手続き 3: 座標変換

最後に，各粒子に対して

\begin{align} \begin{bmatrix} x_i \\ x_i^\prime \end{bmatrix} = \boldsymbol{L}_{tw} \begin{bmatrix} u_i \\ v_i \end{bmatrix} \label{eq:final_transform} \end{align}

を計算します．成分で書けば，

\begin{align} x_i &= \sqrt{\beta\varepsilon_{rms}}\,u_i \label{eq:generated_x} \\ x_i^\prime &= \sqrt{\cfrac{\varepsilon_{rms}}{\beta}}\,(v_i-\alpha u_i) \label{eq:generated_xp} \end{align}

です．平均位置 \(x_0\) と平均傾き \(x_0^\prime\) をもたせたい場合は，最後に

\begin{align} x_i \leftarrow x_i+x_0, \qquad x_i^\prime \leftarrow x_i^\prime+x_0^\prime \label{eq:add_centroid} \end{align}

と平行移動します．この平行移動は分布の中心を変えるだけなので，中心化した二次モーメントや RMS エミッタンスは変わりません．

生成後の確認

生成された分布の二次モーメントは，式\eqref{eq:generated_x}と式\eqref{eq:generated_xp}から確認できます．\(u_i\) と \(v_i\) は独立で，平均 0，分散 1 なので，理想的には

\begin{align} \langle x^2 \rangle &= \beta\varepsilon_{rms} \\ \langle xx^\prime \rangle &= -\alpha\varepsilon_{rms} \\ \langle x^{\prime 2} \rangle &= \cfrac{1+\alpha^2}{\beta}\varepsilon_{rms} = \gamma\varepsilon_{rms} \label{eq:generated_moments} \end{align}

となります．したがって，共分散行列は確かに式\eqref{eq:twiss_covariance}になります．散布図を描いたときの RMS 楕円は，入力したツイスパラメーターで決まる楕円になります．

実装だけが必要なら，次のような Python プログラムで十分です．この例では numpy を使い，生成した粒子を particles.txt に保存します．

#!/usr/bin/env python3
import numpy as np

def sample_twiss(alpha, beta, eps_rms, n, x0=0.0, xp0=0.0, seed=1234):
    if beta <= 0.0:
        raise ValueError("beta must be positive.")
    if eps_rms <= 0.0:
        raise ValueError("eps_rms must be positive.")

    rng = np.random.default_rng(seed)
    u = rng.standard_normal(n)
    v = rng.standard_normal(n)

    x = np.sqrt(beta * eps_rms) * u + x0
    xp = np.sqrt(eps_rms / beta) * (v - alpha * u) + xp0
    return np.column_stack((x, xp))


def calc_twiss(particles):
    x = particles[:, 0] - np.mean(particles[:, 0])
    xp = particles[:, 1] - np.mean(particles[:, 1])

    xx = np.mean(x * x)
    xxp = np.mean(x * xp)
    xpxp = np.mean(xp * xp)
    eps_rms = np.sqrt(xx * xpxp - xxp ** 2)

    alpha = -xxp / eps_rms
    beta = xx / eps_rms
    gamma = xpxp / eps_rms
    return alpha, beta, gamma, eps_rms


def main():
    alpha = 1.65
    beta = 0.051
    eps_rms = 2.894352e-5
    n = 10000

    particles = sample_twiss(alpha, beta, eps_rms, n)
    np.savetxt("particles.txt", particles, fmt="%.6e", header="x xp")

    alpha2, beta2, gamma2, eps2 = calc_twiss(particles)
    print(f"alpha   = {alpha2:.6f}")
    print(f"beta    = {beta2:.6f}")
    print(f"gamma   = {gamma2:.6f}")
    print(f"eps_rms = {eps2:.6e}")
    print("particles.txt written")


if __name__ == "__main__":
    main()

この例では，particles.txt の 1 列目が \(x\)，2 列目が \(x^\prime\) です．平均位置 \(x_0\) や平均傾き \(x_0^\prime\) をもたせたい場合は，sample_twiss() の引数に値を与えます．生成後に calc_twiss() で統計量を計算し直すと，入力値に近い \(\alpha\)，\(\beta\)，\(\gamma\)，\(\varepsilon_{rms}\) が得られることを確認できます．

実装上の注意

• 生成されるのは楕円の境界線ではありません．得られるのは楕円状のガウス分布です．RMS 楕円は，その統計量を表す補助線です．
• 粒子数が有限だと統計量は揺らぎます．\(N\) が小さいと，生成後に計算し直した \(\alpha\)，\(\beta\)，\(\varepsilon_{rms}\) は入力値と完全には一致しません．粒子数を増やすとずれは小さくなります．
• 入力の整合性を確認します．少なくとも \(\beta>0\) と \(\varepsilon_{rms}>0\) は必須です．\(\alpha\)，\(\beta\)，\(\gamma\) を別々に与える場合は，\(\beta\gamma-\alpha^2=1\) も確認します．
• 独立な 3 平面なら同じ手続きを 3 回行います．\(x\)，\(y\)，\((\phi,\,\Delta E)\) を独立とみなす場合は，各平面ごとに同じ計算を行います．平面間の相関まで入れたい場合は，6 次元の共分散行列を直接扱います．
• 正規分布以外を作りたい場合は別の乱数を使います．このページの方法は，標準正規乱数を線形変換するので，生成される分布も正規分布です．一様分布や切断分布を作る場合は，出発点となる乱数の作り方を変える必要があります．

ツイスパラメーターから二次元正規分布を作る手順は，次の 3 段階です．

1. 独立な標準正規乱数から，原点中心で相関のない二次元分布を作ります．
2. ツイスパラメーターが表す共分散行列をコレスキー分解し，下三角行列 \(\boldsymbol{L}_{tw}\) を求めます．
3. \(\boldsymbol{w}_i=\boldsymbol{L}_{tw}\boldsymbol{z}_i\) により座標変換し，必要なら最後に平均値を加えます．

この方法の利点は，理屈が明快で，実装が短く，生成後の統計量も確認しやすいことです．より詳しい導出や，実際のプログラムと入出力ファイルを含めた説明は，同じディレクトリに置いた PDF を参照してください．

参考資料

更新履歴