Raynal-Revai coefficients for a Three-body   problem

By N. Shulgina (shulgina@fy.chalmers.se)

This Notebook calculates Raynal-Revai coefficients, which allows to express hyperspherical harmonics (HH), determined in one reference Jacobi system in terms of HH into another Jacobi system and vise versa. The reference system is the Jacobi system in which the dynamical problem should be solved. This system is defined by the following equations:
eqq:={(Subscript[x, 1,k])==(Sqrt[(((Subscript[m, 1])*(Subscript[m, 2]))/((Subscript[m, 1])+(Subscript[m, 2])))])*((Subscript[r, 1,k])-(Subscript[r, 2,k])),
(Subscript[x, 2,k])==(Sqrt[(((Subscript[m, 3])*((Subscript[m, 1])+(Subscript[m, 2])))/((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3])))])*((Subscript[r, 3,k])-(((Subscript[m, 1])*(Subscript[r, 1,k])+(Subscript[m, 2])*(Subscript[r, 2,k]))/((Subscript[m, 1])+(Subscript[m, 2])))),(((Subscript[m, 1])*(Subscript[r, 1,k])+(Subscript[m, 2])*(Subscript[r, 2,k])+(Subscript[m, 3])*(Subscript[r, 3,k]))/((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3])))==(Subscript[rcm, k])}

Input

Give here the input values, the value of the hypermomentum, and the total angular momentum L value, m1in, m2in and m3in are the mass numbers for the two particles and the core:

Kmax=4;
Lmax=0;
m1in=1;
m2in=1;
m3in=8;

Clear[hh3coupledx,hh3coupled]

Rules

Here the  reference Jacobi coordinates (x) are defined through the particle coordinates (r) with the transformation from x to r. The r coordinates have two indexes, the first one defines the the particle and, the second one its the cartesian projection (1= x, 2 = y, 3 = z)

rulextor=Flatten[Table[{(Subscript[x, 1,k])->(Sqrt[(((Subscript[m, 1])*(Subscript[m, 2]))/((Subscript[m, 1])+(Subscript[m, 2])))])*((Subscript[r, 1,k])-(Subscript[r, 2,k])),
(Subscript[x, 2,k])->(Sqrt[(((Subscript[m, 3])*((Subscript[m, 1])+(Subscript[m, 2])))/((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3])))])*((Subscript[r, 3,k])-(((Subscript[m, 1])*(Subscript[r, 1,k])+(Subscript[m, 2])*(Subscript[r, 2,k]))/((Subscript[m, 1])+(Subscript[m, 2]))))},{k,1,3}]]

Here the  particle coordinates (r) are defined in terms of Jacobi coordinates (x) and the center of mass coordinate rcm - r to x

eqq:={(Subscript[x, 1,k])==(Sqrt[(((Subscript[m, 1])*(Subscript[m, 2]))/((Subscript[m, 1])+(Subscript[m, 2])))])*((Subscript[r, 1,k])-(Subscript[r, 2,k])),
(Subscript[x, 2,k])==(Sqrt[(((Subscript[m, 3])*((Subscript[m, 1])+(Subscript[m, 2])))/((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3])))])*((Subscript[r, 3,k])-(((Subscript[m, 1])*(Subscript[r, 1,k])+(Subscript[m, 2])*(Subscript[r, 2,k]))/((Subscript[m, 1])+(Subscript[m, 2])))),(((Subscript[m, 1])*(Subscript[r, 1,k])+(Subscript[m, 2])*(Subscript[r, 2,k])+(Subscript[m, 3])*(Subscript[r, 3,k]))/((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3])))==(Subscript[rcm, k])}

solr=Solve[eqq,{(Subscript[r, 1,k]),(Subscript[r, 2,k]),(Subscript[r, 3,k])}]//Simplify

If the particle coordinates are calculated relatively to center of mass

(Subscript[r, 3,k])- (Subscript[rcm, k])/.solr//Simplify

the rule takes the form

