mini-symposium
hctak について
青木 保夫
重イオン用 CDCC を中心としたプログラム群

Fortran 90 を使用
zero-base で書いた
http://www.tac.tsukuba.ac.jp/∼yaoki
にて公開中

1

代表的なプログラムとその機能
opm: 中心力光学ポテンシャル
微分断面積/ポテンシャルパラメータ
wsf: WS ポテンシャル下の波動関数
束縛及び散乱状態
共鳴状態探索束縛状態波動関数の Fourier 展開
hctak: CDCC 方程式を解きS 行列を計算する
J π 毎の並列計算が可能 (MPI)
elx: 上記 S 行列 → 弾性散乱微分断面積
trix: 上記 S 行列 → 三重微分断面積
三重微分断面積を多重立体角積分
umtab: 核種 (基底状態) のデータベース
nucspe: 反応記号の解読(p+22Mg)+12C, E(Al)=1702 MeV
hctak.pdf: 利用マニュアル2

projectile

1r

R

2

3 targetint.

pot. typeV 12

real pot.

bound/ scatt.

V 13

OP pot.

scatt.

V 23

OP pot.

scatt.

3statesHamiltonian H:
H = Tr + V12(r)+ TR + V13(r13) + V23(r23)wave function φ(r) be defined by{Tr + V12(r) − Ec} φc(r) = 0
scatt. state wf. at large rr2 sin(k r − lcπ/2 + δc k )φc(k, r) →πrklc
δc kwave numberorb. ang. mom. for ch. cnucl. phase shift

if charge and spin are neglected4

truncation
discretization

k and spin space
k spacehenceContinuum Discretized
Coupled Channels (CDCC)
Z k1jφ̂c j = pφc(k, r) dkkj − kj−1 kj−1
and is orthonormalized ashφ̂c j | φ̂c0 j 0 ir = δc c0 δj j 0suffix
c, c’j, j’meaningspins
wave numbers

5eigen state of H withangular momentum J , M
ΨJ M =1X
RχJL (R)cjc, j× [φ̂c j (r) ilc Ylc (r̂) iLc j YLc j (R̂)]J M[ . . . ] for ang. mom. couplingχJL (R) : motion of (1-2) and 3,cj

satisfies the CDCC eq.
(TR − Ec j )χJL (R)cj=−X

h[c j]J |(V13 + V23)|[c0 j 0]J i

c0 j 0× χJL0c0 j 06(R)CDCC eq. is solved numericallywith usual boundary cond.
χJLc j (R) →ILc0 δc,c0 δLc j ,Lc0s
Kc 0 J−
Sc j, c0 Lc Oc j0Kc j
where c0 stands for inc. ch.
IL (OL) are usual Coul. wf.

7

Elastic cross sec
X

σel ∝ |fc+(geom. factor) ei(σL0 +σL)³´J2× SL,L−δY(R̂)|L L0LM0

triple diff. cross sec.
d3σdΩ1 dΩ2 dE1T ∝∝ (final state dens.)×|T |2ÃX( geom. factor)ei(δc k +σc k )JMc³´J× ei(σL0 +σLc ) SL(k)c L0× [Ylc (r̂) ⊗ YLc (R̂)]J M

8k!

歴史的な流れ八尋・井芹・桜木さん等物理学会での講演
RCNP 偏極国際会議での上村氏の講演
上村氏から 青木へソース提供
使用開始 < −− 井芹さんの協力

CDCC(d,p) ?
(−)

(+)

T (d, p) = hχp ΦB |Vnp|ΨCDCC ΦAi
ΨCDCC の基底状態重心運動部分をTWOFNR へ読み込ませる(テンソル力を受け付ける)

9

Reid soft core pot. での散乱状態
井芹 3S1-3D1青木 3S1-3D1, 3D2

3P , 3P , 3P -3F
0122‘-’ で繋いだ状態は、結合状態
3D 状態は Gauss int.
3重陽子専用 tak.f 開発
上村版より高速、高精度

f (CE; c, c0; λ) ∝ Rλ,相反性 Sa b = Sb a

