2017-09-26 KEKセミナー

Bulk-To-Point 法を用いた極低エミッタンスビーム評価用トラッキングシミュレーションコード

(公財)高輝度光科学研究センター

水野 明彦

昨今のERL, XFEL等では極低エミッタンス電子ビーム源が要である↓
その設計には,低エネルギー電子用のビームトラッキングシミュレーションコードが必須である↓（軽い気持ちで） 自分でコードを作ってみた↓（意外に） 他のコードと一致しない ⇐ 非相対論領域の空間電荷効果↓
真面目にコード間で比較をしてみる ⇒ 全然一致しない。これは大変だ！

Contents
0. Point-To-Point コードと、その精度
Point-To-PointコードからBulk-To-Pointコードへ

1. Multiple Beam Envelope Equations
エンベロープ方程式について
multiple envelope equations

2. 2D Bulk-To-Point Tracking Codeエンベロープ方程式からシミュレーションコードへ

3. 3D Bulk-To-Point Tracking Code
3次元への拡張

4. DC continuous beamのシミュレーション
シミュレーションモデルのビームダイナミクス
精度検証

5. GUI for Bulk-To-Point Tracking Code

0. Point-To-Point コードと、その精度

これまでは ３次元Point-To-Point Code を開発してきた
Point-To-Point == 空間電荷相互作用として 点電荷である各マクロ粒子間の相互作用を全て計算
遅延ポテンシャルを用いて実験室系で計算粒子の相互作用vB,
B ee-r

EA 1

 er

3224 0 2 vB r 2r 
2c



1BA  2 vB  EAc

ただし、マクロ粒子Bは等速直線運動であることを仮定

AvAFA  eE A  v A  B A 

運動方程式の扱い

エミッタンスの計算d v dP e v  B  E  m0dtdtdve  v  E  v  B  E  2 vdt
m0 c𝜀𝑥 = 𝛾𝛽𝑥 2 𝑥′2 − 𝑥 ∙ 𝑥′ 2４次精度のRunge-Kutta method
A. Mizuno et al.,Nucl. Instrum. Methods Phys. Res., Sect. A528, 387 (2004)

Point-To-Point Codeによるフリースペースでのトラッキング(100mm)••••

初期バンチ半径 : φ1mm uniform
初期バンチ長: 3mm uniform
電荷量: 1nC/Bunch
ビームエネルギー： 1MeV & 6MeV

計算時間：
18時間37分@20万個
1 MeV

6 MeV

Intel Xeon W5590 3.33GHz
8 coreGPTで採用

•••

ベンチマークテストの状況 （フリースペース：300mm）

1 MeV

0.5 MeV

6 MeV

初期バンチ半径 : φ1mm uniform
初期バンチ長: 3mm uniform
電荷量: 1nC/Bunch

水野コード
Paｒmela(菅)山本尚人コード
GPT(山本昌志)
KUBLAI(増田): 3次元 Point-To-Point : n=2×105: 2次元 : n=5001（参照粒子を含む）: 3次元 Point-To-Point : n=5000: 3次元 PIC : n=1×105: 2次元 PIC : n=8704

どれが正しいの ？？

Point-To-Point法は、任意のビーム形状を計算でき有用であるが、点電荷からの電磁場を計算する

点電荷同士が接近した場合、非常に大きな相互作用が発生し、計算精度が悪化する

粒子数を増大させれば精度が良くなるものの、計算時間が増大して現実的では無い

点電荷ではなく、バルク電荷からの電磁場を計算した方が良さそう1. Multiple Beam Envelope Equations

Multiple Beam Envelope Equations （2次元コード）
Longitudinal Space Charge Fields• z軸上のfieldのみ解析的に計算• スライスは歪まない• 各スライス内の電子密度分布は均一• Longitudinal: m, Transverse: n に分割, 各交点に電子を配置• 各スライスからの Space Charge Fields を計算• 各電子に対応するエンベロープ方程式を連立して解く

ビームエンベロープ方程式 （連続ビーム）
電子の運動方程式