rulertox=Flatten[Table[{{(Subscript[r, 3,k])->(((Sqrt[((((Subscript[m, 1])+(Subscript[m, 2])) (Subscript[m, 3]))/((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3])))]) (Subscript[x, 2,k]))/(Subscript[m, 3])),(Subscript[r, 1,k])->(((Subscript[m, 2]) (Sqrt[(((Subscript[m, 1]) (Subscript[m, 2]))/((Subscript[m, 1])+(Subscript[m, 2])))]) (Subscript[x, 1,k])+(Subscript[m, 1]) ((Sqrt[(((Subscript[m, 1]) (Subscript[m, 2]))/((Subscript[m, 1])+(Subscript[m, 2])))]) (Subscript[x, 1,k])-(Sqrt[((((Subscript[m, 1])+(Subscript[m, 2])) (Subscript[m, 3]))/((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3])))]) (Subscript[x, 2,k])))/((Subscript[m, 1]) ((Subscript[m, 1])+(Subscript[m, 2])))),(Subscript[r, 2,k])->-(((Subscript[m, 1]) (Sqrt[(((Subscript[m, 1]) (Subscript[m, 2]))/((Subscript[m, 1])+(Subscript[m, 2])))]) (Subscript[x, 1,k])+(Subscript[m, 2]) ((Sqrt[(((Subscript[m, 1]) (Subscript[m, 2]))/((Subscript[m, 1])+(Subscript[m, 2])))]) (Subscript[x, 1,k])+(Sqrt[((((Subscript[m, 1])+(Subscript[m, 2])) (Subscript[m, 3]))/((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3])))]) (Subscript[x, 2,k])))/((Subscript[m, 2]) ((Subscript[m, 1])+(Subscript[m, 2]))))}},{k,1,3}]]

Transformation to other Jacobi systems can be done by permutations of the particle coordinates (r) in the reference system. As usual even permutations P123 and P132 are applied. However, odd permutations P13 and P23 can also be applied, and in this case appropriate chages in the rules have to be done.

rulepijk[i_,j_,l_]:={(Subscript[r, i,k_])->(Subscript[r, j,k]),(Subscript[r, j,k_])->(Subscript[r, l,k]),(Subscript[r, l,k_])->(Subscript[r, i,k]),(Subscript[m, i])->(Subscript[m, j]),(Subscript[m, j])->(Subscript[m, l]),(Subscript[m, l])->(Subscript[m, i])}

rulepij[i_,j_]:={(Subscript[r, i,k_])->(Subscript[r, j,k]),(Subscript[r, j,k_])->(Subscript[r, i,k]),(Subscript[m, i])->(Subscript[m, j]),(Subscript[m, j])->(Subscript[m, i])}

Here an illustration of how different permutations acts at on the Jacibi coordinates in the reference system is given. Note the different signs of the vectors under different permutations

(Subscript[x, 1,1])/.rulextor/.rulepijk[1,2,3]
(Subscript[x, 2,1])/.rulextor/.rulepijk[1,2,3]

(Subscript[x, 1,1])/.rulextor/.rulepijk[1,3,2]
(Subscript[x, 2,1])/.rulextor/.rulepijk[1,3,2]

(Subscript[x, 1,1])/.rulextor/.rulepij[1,3]
(Subscript[x, 2,1])/.rulextor/.rulepij[1,3]

(Subscript[x, 1,1])/.rulextor/.rulepij[2,3]
(Subscript[x, 2,1])/.rulextor/.rulepij[2,3]

This rule just substitutes real input masses

rulem:={(Subscript[m, 1])->m1in,(Subscript[m, 2])->m2in,(Subscript[m, 3])->m3in}

Below is quantum number generator, which generates complete HH basis for fixed K and L

Quantum numbers for fixed hypermoment K1 and fixed total  angular moment Lf. The Raynal Revai coefficient do not depend on total Lf projection (MLf) and therefore MLf is set to zero here.