R→0

３次スプラインによる行列要素の内挿
Stömer の８点法
初期条件設定による uJc の発散回避
S 行列計算の次元低下と条件数評価

10

Reid soft core pot.
原点付近での振舞
湯川理論だと 1/r 3(1)h2 (ix) =

µ−xex

1+

33

¶+ 2x x

ここで、x = µ r ∼ 0.7 r補正項を加えて、1/r としている。

Schrödinger eq. は確定特異点を持ち、原点付近では非常に硬い (stiff )級数展開項と 対数項とで二つの独立解を書き下せる。

Lattice QCD の核力
テンソル部分は原点で０湯川理論と矛盾している！！
singularity は簡単に消せる？

11

非常にテンソル力が強いため数値計算で α 状態と β(γ) 状態の弁別するのが困難

d2dx2

Ã
uw!Ã

Ã!

S TT D=!
uw

原点付近で初期値を(u = ², w = 0)
(u = 0, w = ²)と与えて数値的に解いても二つの解の直交性が極端に悪い。

12

V =VC + VT S1 2 + VLS L · S
VC = − 10.463 Y (x) + 210.936Y (2 x)− 12715.2 Y (4 x) + 59545.8 Y (6 x)
VT = − 10.463{(1 + 3/x + 3/x2)Y (x)− (48/x + 12/x2)Y (4x)}+ 1407.8Y (4x) − 10041Y (6x)VLS =2835.64Y (4x) − 16278.6Y (6x)
ここで、Y (x) ≡ e−x/x13

|pot| [MeV]

1.0E+05
1.0E+04

Reid soft core pot

1.0E+0333

S1- D1 potential

1.0E+02
1.0E+01
1.0E+00-VT(x)

1.0E-01-Vc(x)

1.0E-02
1.0E-03+VLS(x)
00.5

1

1.5

40.0

2x

2.5

3

3.5

43.5

43

S1-3D1 potential

VcentralReid soft core pot

20.0
[MeV]

Vspin-orb
0.0

E=-2.22 MeVVtensor-20.0

-40.0
00.5

11.52x142.53発散回避
大きな軌道角運動量の部分波 波動関数の振幅は、積分中に巨大振幅になる
積分の立ち上げ半径を x jL(x) の振幅が１ st ピークでの振幅の10−8 となる半径からとする

15

S 行列計算の次元低下と条件数評価
CDCC eq. の独立解u~ i, (i = 1, 2 · · · n) の重ね合わせU ≡ (~u1 u~2 · · · u~ n)
漸近領域での接続条件

Ã!Ã

U O
U 0 O0!

WS

Ã
=!
II0

W は重ね合わせ係数、S は S 行列
O (I) は外 (内) 向きクーロン波動関数を対角要素とする行列(O 0 U − O U 0) W = 2 i K
K は波数を対角要素とする行列次元数が半分になった行列 (O 0 U − O U 0) の条件数 (１ノルム) を評価する→ 計算精度の確認

16

ctak.f へ拡張

クーロン分解過程を考慮
真空偏極 (Uehling) ポテンシャル
Pong-Austern 流の反対称化
三重微分断面積 ← 井芹さんの協力

応用
重陽子：弾性散乱：０度での三重微分断面積
重陽子の polarizability

CDCC approach to
Oppenheimer andPhillips process17微分方程式を解きながら直交化

uJc (R) ∝(kR)J +1(2J + 1)!!

kc0 /kc = 1/2,µkc0

J = 100 の時、¶100kc= 10−30

１５桁の有効数字では取り扱えない！

CDCC 方程式は線形だから、解きながら直交化する
この手法は原理的には有効だが、効率的なアイデアは知らない

18真空偏極: Uehling のポテンシャル
後方でのクーロン散乱断面積の増し分は陽子・陽子では の約 0.5 ％重陽子・原子核では約 0.2 - 0.3 ％δV (r) = −2α3πZ ∞××Z1d3x ρ(x)µ¶p112dt t − 1 2 + 4t2te−2ctRR