ペンシルビーム

空間電荷効果外部電磁場
連続ビーム、またはビームエネルギー無限大の場合
アンペールの法則

Lawsonのエンベロープ方程式

ビームエンベロープ方程式 （単バンチビーム transverse方向）

パルスのペンシルビーム

円柱ビームが誘起する電磁場

5変数

3変数

3次元データファイルマップを用意しておいてfieldを計算

ビームエンベロープ方程式 （単バンチビーム longitudinal方向ビーム軸上）

解析的に計算可能

エンベロープ方程式

空間電荷効果外部電磁場

Multiple Beam Envelope Equations

For transverse direction
For transverse direction
For longitudinal directionFor energy

各微分方程式を連立して Octave で解く

A. Mizuno et al., Phys. Rev. ST Accel. Beams 15, 064201(2012)

エミッタンスの計算方法

各期待値、𝑟 −エミッタンス : 𝜀𝑟 = 𝛾𝛽𝑥 −エミッタンス : 𝜀𝑥 = 𝜀𝑟 /2

𝑟 2 𝑟′2 − 𝑥 ∙ 𝑥′ 2

Simulation cases #1: BNL RFgun Cavity and Free space
Multiple Beam Envelope Equations (0.1 nC, 𝑚 = 𝑛 = 10)

Initial Parameters:
レーザースポットサイズ
レーザー幅
電荷量
エネルギー@cavity exit

：f1.3 mm：10 ps：0.1 nC：5.75 MeV

Simulation cases #1: BNL RFgun Cavity and Free space
Multiple Beam Envelope Equations (1.0 nC, 𝑚 = 𝑛 = 10)

Initial Parameters:
レーザースポットサイズ
レーザー幅
電荷量
エネルギー@cavity exit•••

：f1.3 mm：10 ps：1.0 nC：5.75 MeV

z軸上のfieldのみ解析的に計算
スライスは歪まない
各スライス内の電子密度分布は均一

Simulation cases #2: Spring-8 RFgun Cavity+ Solenoid + Acc, overfocus
Multiple Beam Envelope Equations (50.0 pC, 𝑚 = 𝑛 = 10)••

スライスは歪まない
各スライス内の電子密度分布は均一

Initial Parameters:
レーザースポットサイズ
レーザー幅
電荷量
エネルギー@cavity exit
エネルギー@Acc. Exit

：f1.2 mm
：20 ps：50.0 pC
：3.7 MeV：30.0 MeV

2. 2D Bulk-To-Point Tracking Code

2. 2D Bulk–To–Point Tracking Code （2次元コード）
新たにOff-axisのLongitudinal Fieldsも計算• 複数のドーナツで構成されるスライスは歪んでも良い• 各ドーナツ内の電子密度分布が均一• 複数のドーナツにより、径方向電子密度分布を表現可能
径方向Off-axis電場の計算• Longitudinal: m, Transverse: n に分割, 各交点に電子を配置• 各ドーナツからの Space charge field を計算• 各電子の運動方程式を数値的（ルンゲクッタ法）に解く

A. Mizuno, Phys. Rev. Accel. Beams 19, 024201 (2016)

Simulation cases #1: BNL Cavity and Free space (1.0 nC)
2D Bulk-To-Point Tracking Code (1.0 nC, 𝑚 = 𝑛 = 20)

Multiple Beam Envelope Equations

Initial Parameters:
レーザースポットサイズ
レーザー幅電荷量エネルギー@cavity exit

2D Bulk-To-Point Tracking Code：f1.3 mm：10 ps：1.0 nC：5.75 MeVMultiple Beam Envelope Equations

2D Bulk-To-Point Tracking Code

Simulation cases #2: Spring-8 RFgun Cavity+ Solenoid + Acc, overfocus

2D Bulk-To-Point Tracking Code (50.0 pC, 𝑚 = 𝑛 = 20)

Initial Parameters:
レーザースポットサイズ
レーザー幅
電荷量
エネルギー@cavity exit
エネルギー@Acc. Exit