qt[k1_,lf_]:=DeleteCases[Flatten[Table[If[k1≥l1+l2&&EvenQ[k1-l1-l2]&&lf≤l1+l2&&lf≥Abs[l1-l2],hh3coupledx[k1,l1,l2,lf,0, (Subscript[x, 1,1]), (Subscript[x, 1,2]), (Subscript[x, 1,3]), (Subscript[x, 2,1]), (Subscript[x, 2,2]), (Subscript[x, 2,3])],0],{l2,0,k1},{l1,0,k1}]],0]

qtaref[k1_,lf_]:=DeleteCases[Flatten[Table[If[k1≥l1+l2&&EvenQ[k1-l1-l2]&&lf≤l1+l2&&lf≥Abs[l1-l2],hh3coupledref[k1,l1,l2,lf,0, α,θ1,φ1, θ2, φ2],0],{l2,0,k1},{l1,0,k1}]],0]

qtarot[k1_,lf_]:=DeleteCases[Flatten[Table[If[k1≥l1+l2&&EvenQ[k1-l1-l2]&&lf≤l1+l2&&lf≥Abs[l1-l2],hh3coupledrot[k1,l1,l2,lf,0, α,θ1,φ1, θ2, φ2],0],{l2,0,k1},{l1,0,k1}]],0]

Here are lists of hyperharmonics with fixed K and L in A=3

listhh2lx=qt[Kmax,Lmax]

lhh=Length[listhh2lx]

listhh2laref=qtaref[Kmax,Lmax]

listhh2larot=qtarot[Kmax,Lmax]

mat=Table[c[i,j],{i,1,lhh},{j,1,lhh}]

Below is a represention of coupled HH in terms of Jacobi coordinates (x1 and x2)

myspher[l1_,m1_,x1x_,x1y_,x1z_]:=SphericalHarmonicY[l1,m1,θ,ϕ]/.{θ->ArcCos[x1z/Sqrt[x1x^2+x1y^2+x1z^2]],Exp[I*m1*ϕ]->(((x1x+Sign[m1]*I*x1y)/Sqrt[x1x^2+x1y^2+x1z^2])/(Sqrt[1-x1z^2/(x1x^2+x1y^2+x1z^2)]))^Abs[m1]}

myclebs[j1_,mj1_,j2_,mj2_,j_,mj_]:=If[(j1+j2)≥j&&Abs[j1-j2]≤j&&(mj1+mj2)==mj&&Abs[mj]≤j,ClebschGordan[{j1,mj1},{j2,mj2},{j,mj}],0]

sphcoupled2x[l1_,l2_,l_,ml_,x1x_,x1y_,x1z_,x2x_,x2y_,x2z_]:=Sum[myspher[l1,ml1,x1x,x1y,x1z]*
myspher[l2,ml2,x2x,x2y,x2z]*myclebs[l1,ml1,l2,ml2,l,ml],{ml2,-l2,l2},{ml1,-l1,l1}]
hh3coupledx[k1_,l1_,l2_,l_,ml_,x1x_,x1y_,x1z_,x2x_,x2y_,x2z_]:=√((2((k1-l2-l1)/2)!(k1+2)((k1+l2+l1)/2+1)!)/(Gamma[(k1-l2+l1+3)/2]Gamma[(k1+l2-l1+3)/2]))*Sqrt[(x1x^2+x1y^2+x1z^2)/(x1x^2+x1y^2+x1z^2+x2x^2+x2y^2+x2z^2)]^l1*Sqrt[(x2x^2+x2y^2+x2z^2)/(x1x^2+x1y^2+x1z^2+x2x^2+x2y^2+x2z^2)]^l2*JacobiP[(k1-l1-l2)/2,l1+1/2,l2+1/2,(x2x^2+x2y^2+x2z^2-x1x^2-x1y^2-x1z^2)/(x1x^2+x1y^2+x1z^2+x2x^2+x2y^2+x2z^2)]*sphcoupled2x[l1,l2,l,ml,x1x,x1y,x1z,x2x,x2y,x2z]*Sqrt[((x1x^2)+(x1y^2)+(x1z^2)+(x2x^2)+(x2y^2)+(x2z^2))]^k1