電荷密度 ρ はZd3x ρ(x) = Z19ここで、R = |~r−~x|c は電子のコンプトン波長 (約 240 fm)の逆数
α は微細構造定数t = 1 で微分が発散

Euler-Mclaurin の和公式を利用

積分＝台形公式＋(両端での奇数階微分の差)× 荷重

y = log(t − 1) と変数変換
両端点での微分が０2520Coul. pot
208

Pb+d[MeV]

1510Coul. pot
Uehl. pot. x 200

50051015202025R [fm]3035404550

重陽子 folded potential の反対称化
Pong and Austern を参照
核子光学ポテンシャルを利用する。
→ 入射核子と核内核子の反対称化は現象論的に取り込まれている
重陽子内核子と核内核子の反対称化
核内核子： 密度とフェルミガスモデル
重陽子： Hulthén 波動関数
非局所ポテンシャルを導く
局所等価ポテンシャルを作る
局所運動量近似を利用
実部は浅く、虚部は深くなる
中心部には近付くな！！
中心部には来たら消す！！

21

40
real corr. x10

20

imag corr. x100

[MeV]

imag. pot.
-20

208Pb+d, Ed=56 MeV
Anti. symm. correction-40-60real pot.
-80-10002468R [fm]101214２０ MeV 重陽子の鉛核による弾性散乱：
微分断面積の山谷が約５度程度後方にずれて、実験値の再現性が改良

22

(0)

(0)

{E − E (0) − tp − tn − Up − Un }ψ
= (1 − Pp)(1 − Pn)Vpnψ
Singly proj. app:()
24~~ R~ 0) =~U S (R,(β 2 − α2)ρN (R)
M~ −R~ 0|)× H(|RNonlocal kernel:
2αβ(α + β)H(s) =(α − β)2½¾
1 φd(2s)
exp (−2βs)
×−β + α φd(0)2β3j1(kf s)
×kf s∼ A exp (−Bs)∼ AH exp (−BH s2)
kf by local Fermi gas model.
23

S (R):
Equivalent local potential UL
S (R) =UL

4~2(β 2 − α2)Mµ¶3/2π× ρ¯N (R)BH¾½MS − E)(U+U× expdL
B H ~2nucleon density:
208
ρcharge(R)ρ¯N (R) =8224

松岡さんの実験データ
1208

Pb(d,d)σEL/σRuth

Ed= 56 MeV

0.1

0.01

0204060
θCM (deg)

80

100

120

0.8
0.60.4Ay

0.2
0-0.2-0.4
-0.6
-0.8

0204060
θCM (deg)80

緑線はクーロン分解ありの計算
25100120

$d3sigma/dEp,dOmegap,dOmegan (mb/ MeV sr2)$

岡村さんの実験データ
500
450$muI=pm 1$

400
350

$208$Pb(d,pn)
E$d$= 56 MeVunpol. d300$muI=0$250
200
1501005001015202530
35
E$p$ (MeV)40

計算は、偏極ビーム使用。
これをスピン平均すると青線になる

2645501050-5Prot. ang. [deg]15-10-20 -15 -10-50
510Neut. ang. [deg]15

クーロン軌道より、陽子は外へ
中性子は内へ但し、陽子と中性子は同方向へ出る
2720-15両粒子が０度方向に放出される確率は小さい20Prot. [deg]4-2-4-4-202Neut. [deg]4near/far comp. が０度で干渉する
０度での幅と角運動量は不確定性を満足δθ L = 1
L = kR とすると, 分解半径が推定可28

重陽子の電気的分極率~ の下でのRamsay 達は一様電場 E
エネルギー W を考えた

W =−=−αE 22X (0|ezE|n)(n|ezE|0)4(En − E0)n>0
状態 |n) は 3P が主寄与項
８ MeV d+209Bi 弾性散乱
後方で約３％の寄与と評価

29(1)hctak での簡単な計算
1.05
1folding

0.95

dσEL/dσRuth

0.9

full cc

0.85
0.8

209

0.75

BE(d)= 2.0 MeVBi+d, Ed= 8 MeV

0.7

No cont-cont coupling