Multiple Beam Envelope Equations

2D Bulk-To-Point Tracking Code：f1.2 mm：20 ps：50.0 pC
：3.7 MeV：30.0 MeVMultiple Beam Envelope Equations

2D Bulk-To-Point Tracking Code

Simulation cases #3: Spring-8 RFgun Cavity+ Solenoid + Acc, 計算精度と計算時間
2D Bulk-To-Point Tracking Code (50 pC)
Initial Parameters:
レーザースポットサイズ
レーザー幅
電荷量
エネルギー@cavity exit
エネルギー@Acc. Exit𝑛 = 𝑚 = 100𝑛 = 𝑚 = 20

：f1.4 mm
：9 ps
：50.0 pC
：3.7 MeV：30.0 MeV𝑛 = 2 × 105

𝑛 = 5 × 103

計算時間 by 8-core Xeon W5590 3.33 GHz• Point-To-Point Code (𝑛 = 2 × 105 ): 3 weeks• Bulk-To-Point Tracking Code ( 𝑚 = 𝑛 = 100) : 24 hours• Bulk-To-Point tracking Code ( 𝑚 = 𝑛 = 20) : 3 min3. 3D Bulk-To-Point Tracking Code3. 3D-Bulk–To–Point Tracking Code （準3次元コード）• 電子を 𝑥軸上, および 𝑦軸上に配置• 各楕円ドーナツからの Space charge field を計算• ２Dコードと同様に、各電子の運動方程式を数値的に解く

アスペクト比 𝛼 =𝑥を追加して4次元マッピングデータが必要？
𝑦𝑥𝑦のアスペクト比: a4次元フィールドマップ ⇒ 3次元フィールドマップ

A,B,C,Dでフィールドをフィッティング

A. Mizuno and H. Hanaki, Proceedings of Linac2016, p120-122, East Lansing, USA (2016)

Parameters:
Simulation cases #4: Spring-8 RFgun Cavity Initial
レーザースポットサイズ+ Q-Doublet + Acc レーザー幅

3D Bulk-To-Point Tracking Code (0.2 nC)

電荷量エネルギー@cavity exit
エネルギー@Acc. Exit

：f1.2 mm：20 ps：0.2 nC：3.7 MeV：30.0 MeV

Point-To-Point Code :𝑛 = 6 × 104
3D Bulk-To-Point Code :𝑚 = 𝑛 = 20

4. DC continuous beamのシミュレーション
空間電荷効果によるエミッタンスリダクション

スライスエミッタンスリダクションの事例

径方向： Φ2mm 一様分布
長手方向
： 40A/cm2 一様な連続ビーム
初期熱エミッタンス ： 0 mm mrad
理想的な平面カソード
カソード鏡像効果を考慮
一様電界5MV/mで引き出して加速し続ける

simulated by KUAD2
京大エネルギー理工学研究所
増田准教授

ある位置zを通過した粒子の積算によりzでのエミッタンスを算出

rms エミッタンス ≠ Liouville’s volume

プロジェクションエミッタンスではなく、スライスエミッタンスに近いエミッタンスが低減しても理論的矛盾はないrmsエミッタンスの保存 (𝑟 − 𝑟 ′ 空間)
四角形Oacbの面積が保存される条件は、四角形の面積面積の増分
Linear forceが働く場合はO面積(エミッタンス)が保存

Linear force以外の作用がある場合は、エミッタンスは保存せず、したがって、減少する場合があるかもしれないエミッタンスの増減は、W>0としての正負で判断する

常に 𝑊 = 𝑎 × 𝑏 = 𝑟𝑎 𝑟𝑏′ − 𝑟𝑎′ 𝑟𝑏 > 0となるようにベクトル𝑎, 𝑏を配置 ⇒ 常に𝜃𝑎 < 𝜃𝑏

1. 𝑟𝑎 > 𝑟𝑏 の場合（上に凸）電子が受ける力を2. 𝑟𝑎 < 𝑟𝑏 の場合（下に凸）と表現したとき、それぞれで表される作用を受けて生成された位相空間分布