sphcoupled2[l1_,l2_,l_,ml_,teta1_,phi1_, teta2_,phi2_]:=Sum[SphericalHarmonicY[l1,ml1,teta1,phi1]*
SphericalHarmonicY[l2,ml2,teta2,phi2]*myclebs[l1,ml1,l2,ml2,l,ml],{ml2,-l2,l2},{ml1,-l1,l1}]
sphcoupled2con[l1_,l2_,l_,ml_,teta1_,phi1_, teta2_,phi2_]:=Sum[SphericalHarmonicY[l1,-ml1,teta1,phi1]*(-1)^(ml1+ml2)*
SphericalHarmonicY[l2,-ml2,teta2,phi2]*myclebs[l1,ml1,l2,ml2,l,ml],{ml2,-l2,l2},{ml1,-l1,l1}]
hh3coupled[k1_,l1_,l2_,l_,ml_,al1_,teta1_,phi1_,teta2_,phi2_]:=√((2((k1-l2-l1)/2)!(k1+2)((k1+l2+l1)/2+1)!)/(Gamma[(k1-l2+l1+3)/2]Gamma[(k1+l2-l1+3)/2]))*Sin[al1]^l1*Cos[al1]^l2*JacobiP[(k1-l1-l2)/2,l1+1/2,l2+1/2,Cos[2*al1]]*sphcoupled2[l1,l2,l,ml,teta1,phi1, teta2,phi2]
hh3coupledcon[k1_,l1_,l2_,l_,ml_,al1_,teta1_,phi1_,teta2_,phi2_]:=√((2((k1-l2-l1)/2)!(k1+2)((k1+l2+l1)/2+1)!)/(Gamma[(k1-l2+l1+3)/2]Gamma[(k1+l2-l1+3)/2]))*Sin[al1]^l1*Cos[al1]^l2*JacobiP[(k1-l1-l2)/2,l1+1/2,l2+1/2,Cos[al1]^2-Sin[al1]^2]*sphcoupled2con[l1,l2,l,ml,teta1,phi1, teta2,phi2]

Raynal-Revai Coefficients Calculation

vhod2=Together[listhh2lx];
vhod1=Simplify[vhod2,(Subscript[x, 2,1])^2+ (Subscript[x, 2,2])^2+ (Subscript[x, 2,3])^2>0&&(Subscript[x, 1,1])^2+ (Subscript[x, 1,2])^2+ (Subscript[x, 1,3])^2>0&&(Subscript[x, 2,1])^2+ (Subscript[x, 2,2])^2+ (Subscript[x, 2,3])^2+(Subscript[x, 1,1])^2+ (Subscript[x, 1,2])^2+ (Subscript[x, 1,3])^2>0]

Direct transformation (function in rotated system expressed in terms of functions  in reference system)

Clear[ev,p123,p132,p13,p23]

ev[i_]:=Extract[vhod1,i]

p13[i_]:=ev[i]/.rulextor/.rulepij[1,3]/.rulertox/.rulem//Simplify

p23[i_]:=ev[i]/.rulextor/.rulepij[2,3]/.rulertox/.rulem//Simplify

p123[i_]:=ev[i]/.rulextor/.rulepijk[1,2,3]/.rulertox/.rulem//Simplify

p132[i_]:=ev[i]/.rulextor/.rulepijk[1,3,2]/.rulertox/.rulem//Simplify

eq13:=Flatten[Table[p13[i]==Sum[c[i,j]*ev[j],{j,1,lhh}],{i,1,lhh}]]

eq23:=Flatten[Table[p23[i]==Sum[c[i,j]*ev[j],{j,1,lhh}],{i,1,lhh}]]

eq123:=Flatten[Table[p123[i]==Sum[c[i,j]*ev[j],{j,1,lhh}],{i,1,lhh}]]

eq132:=Flatten[Table[p132[i]==Sum[c[i,j]*ev[j],{j,1,lhh}],{i,1,lhh}]]