0.65
0.6
0.55

020406080100
θCM [deg]120140160180

No cont-cont coupling と書いたのは、以下の連立方程式を解いている。

D 2 u0 =X

V0 i uii

D 2uj = Vj j uj + Vj 0u0,30

(j 6= 0)

個人的な考え
Ramsay の評価は、クーロン分解の２次項分解過程はクーロン分解の１次項
従って、弾性散乱からは評価出来ない！
後知恵式 (1) を使用し CDCC のアイデアで直接評価出来そうだ31

OP process
Lawrence 達は、p, α 誘起反応断面積の‘入射エネルギー依存性’ は
Gamow 理論で OK
d 誘起反応断面積はガモフ理論よりも、穏やかなエネルギー依存性

Oppenheimer と Philips はクーロン分解後の中性子吸収として、(d,p) 反応断面積の Ed 依存性を説明

32

OP proc. への CDCC approach
CDCC 波動関数から、重心運動の流れの密度を大きな球面で立体角積分する。
分解チャネルに対し
c J |2 ∝|SLL0Z×Xc0 L0³dR Im0J
∗cJcχL L0 Vc c0 χL0 L0

被積分関数は、湧き口密度を与える
右辺対角要素は虚部ポテンシャルの対角要素だから負吸い込み過程を表す
虚部の中性子 (陽子) 成分は分解後の中性子 (陽子) 吸収断面積σ(d, p) 、σ(d, n)33´

118Sn+d の計算値
1.0e+01
1.0e+00
Nuc. + Coul. BU

1.0e-01

Xsec(au)

1.0e-02118

1.0e-031.0e-04

pure Coul. BU

1.0e-05

Sn+dneutron absorption

pure Nuc. BU

1.0e-06
1.0e-071.0e-0834567
Ed (MeV)8

核力だけの分解過程よりも、かなり穏やかなエネルギー依存性

34910CDCC(d,p)
反応の T 行列を次式と近似する。
(−)

Tβ α =< χb(+)ΦB | Vb − Ub|ΨaΦA >広義重陽子の内部構造
nnnn
lΦn(l sd Ia ma ; r) = i Yln sd Ian,mna φn(r)システム１の波動関数< IB MB |IA MA >X (1)=SB A l j ul1 j (r1)l1 j1

× (IA MA j MB − MA| IB MB )× i−l1 [Yl1 (r̂1) ⊗ φsx ]∗j M −M
B35A

システム２関連の成分
Vb − Ub ∼ Vnp と近似する
n s I n , mn ; r )Vnp(r2) Φn(ld aa 2X a=Dln2 (r2) il2 Yl2 sd Ian,mna (r̂2)
l2

Dln2 (r2) ≡Xlnφln Ian (r2)nll2
× < i Yl2 sd Ian |Vnp|i Yln sd Ian >l2 は重陽子 D 状態及び Vnp のテンソル力に起因< sb mb|Vnp(r2)|Ian mna >X (2)=Sab l s(sb mb s ms| Ian mna)
2l2 s× Dln2 (r2) il2 [Yl2 (r̂2) ⊗ φsx ]s ms

この式で移行スピン s が定義された。
36

核行列要素DWBA の標準的な行列要素に持ち込む。

< IB MB , sb mb|Vnp|IA MA, Ian mna >X X (1)(2)=
SB A l j Sa b l s1l s j l1 l 2

2× (幾何学的因子)
l jl s× i−l fl 1m 2 (r1, r2)
l1 j l 2 s
fl m(r1, r2) = ul1 j (r1) Dln2 (r2)× [Yl1 (r̂1) ⊗ Yl2 (r̂2)]∗l ml jl s
fl 1m 2 は束縛状態 ul1 j (r1) と Dln (r2) の積で、2有効 range がきまる

37

図は、Dln (r2) を色々な” 重陽子 “状態に2対して描いたものである。

Vvp φ

1.0E+02

Vn,n φ Ed=10 MeV

1.0E+01

bound state: L=0

1.0E+00

bound state: L=2

1.0E-01
1.0E-02
1.0E-03
1.0E-040123rn p (fm)456

通常の重陽子の基底状態に対応する場合が大きな振幅を持つ事は明らかである。
いわゆる重陽子 D 状態さえ、大きな寄与はしない。

38

重イオン用の開発 → hctak.f
クーロン波動関数の副プログラム
Virmingham アルゴリズム
３項漸化式の連分数表現

3j, 6j 係数副プログラム
階乗計算の排除
今後の機能追加
１) 区分的摂動近似 (Jost func?)
２) E-CDCC
３) inelastic break up
考えられる利用法
１) 複合核形成断面積の計算
２) 分極率３) 触媒核融合 ？