1. 𝑟𝑎 > 𝑟𝑏 の場合（上に凸）で表される作用を受けて生成された位相空間分布
ビーム中心の電荷密度が大きく、外側が小さい一般的な場合に相当

粒子a、bに対し、それぞれA. Kb > Kaの場合

エミッタンス増加すなわち、上に凸の形で表される力を受け続ける場合に相当の力が作用が作用する場合を考える

B. Kb < Kaの場合

エミッタンス減少

すなわち、下に凸の形で表される力を新たに受ける場合に相当

中心の密度が大きいビームが、外側の密度が大きな状態に移行することを意味する。
空間電荷効果によっては、通常は起こり得ないと考えられる。

2. 𝑟𝑎 < 𝑟𝑏 の場合（下に凸）で表される作用を受けて生成された位相空間分布
ビーム中心の電荷密度が小さく、外側が大きいような場合。
又は、電荷密度が一定であっても、カソード近傍で実現。

ra > rb の場合と同じく、であるので

A. Kb > Kaの場合

エミッタンス増加

すなわち、下に凸の形で表される力を受け続ける場合に相当

B. Kb < Kaの場合

エミッタンス減少

すなわち、上に凸の形で表される力を新たに受ける場合に相当

この状態は、カソード近傍で起こりえる。

エミッタンスが減少する条件
位相空間作用する径方向の力• rmsエミッタンスは変化• Liouvilleの面積 = 0, Consatnt
ororReduction of transverse emittance in electron injectors
caused by space charge effects
Mizuno, A., Masuda, K., Yamamoto, M.
Nucl. Instrum. Methods Phys. Res.,Sect. A774,51(2015)

𝑟 −エミッタンス

𝑥 −エミッタンス

空間電荷効果が作用する系では、𝑥 − 𝑥 ′ 空間ではConstantなハミルトニアンが定義できない⇒ 𝑥 − 𝑥 ′ 空間の面積も保存しない

Ｂ、Ｃ、Ｄ位置での𝑟 − 𝑟 ′ 空間の曲がり具合は僅か

非常に高精度な計算が必要

Simulation cases #5:
Accelerated DC Continuous Beam

In GPT

Using Spacecharge2Dcircle element
No Image charge on the cathode

径方向： Φ2mm 一様分布
長手方向： 40A/cm2 一様な連続ビーム
初期熱エミッタンス ： 0 mm mrad
理想的な平面カソード
カソード鏡像効果を考慮
一様電界5MV/mで引き出して加速し続ける

M100N10 : 4分
M200N20 : 46分
M400N40 : 12時間

Bulk-To-Point Codeメリット• パルスビームと連続ビームの両方が計算可能• 計算速度が速い• 広い領域でも計算可能
デメリット• Wake Field 等、解析的に計算できない電磁場を扱えない

PIC Code (KUAD2)
メリット• 有限要素法を用いた電磁場計算の正確性• Wake Field 等も扱うことが可能
デメリット• 連続ビームのみ計算可能。パルスビームは不可能• 広い計算領域では計算時間が増大し現実的では無い

5. GUI for Bulk-To-Point Tracking Code

Linuxで稼働
GUIはTcl/Tkで作成

まとめ• Bulk-To-Point法を用いて空間電荷効果を計算するビームトラッキングコードを開発した• エンベロープ方程式を連立して解くコード• スライスが歪んではならないために• 空間電荷効果が小さな場合は高精度• 空間電荷効果が大きい場合、ビームをオーバーフォーカスした場合等では、精度が悪い• バンチをドーナツに分けて各電子を数値的にトラッキングする2次元Bulk-To-Pointコード• 歪んだスライスもドーナツで表現可能• 空間電荷効果が大きい場合や、ビームを絞った場合でも高精度、短時間計算を実現した• 2次元Bulk-To-Pointコードを3次元に拡張した• DC Continuous Beamのエミッタンスリダクションを例に、Bulk-To-Pointコードの高精度を検証した• Bulk-To-Point コードについて、GUIを作成した