sol13=Flatten[SolveAlways[eq13,{(Subscript[x, 1,1]),(Subscript[x, 1,2]),(Subscript[x, 1,3]),(Subscript[x, 2,1]),(Subscript[x, 2,2]),(Subscript[x, 2,3])}]]

sol23=Flatten[SolveAlways[eq23,{(Subscript[x, 1,1]),(Subscript[x, 1,2]),(Subscript[x, 1,3]),(Subscript[x, 2,1]),(Subscript[x, 2,2]),(Subscript[x, 2,3])}]]

sol123=Flatten[SolveAlways[eq123,{(Subscript[x, 1,1]),(Subscript[x, 1,2]),(Subscript[x, 1,3]),(Subscript[x, 2,1]),(Subscript[x, 2,2]),(Subscript[x, 2,3])}]]

sol132=Flatten[SolveAlways[eq132,{(Subscript[x, 1,1]),(Subscript[x, 1,2]),(Subscript[x, 1,3]),(Subscript[x, 2,1]),(Subscript[x, 2,2]),(Subscript[x, 2,3])}]]

Inverse transformation (function in reference system experssed in terms of functions in rotated system)

eq13inv:=Flatten[Table[ev[i]==Sum[c[i,j]*p13[j],{j,1,lhh}],{i,1,lhh}]]

eq23inv:=Flatten[Table[ev[i]==Sum[c[i,j]*p23[j],{j,1,lhh}],{i,1,lhh}]]

eq123inv:=Flatten[Table[ev[i]==Sum[c[i,j]*p123[j],{j,1,lhh}],{i,1,lhh}]]

eq132inv:=Flatten[Table[ev[i]==Sum[c[i,j]*p132[j],{j,1,lhh}],{i,1,lhh}]]

sol13inv=Flatten[SolveAlways[eq13inv,{(Subscript[x, 1,1]),(Subscript[x, 1,2]),(Subscript[x, 1,3]),(Subscript[x, 2,1]),(Subscript[x, 2,2]),(Subscript[x, 2,3])}]]

sol23inv=Flatten[SolveAlways[eq23inv,{(Subscript[x, 1,1]),(Subscript[x, 1,2]),(Subscript[x, 1,3]),(Subscript[x, 2,1]),(Subscript[x, 2,2]),(Subscript[x, 2,3])}]]

sol123inv=Flatten[SolveAlways[eq123inv,{(Subscript[x, 1,1]),(Subscript[x, 1,2]),(Subscript[x, 1,3]),(Subscript[x, 2,1]),(Subscript[x, 2,2]),(Subscript[x, 2,3])}]]

sol132inv=Flatten[SolveAlways[eq132inv,{(Subscript[x, 1,1]),(Subscript[x, 1,2]),(Subscript[x, 1,3]),(Subscript[x, 2,1]),(Subscript[x, 2,2]),(Subscript[x, 2,3])}]]

The solution has been already obtained by now. The values of C[i,j] are the Raynal Revai coefficients. Below follows solution in various forms, suitable for different applications

Clear [hh3coupledx,hh3coupled]

SolR13=mat.listhh2laref/.sol13

SolR23=mat.listhh2laref/.sol23

SolR123=mat.listhh2laref/.sol123

SolR132=mat.listhh2laref/.sol132

SolR13inv=mat.listhh2larot/.sol13inv

SolR23inv=mat.listhh2larot/.sol23inv

SolR123inv=mat.listhh2larot/.sol123inv

SolR132inv=mat.listhh2larot/.sol132inv

Outputs

As rules
HHrot[i]->∑Cdirect[i,j]HHref[j] -direct transformation
HHref[i]->∑Cinverse[i,j] HHrot[j] -inverse transformation

Direct

Rule13rot=Table[Extract[listhh2larot,i]->Extract[SolR13,i],{i,1,lhh}]

Rule23rot=Table[Extract[listhh2larot,i]->Extract[SolR23,i],{i,1,lhh}]

Rule123rot=Table[Extract[listhh2larot,i]->Extract[SolR123,i],{i,1,lhh}]