39

クーロン波動関数
これまでは、1/ρ の漸近展開使用が幾らか難しく利用範囲は制限される
３項漸化式

αFl+1 + βFl + γFl−1 = 0
反復利用で、Fl−1 は発散？

⇓Ã

αFl

µ!+β+γ

Fl−1Fl¶=0Fl+1比で書き、連分数計算に持ち込むと安定的に使用出来る場合がある

40

３ J 係数j~3 = j~1 + j~2, m3 = m1 + m2 とする
５自由度があるj1, j2, j3, m3 を固定した３項漸化式
2~3J 定義式に j3 = (j~1 + j~2)2 を作用

3J (m2 − 1) = p 3J (m2)+q 3J (m2 + 1)
j1, j2, m1, m2 を固定した３項漸化式
3J 定義式に j1z を作用
3J (j3 − 1) = a 3J (j3) + b 3J (j3 + 1)

41

６ J 係数

Elliot-Biedenharn の恒等式()()

j1 j2 j3
j4 j5 j6

j1 j2 j3
j7 j8 j9
P
X=(2J + 1) (−1)J + i ji
J(×)(? ? ?
? ? ?

)(? ? ?
? ? ?

)
? ? ?
? ? ?

j7 = j2, j8 = j1, j9 = 1 を代入
変数に１を含む場合は、顕な式を使用

j6 をパラメータとする、３項漸化式
6J {j6 + 1} = α 6J {j6}+β 6J {j6 − 1}42

漸化式の利用法
3J (j3) を例にとる
先ず jmax に適用。
3J (jmax + 1) = 0 だから、3J (jmax − 1) = C1 × 3J (jmax)
次に、jmax − 1 に適用。etc.
規格化の式

X

(2j + 1) 3J (j)2 = 1

比例係数 Cn の平方和が決まる

3J (jmax) の絶対値は確定
Condon-Shortley の位相規約直交性を用いて、精度チェック
43触媒核融合分解反応の逆過程を取り上げる
粒子１と２が、標的３の近くに来る粒子３を媒介して、１と２が融合する
見掛けの三体力が発生している

γ 線を出さずに、粒子３の運動で保存則を保証する
散乱状態にある粒子１・２が始状態
１・２の束縛状態が終状態
この遷移の S 行列は CDCC で計算可

44弱結合 CDCC 方程式の解法
R. G. Gordon, J. Chem. Phys. 51 14(1969)

d2χc
dr 2

= Vc c χc +X

Vc c0 χc0(2)c06=c

区間 [r0, r1] での解は既知
区間 [r1, r2] での解は？
補助関数を導入する

d2
dr 2

Ã!

FcGc

Ã
= Vc c!
Fc
Gc(3)

Wronskian は一定値 (1)W (r) ≡ Fc0 Gc − Fc G0c = 1(4)(2) 式の解を Fc, Gc で展開
χc(r) = αc(r) Fc(r) + βc(r) Gc(r)(5)αc, βc は、穏やかに変化する r の関数45

Fc(r) α0c(r)+Gc(r) βc0 (r) = 0
Fc0(r) α0c(r)+G0c(r) βc0 (r)X=
Vc c0 (r) χc0 (r)c06= cα0, β 0 について解き、これらを区間 [r1, r2] で積分する
必要に応じ、次の区間へ適用する
激しく振動する χc(r) が、穏やかに変化する α(r), β(r) に変換
閉チャネルにも、適用可能

46