Rule132rot=Table[Extract[listhh2larot,i]->Extract[SolR132,i],{i,1,lhh}]

Inverse

Rule13ref=Table[Extract[listhh2laref,i]->Extract[SolR13inv,i],{i,1,lhh}]

Rule23ref=Table[Extract[listhh2laref,i]->Extract[SolR23inv,i],{i,1,lhh}]

Rule123ref=Table[Extract[listhh2laref,i]->Extract[SolR123inv,i],{i,1,lhh}]

Rule132ref=Table[Extract[listhh2laref,i]->Extract[SolR132inv,i],{i,1,lhh}]

As matrices
HHrot =Mat x HHref - direct transformation
HHref= InverseMat x HHrot - inverse transformation

Direct

Mat13=mat/.sol13//MatrixForm

Mat23=mat/.sol23//MatrixForm

Mat123=mat/.sol123//MatrixForm

Mat132=mat/.sol132//MatrixForm

Inverse

InverseMat13=Inverse[mat/.sol13]//MatrixForm

InverseMat23=Inverse[mat/.sol23]//MatrixForm

InverseMat123=Inverse[mat/.sol123]//MatrixForm

InverseMat132=Inverse[mat/.sol132]//MatrixForm

In table form

lcoerr=Partition[DeleteCases[Flatten[Table[If[Kmax≥l1+l2&&EvenQ[Kmax-l1-l2]&&Lmax≤l1+l2&&Lmax≥Abs[l1-l2],{Kmax,Lmax,l1,l2},net],{l2,0,Kmax},{l1,0,Kmax}]],net],4]

from (1-2)-3 to (3-2)-1  by transposition P13

out13d = Partition[Flatten[Table[{"rot", Extract[lcoerr, i], "ref", Extract[lcoerr, j], c[i, j] /. sol13}, {i, 1, lhh}, {j, 1, lhh}]], 11] // TableForm

from (1-2)-3 to (1-3)-2  by transposition P23

out23d = Partition[Flatten[Table[{"rot", Extract[lcoerr, i], "ref", Extract[lcoerr, j], c[i, j] /. sol23}, {i, 1, lhh}, {j, 1, lhh}]], 11] // TableForm

from (1-2)-3 to (2-3)-1  by transposition P123

out123d = Partition[Flatten[Table[{"rot", Extract[lcoerr, i], "ref", Extract[lcoerr, j], c[i, j] /. sol123}, {i, 1, lhh}, {j, 1, lhh}]], 11] // TableForm

from (1-2)-3 to (3-1)-2  by transposition P132

out132d = Partition[Flatten[Table[{"rot", Extract[lcoerr, i], "ref", Extract[lcoerr, j], c[i, j] /. sol132}, {i, 1, lhh}, {j, 1, lhh}]], 11] // TableForm

from (3-2)-1 to (1-2)-3  by transposition P13

out13inv=Partition[Flatten[Table[{"ref",Extract[lcoerr,i],"rot",Extract[lcoerr,j],c[i,j]/.sol13inv},{i,1,lhh},{j,1,lhh}]],11]//TableForm

from (1-3)-2 to (1-2)-3  by transposition P23

out23inv=Partition[Flatten[Table[{"ref",Extract[lcoerr,i],"rot",Extract[lcoerr,j],c[i,j]/.sol23inv},{i,1,lhh},{j,1,lhh}]],11]//TableForm

from (2-3)-1 to (1-2)-3  by transposition P123

out123inv=Partition[Flatten[Table[{"ref",Extract[lcoerr,i],"rot",Extract[lcoerr,j],c[i,j]/.sol123inv},{i,1,lhh},{j,1,lhh}]],11]//TableForm

from (3-1)-2 to (1-2)-3  by transposition P132

out132inv=Partition[Flatten[Table[{"ref",Extract[lcoerr,i],"rot",Extract[lcoerr,j],c[i,j]/.sol132inv},{i,1,lhh},{j,1,lhh}]],11]//TableForm