Harmonics

Tensor harmonics on the 2-sphere

David Brizuela, José M. Martín-García and Guillermo A. Mena Marugán
Instituto de Estructura de la Materia, CSIC, Madrid

(c) 2005-2020, under GPL

jose@xact.es

http://www.xact.es/

http://groups.google.com/group/xAct

October 2005 -- May 2006.
August 2007: Notebook adapted for version 0.9.2 of xAct.
November 2007: Notebook adapted for version 0.9.3 of xAct.
July 2009: Notebook adapted for version 0.9.9 of xAct.

This notebook defines different forms of the tensor harmonics on the sphere S^2. It requires the package xTensor for abstract computations and its twin package xCoba for component computations.

The package Harmonics is distributed under the GNU GPL License

Date[]

xAct`Harmonics`$Version={"0.7.2",{2014,2,15}}

1. Initialization

1.1. GPL

(* Harmonics, a free package for tensor computations with tensor harmonics in Mathematica *)

(* Copyright (C) 2005-2020 David Brizuela, Jose M. Martin-Garcia and Guillermo A. Mena Marugan *)

(* This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License,or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place-Suite 330, Boston, MA 02111-1307, USA. 
*)

1.2. Info package

(* :Title: Harmonics *)

(* :Author: David Brizuela, Jose M. Martin-Garcia and Guillermo A. Mena Marugan *)

(* :Summary: Free package for computations with tensor harmonics *)

(* :Brief Discussion:
   - TO FILL IN
*)
  
(* :Context: xAct`Harmonics` *)

(* :Package Version: 0.7.2 *)

(* :Copyright: David Brizuela, Jose M. Martin-Garcia and Guillermo A. Mena Marugan (2005-2020) *)

(* :History: see Harmonics.History *)

(* :Keywords: *)

(* :Source: Harmonics.nb *)

(* :Warning: *)

(* :Mathematica Version: 5.0 and later *)

(* :Limitations: *)

1.3. BeginPackage

Protect against multiple loading of the package:

With[{xAct`Harmonics`Private`HarmonicsSymbols=DeleteCases[Join[Names["xAct`Harmonics`*"],Names["xAct`Harmonics`Private`*"]],"$Version"|"xAct`Harmonics`$Version"|"$HarmonicsVersionExpected"|"xAct`Harmonics`$xTensorVersionExpected"]},
Unprotect/@xAct`Harmonics`Private`HarmonicsSymbols;
Clear/@xAct`Harmonics`Private`HarmonicsSymbols;
]

If[Unevaluated[xAct`xCore`Private`$LastPackage]===xAct`xCore`Private`$LastPackage,xAct`xCore`Private`$LastPackage="xAct`Harmonics`"];

xAct`xCore`Private`$LastPackage

Explicit (not hidden) import of xCoba, xTensor, xPerm and xCore:

BeginPackage["xAct`Harmonics`",{"xAct`xCoba`","xAct`xTensor`","xAct`xPerm`","xAct`xCore`"}]

Welcome message:

Print[xAct`xCore`Private`bars];
Print["Package xAct`Harmonics`  version ",$Version[[1]],", ",$Version[[2]]];
Print["CopyRight (C) 2005-2020, David Brizuela, Jose M. Martin-Garcia and Guillermo A. Mena Marugan, under the General Public License."];

We specify the context xAct`Harmonics` to avoid overriding the Disclaimer of the other packages. However we need to turn off the message General:shdw temporarily:

Off[General::shdw]
xAct`Harmonics`Disclaimer[]:=Print["These are points 11 and 12 of the General Public License:\n\nBECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM `AS IS´ WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.\n\nIN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES."]
On[General::shdw]

If[xAct`xCore`Private`$LastPackage==="xAct`Harmonics`",
Print[xAct`xCore`Private`bars];
Print["These packages come with ABSOLUTELY NO WARRANTY; for details type Disclaimer[]. This is free software, and you are welcome to redistribute it under certain conditions. See the General Public License for details."];
Print[xAct`xCore`Private`bars]];

Eliminate annoying spell messages (this is now included in xCore):

Off[General::spell]
Off[General::spell1]

1.4. Set up

Screen all dollar indicess:

ReportSet[$PrePrint,ScreenDollarIndices];

Postfix format for derivatives:

ReportSet[$CovDFormat,"Postfix"];

In this notebooks we shall allow metric contractions everywhere:

ReportSetOption[ContractMetric,AllowUpperDerivatives->True];

ReportSetOption[MakeRule,MetricOn->All];

ReportSetOption[MakeRule,ContractMetrics->True];

1.5. Usage messages

This is the only angular function we really need:

calDEdmonds::usage="calDEdmonds[j, mp, m][α, β, γ] computes the Wigner matrix with labels j, mp and m for the SO(3) rotation described by Euler angles α, β, γ, as given by Edmonds. Note that other authors use a different convention for this function.";
calD::usage="calD[LI[l], LI[mp], LI[m]] is the formal scalar denoting the Wigner matrix. It can be transformed into actual functions of the spherical angles using the rule calDRule.";
calDRule::usage="calDRule transforms the calD scalar into the function calDEdmonds with Euler angles α=0, β=θ, γ=ϕ.";

The manifold:

S2::usage="The package Harmonics constructs the 2-sphere S2 with indices IndexList[a, b, c, d, e, f] and ultraindex uS2, and its tangent vbundle TangentS2.";
TangentS2::usage=S2::usage;
a::usage=S2::usage;
b::usage=S2::usage;
c::usage=S2::usage;
d::usage=S2::usage;
e::usage=S2::usage;
f::usage=S2::usage;
uS2::usage=S2::usage;
etaUpS2::usage=S2::usage;
etaDownS2::usage=S2::usage;

The metric on the manifold:

gamma::usage="The package Harmonics defines the metric gamma (prints as γ) on the 2-sphere S2. Its Levi-Civita connection is cd, with associated tensors Riemanncd, Riccicd, RicciScalarcd = 2, and vanishing Torsioncd, Einsteincd and Weylcd. There is also the Christoffelcd connection cd with PD, the antisymmetric tensor epsilongamma, and the absolute value of the determinant AbsDetgamma in the DefaultChart.";
cd::usage=gamma::usage;
Riemanncd::usage=gamma::usage;
Riccicd::usage=gamma::usage;
RicciScalarcd::usage=gamma::usage;
Torsioncd::usage=gamma::usage;
Einsteincd::usage=gamma::usage;
Weylcd::usage=gamma::usage;
Christoffelcd::usage=gamma::usage;
epsilongamma::usage=gamma::usage;
AbsDetgamma::usage=gamma::usage;

The spherical chart:

spherical::usage="The package Harmonics defines the chart spherical on the 2-sphere S2 and assigns values to the components of the metric tensors in this chart. The coordinate fields are called θ[] and ϕ[]. The associated parallel (partial) derivative is called PDspherical. The associated";
θ::usage=spherical::usage;
ϕ::usage=spherical::usage;
PDspherical::usage=spherical::usage;
TorsionPDspherical::usage=spherical::usage;
ChristoffelPDspherical::usage=spherical::usage;
RiemannPDspherical::usage=spherical::usage;
RicciPDspherical::usage=spherical::usage;

Complex basis vectors:

m::usage="m[a] and mbar[a] are the complex Newman-Penrose vectors on the 2-sphere. They are null vectors.";
mbar::usage=m::usage;
mmbarRule::usage="mmbarRule transforms the product of one m[a] and one mbar[b] vectors into a sum of a gamma[a,b] object and epsilongamma[a,b] with suitable coefficients.";

STF harmonics:

Z::usage="Z[LI[l], LI[m], a1, a2, ..., as] represents the polar generalized-GS tensor harmonics with labels l, m, and s indices.";
X::usage="X[LI[l], LI[m], a1, a2, ..., as] represents the axial generalized-GS tensor harmonics with labels l, m, and s indices.";
calZ::usage="calZ[LI[sign], LI[l], LI[m], a1, a2, ..., as] is a combined notation for the generalized-GS tensor harmonics. For sign=1 calZ = Z and for sign=-1 calZ = -I*X.";
ZToX::usage="ZToX and XToZ are rules which transform the Z tensor harmonics into X harmonics and viceversa, respectively.";
XToZ::usage=ZToX::usage;
Dropepsilon::usage="Dropepsilon is a rule which contracts epsilongamma with the tensor harmonics Z and X.";
coeffK::usage="coeffK[l, s] is the normalization factor for pure-spin harmonics.";
PureSpin::usage="PureSpin[l, sign, m][a1, ..., as] returns the pure-spin harmonic as a product of coeffK, calD and a number of complex NP vectors.";
ZXToY::usage="ZXToY is a rule which transforms the Z and X harmonics into PureSpin harmonics.";
Orbitalt::usage="Orbitalt[LI[m], a] are the complex Cartesian vectors. Orbitalt[LI[m], a1, ..., as] are the complex Cartesian tensors obtained by combination of the former.";
PureOrbital::usage="PureOrbital[j, l, m][a1, ..., as] returns the pure-orbital harmonic with labels j, l, m and s indices.";

Computation of components:

HarmonicComponent::usage="HarmonicComponent[expr] returns the scalar expression as a function of the coordinates θ[] and ϕ[].";

Products of harmonics:

coeffE::usage="coeffE[{l1, s1, m1}, {l2, s2, m2}, l] is the formal coefficient E in the expansion of a product of two tensor harmonics. It can be evaluated using the coeffERule.";
coeffERule::usage="coeffERule is a rule which evaluates the numerical value of the E coefficients.";
ProductRule::usage="ProductRule is a rule which expands the product of two tensor harmonics with integer labels.";
FormalProductRule::usage="FormalProductRule[l] is a rule which expresses the product of two tensor harmonics with non-integer labels, using l as last argument of the coefficients E.";

1.6. Begin private

Begin["`Private`"]

Names["xAct`Harmonics`*"]

%//Length

2. Conventions for the angular functions

There are several conventions in the literature for the functions which are used in the description of the rotation group.

The choice of spherical harmonics is now almost universal. From Mathematica's help:

For l≥0, (Subscript[Y, l]^m)(θ,ϕ)=(Sqrt[(2l+1)/(4π)])(Sqrt[(l-m)!/(l+m)!])(Subscript[P, l]^m)(cos(θ))(e^(imϕ)) where (Subscript[P, l]^m) is the associated Legendre function: (Subscript[P, l]^m)(x)=(((-1))^m)/((2^l)l!)(((1-(x^2)))^(m/2))((d^(l+m))/d(x^(l+m)))((((x^2)-1))^l)

Note the phase (((-1))^m) (so-called Condon-Shortley phase) in the Legendre polynomial. Edmonds does not use that phase (Edmonds 2.5.17), but later on he compensates that with the definition of spherical harmonics (Edmonds 2.5.29). Therefore the definitions for (Subscript[Y, l]^m)of Edmonds, Galindo-Pascual and Mathematica all coincide.

For the representation matrices we use Edmonds' convention. Note that α=0, β=θ and γ=ϕ. Note also that his equation (4.1.12) contains a mistake: angles α and γ must be exchanged (this error has been previously pointed our for example in PhysRev A 18 (1978) 392, appendix B):

dEdmonds[j_,mp_,m_][β_]:=0/;Abs[mp]>j||Abs[m]>j
dEdmonds[j_,mp_,m_][β_]:=Sqrt[(j+mp)!(j-mp)!(j+m)!(j-m)!]Sum[(-1)^(j-mp-σ)/((j-mp-σ)!(j-m-σ)!(mp+m+σ)!σ!)Sin[β/2]^(2j-mp-m-2σ)Cos[β/2]^(mp+m+2σ),{σ,Max[-mp-m,0],Min[j-mp,j-m]}]

calDEdmonds[j_,mp_,m_][α_,β_,γ_]:=Exp[I mp α]dEdmonds[j,mp,m][β]Exp[I m γ]

The the formal Wigner matrices calD will be defined below and replaced by their values with the calDRule:

calDRule=calD[LI[l_],LI[mp_],LI[m_]]:>calDEdmonds[l,mp,m][0,θ[],ϕ[]]

dEdmonds[7,3,3][β]//Simplify

Plot[%,{β,0,Pi}]

Galindo-Pascual's definition of the coefficient d coincides with that of Goldberg, but they are different from Edmonds' due to a phase (-1)^(m+mp), or equivalently, they have opposite signs of β. The complex exponentials are also different between Edmonds and Galindo-Pascual and (after correcting the mistake referred to above) we see that the difference is just the exchange of angles α and γ. The exponentials of Goldberg coincide with those of Edmonds before correcting the mistake! Therefore, taking Galindo-Pascual as a reference we find that:
	calDEdmonds(α, β, γ) = calDGalindo(-α, -β, -γ)
	calDGoldberg(α, β, γ) = calDGalindo(-γ, -β, -α) = DGalindo(α, β, γ)^(-1)
That means that the three calD's are different!

d can be related to the Jacobi polynomials (interestingly, timings are very similar).

dJacobi[j_,mp_,m_][β_]:=Sqrt[(j+mp)!(j-mp)!/(j+m)!/(j-m)!]Cos[β/2]^(m+mp)Sin[β/2]^(mp-m)JacobiP[j-mp,mp-m,mp+m,Cos[β]]

dEdmonds[3,1,-3][beta]/dJacobi[3,1,-3][beta]//Simplify

The Clebsch-Gordan coefficients are also standarized. They can be given (and are actually computed by Mathematica in this way) in terms of a hypergeometric function:

CG[{j1_,m1_},{j2_,m2_},{j3_,m3_}]:=(-1)^(j1+m2-j3)W[j3-m3,j2-m2,j1-m1][j1+m1-j2+m2,-j1+m1+j2+m2]

W[m_,n1_,n2_][k1_,k2_]:=(-1)^m (n1+n2)!/Gamma[1+k2]Sqrt[(2m+k1+k2+1)Gamma[m+k1+k2+1]Gamma[m+k2+1]/n1!/n2!/m!/(n1+n2-m)!/Gamma[m+k1+1]]Sqrt[Gamma[n1+k1+1]Gamma[n2+k2+1]/Gamma[n1+n2+m+k1+k2+2]]HypergeometricPFQ[{-m,m+k1+k2+1,-n2},{1+k2,-n1-n2},1]

With[{j1=4,m1=-3,j2=5,m2=2,j3=7,m3=-1},{ClebschGordan[{j1,m1},{j2,m2},{j3,m3}],CG[{j1,m1},{j2,m2},{j3,m3}]}]

This expression allows us to find a limit for large j1+j2:

With[{Lp=4,m1=2,m2=-4,lpp=5,m3=-2,lower=7,upper=1000},nums=Table[ClebschGordan[{(L+Lp+lpp)/2,m1},{(L-Lp+lpp)/2,m2},{lpp,m3}],{L,lower,upper,2}];
approx=Table[(-1)^((L-Lp+lpp)/2+m2)Sqrt[(2lpp+1)/(L+lpp+1)]dJacobi[lpp,Lp,m1+m2][ArcCos[(m1-m2)/(L+lpp+1)]],{L,lower,upper,2}]];

ListPlot[nums/approx]

3. Unit sphere

3.0. Unit sphere

DefUnitSphere[
{manifold_Symbol,dim_Integer,indices_List,manifoldoptions___},
{metric_Symbol,covd_Symbol,metricoptions___},
{chart_Symbol,cnumbers_List,scalars_List,chartoptions___}
]:=Module[{tmp},

(* Checks *)
If[dim=!=2,Throw@Message[]];
If[Length[indices]<3,Throw@Message[]];

With[{a=indices[[1]],b=indices[[2]],c=indices[[3]],theta=scalars[[1]]},

(* Define geometric structures *)
DefManifold[manifold,dim,indices,manifoldoptions];
DefChart[chart,manifold,cnumbers,scalars,chartoptions];
DefMetric[1,metric[-a,-b],covd,metricoptions];
With[{rs=GiveSymbol[RicciScalar,covd]},Set[rs[],2]];

(* Unit metric *)
$CVVerbose=False;
metricvalues=DiagonalMatrix[{1,Sin[theta[]]^2}];
MetricInBasis[metric,-chart,metricvalues];
MetricInBasis[metric,chart,Inverse[metricvalues]];
ComponentValue[Determinant[metric[-a,-b],chart],Sqrt[Det[metricvalues]]//PowerExpand];
Remove[metricvalues];

$BMUseValues=All;
$TUseValues=All;
$CVSimplify=Simplify;
$CCSimplify=Simplify;
With[{eps=GiveSymbol[epsilon,metric],det=GiveSymbol[Det,metric,chart]},
ComponentValue[ComponentArray[eps,{-chart,-chart}],$epsilonSign det[]{{0,1},{-1,0}}//ToValues];
ChangeComponents[eps[{a,chart},{b,chart}],eps[{-a,-chart},{-b,-chart}]]];
AllComponentValues[Christoffel[cd,PDspherical][{-a,-chart},{-b,-chart},{-c,-chart}]];

tmp=metric[{-a,-chart},{-b,-chart}];
tmp=SeparateBasis[AIndex][tmp]-tmp;
tmp=covd[-c][tmp];
tmp=tmp//ContractMetric//ContractBasis;
tmp=tmp//ToBasis[chart];
tmp=ComponentArray[tmp]//ToValues//Flatten;
tmp=Solve[Map[#==0&,tmp],First/@Last[TensorValues[GiveSymbol[Christoffel,covd,GiveSymbol[PD,chart]]]]];
tmp/.Rule->ComponentValue;

ChangeComponents[ChristoffelcdPDspherical[{a,spherical},{-b,-spherical},{-c,-spherical}],ChristoffelcdPDspherical[{-a,-spherical},{-b,-spherical},{-c,-spherical}]];

SetAttributes[identity,Orderless];
identity[a_Symbol,-b_Symbol]:=Basis[a,{1,-spherical}]Basis[-b,{1,spherical}]+Basis[a,{2,-spherical}]Basis[-b,{2,spherical}]+Basis[a,{3,-spherical}]Basis[-b,{3,spherical}];
DefTensor[m[a],S2];
DefTensor[mbar[a],S2,PrintAs->"mb"];
m/:m[a_]m[-a_]:=0;
mbar/:mbar[a_]mbar[-a_]:=0;
mbar/:m[a_]mbar[-a_]:=1;
m/:m[-a_]mbar[a_]:=1;
ComponentValue[ComponentArray[m[{a,spherical}]],{0,1/2,I Csc[θ[]]/2}];
ChangeComponents[m[{-a,-spherical}],m[{a,spherical}]];ComponentValue[ComponentArray[mbar[{a,spherical}]],{0,1/2,-I Csc[θ[]]/2}];
ChangeComponents[mbar[{-a,-spherical}],mbar[{a,spherical}]];
m/:m[a_]Basis[-a_,{b_,spherical}]:=m[{b,spherical}];
m/:m[-a_]Basis[a_,{b_,spherical}]:=m[{b,spherical}];
mbar/:mbar[a_]Basis[-a_,{b_,spherical}]:=mbar[{b,spherical}];
mbar/:mbar[-a_]Basis[a_,{b_,spherical}]:=mbar[{b,spherical}];

]

];

DefHarmonicStructure[]:=Module[{},
(* Cartesian structure *)
DefChart[cartesian,S2,{2,3},{y[],z[]},ChartColor->RGBColor[0,0,1]];
DefTensor[x[],S2];
Coordinate[1,cartesian]:=x;
DefTensor[r[],S2];
Coordinate[1,spherical]:=r;
change={x[]->r[]Sin[θ[]]Cos[ϕ[]],y[]->r[]Sin[θ[]]Sin[ϕ[]],z[]->r[]Cos[θ[]]};
cartesian/:CNumbersOf[cartesian,TangentS2]={1,2,3};
spherical/:CNumbersOf[spherical,TangentS2]={1,2,3};
MatrixForm[changeL=ComponentArray[Basis,{cartesian,-spherical}]];MatrixForm[changeLinv=ComponentArray[Basis,{spherical,-cartesian}]];sphericalcoordinates=Table[Coordinate[i,spherical][],{i,3}];cartesiancoordinates=Table[Coordinate[i,cartesian][],{i,3}]/.change;MatrixForm[changeR=Outer[D,cartesiancoordinates,sphericalcoordinates]];
MatrixForm[changeRinv=Simplify[Inverse[changeR]]];
ComponentValue[changeL,changeR];
ComponentValue[changeLinv,changeRinv];
ChangeComponents[m[{a,cartesian}],m[{a,spherical}]];
ChangeComponents[m[{-a,-cartesian}],m[{-a,-spherical}]];
r[]:=1;

mmbarRule=mbar[a_]m[b_]->(gamma[a,b]+I epsilongamma[a,b])/2;
DefConstantSymbol[l];
DefTensor[Z[LI[l],LI[m],AnyIndices[TangentS2]],S2];
DefTensor[X[LI[l],LI[m],AnyIndices[TangentS2]],S2];
DefTensor[calZ[LI[sign],LI[l],LI[m],AnyIndices[TangentS2]],S2,PrintAs->"ℨ"];
calZ[LI[1],LI[l_],LI[m_]]:=2Z[LI[l],LI[m]];
calZ[LI[1],LI[l_],LI[m_],inds___]:=Z[LI[l],LI[m],inds];
calZ[LI[-1],LI[l_],LI[m_],inds___]:=-I X[LI[l],LI[m],inds];
calZ[LI[c_.(-1)^(a_.+2b_)],LI[l_],LI[m_],inds___]:=calZ[LI[c(-1)^a],LI[l],LI[m],inds];
Z/:Z[LI[l_],LI[m_],inds___]:=0/;Length[{inds}]>l;
X/:X[LI[l_],LI[m_],inds___]:=0/;Length[{inds}]>l;
calZ/:calZ[LI[sign_],LI[l_],LI[m_],inds___]:=0/;Length[{inds}]>l;
SymmetryGroupOfTensor[Z[_LI,_LI,indices___]]^:=Symmetric[2+Range@Length@{indices},Cycles];
SymmetryGroupOfTensor[X[_LI,_LI,indices___]]^:=Symmetric[2+Range@Length@{indices},Cycles];
SymmetryGroupOfTensor[calZ[_LI,_LI,_LI,indices___]]^:=Symmetric[3+Range@Length@{indices},Cycles];
Z[_LI,_LI,___,a_,___,-a_,___]:=0;
Z[_LI,_LI,___,-a_,___,a_,___]:=0;
X[_LI,_LI,___,a_,___,-a_,___]:=0;
X[_LI,_LI,___,-a_,___,a_,___]:=0;
calZ[_LI,_LI,_LI,___,a_,___,-a_,___]:=0;
calZ[_LI,_LI,_LI,___,-a_,___,a_,___]:=0;
X[_LI,_LI]:=0;
X/:Z[LI[l_],LI[m_],inds1__?AIndexQ]X[LI[l_],LI[m_],inds2__?AIndexQ]:=0/;Sort[xAct`xTensor`Private`TakePairs[IndexList[inds1,inds2]]]===Sort[IndexList[inds1,inds2]];
ZToX=Z[l_LI,m_LI,a_,b___]:>Module[{c},-Symmetrize[epsilongamma[a,c]X[l,m,-c,b],{a,b}]];XToZ=X[l_LI,m_LI,a_,b___]:>Module[{c},Symmetrize[epsilongamma[a,c]Z[l,m,-c,b],{a,b}]];calZTocalZ=calZ[LI[sign_],l_LI,m_LI,a_,b___]:>Module[{c},-I Symmetrize[epsilongamma[a,c]calZ[LI[-sign],l,m,-c,b],{a,b}]];
Dropepsilon=term:(_. epsilongamma[a_,b_](Z|X|calZ)[__]):>ReplaceAll[term,{ZToX,XToZ,calZTocalZ}]/;IsIndexOf[term,-a]||IsIndexOf[term,-b];
cd[a_][Z[LI[l_],LI[m_]]]:=Z[LI[l],LI[m],a];
cd[b_][Z[LI[l_],LI[m_],a_?AIndexQ]]:=Z[LI[l],LI[m],a,b]-l(l+1)/2gamma[a,b]Z[LI[l],LI[m]];
cd[b_][X[LI[l_],LI[m_],a_?AIndexQ]]:=X[LI[l],LI[m],a,b]-l(l+1)/2epsilongamma[a,b]Z[LI[l],LI[m]];
cd[b_][Z[LI[l_],LI[m_],a1_?AIndexQ,a2_?AIndexQ,a___?AIndexQ]]:=Module[{b$},With[{s=1+1+Length[{a}]},Z[LI[l],LI[m],a1,a2,a,b$]+(l+s)(l-s+1)/2Symmetrize[1/2gamma[a1,a2]Z[LI[l],LI[m],a,b$]-gamma[b$,a1]Z[LI[l],LI[m],a2,a],{a1,a2,a}]]/.b$->b];
cd[b_][X[LI[l_],LI[m_],a1_?AIndexQ,a2_?AIndexQ,a___?AIndexQ]]:=Module[{b$},With[{s=1+1+Length[{a}]},X[LI[l],LI[m],a1,a2,a,b$]+(l+s)(l-s+1)/2Symmetrize[1/2gamma[a1,a2]X[LI[l],LI[m],a,b$]-gamma[b$,a1]X[LI[l],LI[m],a2,a],{a1,a2,a}]]/.b$->b];
coeffK[l_,s_]:=Infinity/;s>l;
coeffK[l_,s_]:=Sqrt[(2l+1)/(4Pi)/2^s (l+s)!/(l-s)!];

DefTensor[calD[LI[l],LI[mp],LI[m]],manifold,PrintAs->"𝔇"];

PureSpin[l_,1,mlabel_][inds___]:=With[{s=Length[{inds}]},coeffK[l,s]calD[LI@l,LI@s,LI@mlabel](-1)^s Times@@(m/@{inds})];
PureSpin[l_,-1,mlabel_][inds___]:=With[{s=Length[{inds}]},coeffK[l,s]calD[LI@l,LI@-s,LI@mlabel]Times@@(mbar/@{inds})];
ZXToY:={Z[LI[l_],LI[m_]]:>PureSpin[l,1,m][],Z[LI[l_],LI[m_],inds__]:>PureSpin[l,1,m][inds]+PureSpin[l,-1,m][inds],X[LI[l_],LI[m_],inds__]:>I(PureSpin[l,1,m][inds]-PureSpin[l,-1,m][inds])};

spherical/:CNumbersOf[spherical,TangentS2]={2,3};
cartesian/:CNumbersOf[cartesian,TangentS2]={2,3};

HarmonicComponent[expr_]:=Simplify[ToValues[expr/.ZXToY/.calDRule]];

Format[coeffE[{s_,l_,m_},{sp_,lp_,mp_},lpp_]]:=Subsuperscript["E",StringJoin[ToString[s],",",ToString[l],",",ToString[m],";",ToString[lpp]],StringJoin[ToString[sp],",",ToString[lp],",",ToString[mp]]];
coeffE/:FindIndices[coeffE[___]]=IndexList[];
_?FirstDerQ[coeffE[__]]^=0;
IsIndexOf[coeffE[___],__]^:=False;
coeffERule=coeffE[{s_,l_,m_},{sp_,lp_,mp_},lpp_]:>coeffK[lp,Abs[sp]]coeffK[l,Abs[s]]/coeffK[lpp,Abs[s+sp]]ClebschGordan[{lp,mp},{l,m},{lpp,m+mp}]ClebschGordan[{lp,sp},{l,s},{lpp,s+sp}];

TTensor[1][ainds_List,binds_List]:=(-1)^Length[ainds]Times@@(mbar/@ainds)Times@@(m/@binds);TTensor[-1][ainds_List,binds_List]:=(-1)^Length[ainds]Times@@(m/@ainds)Times@@(mbar/@binds);

calTTensor[sign_][ainds_List,binds_List]:=1/2(TTensor[-1][ainds,binds]+sign TTensor[1][ainds,binds]);

ProductRule:={rest_. (tp:(Z|X))[LI[lp_],LI[mp_],a___](t:(Z|X))[LI[l_],LI[mlabel_],b___]:>Module[{sp=Length[{a}],s=Length[{b}],ab,aa,sigma,sigmap,coeff,coeffp},{coeffp,sigmap}=Switch[tp,Z,{1,1},X,{I,-1}];
{coeff,sigma}=Switch[t,Z,{1,1},X,{I,-1}];
ab=Sequence@@Drop[{a},s];
aa=Drop[{a},-sp+s];
rest coeffp coeff If[s===0,1/2,1]If[sp===0,1/2,1]Sum[coeffE[{sp,lp,mp},{s,l,mlabel},lpp]calZ[LI[(-1)^(l+lp-lpp)sigma sigmap],LI[lpp],LI[mp+mlabel],a,b]+sigma coeffE[{sp,lp,mp},{-s,l,mlabel},lpp](calTTensor[+1][aa,{b}]calZ[LI[(-1)^(l+lp-lpp)sigma sigmap],LI[lpp],LI[mp+mlabel],ab]+calTTensor[-1][aa,{b}]calZ[LI[-(-1)^(l+lp-lpp)sigma sigmap],LI[lpp],LI[mp+mlabel],ab]),{lpp,Abs[lp-l],lp+l}]]/;Length[{a}]≥Length[{b}],rest_. Z[LI[l_],LI[mlabel_]]^2:>rest Sum[coeffE[{0,l,mlabel},{0,l,mlabel},lpp]Z[LI[lpp],LI[2mlabel]],{lpp,0,2l}]};

FormalProductRule[lc_]:={rest_. (ta:(Z|X))[LI[la_],LI[ma_],a___](tb:(Z|X))[LI[lb_],LI[mb_],b___]:>Module[{sa=Length[{a}],sb=Length[{b}],sigmaa,sigmab,coeffa,coeffb,indicesA,indicesZ},{coeffa,sigmaa}=Switch[ta,Z,{1,1},X,{I,-1}];
{coeffb,sigmab}=Switch[tb,Z,{1,1},X,{I,-1}];
indicesZ=Sequence@@(Drop[{a},sb]);
indicesA=Drop[{a},sb-sa];
rest coeffa coeffb If[sa===0,1/2,1]If[sb===0,1/2,1](coeffE[{sa,la,ma},{sb,lb,mb},lc]calZ[LI[(-1)^(lb+la-lc)sigmaa sigmab],LI[lc],LI[ma+mb],a,b]+sigmab coeffE[{sa,la,ma},{-sb,lb,mb},lc](calTTensor[+1][indicesA,{b}]calZ[LI[(-1)^(lb+la-lc)sigmaa sigmab],LI[lc],LI[ma+mb],indicesZ]+calTTensor[-1][indicesA,{b}]calZ[LI[-(-1)^(lb+la-lc)sigmaa sigmab],LI[lc],LI[ma+mb],indicesZ]))]/;Length[{a}]≥Length[{b}],rest_. Z[LI[la_],LI[ma_]]^2:>rest coeffE[{0,la,ma},{0,la,ma},lc]calZ[LI[(-1)^(-lc)],LI[lc],LI[2ma]]/2};

]

]

3.1. Definition of manifold

Define S2 manifold

DefManifold[S2,2,{a,b,c,d,e,f}]

with round metric

DefMetric[1,gamma[-a,-b],cd,{":","D"},PrintAs->"γ"]

Constant (unit) Gaussian curvature:

RicciScalarcd[]:=2

We can now define the calD scalars:

DefTensor[calD[LI[l],LI[mp],LI[m]],S2,PrintAs->"𝔇"]

We shall use the convention of having three upper label indices

calD[LI[l],LI[m'],LI[m]]

instead of the usual

calD[LI[l],-LI[m'],-LI[m]]

3.2. Coordinate systems

We use a dirty trick to have 3-dimensional coordinate systems on the 2-dimensional manifold S2:

Define the cartesian coordinate system:

DefChart[cartesian,S2,{2,3},{y[],z[]},ChartColor->RGBColor[0,0,1]]

We enlarge the coordinate system with an additional coordinate:

DefTensor[x[],S2]
Coordinate[1,cartesian]:=x;

Define the spherical system:

DefChart[spherical,S2,{2,3},{θ[],ϕ[]}]

Again, We add a third coordinate to the spherical chart:

DefTensor[r[],S2]
Coordinate[1,spherical]:=r;

Now we give the metric in this coordinate system:

$CVVerbose=False;

metricvalues=DiagonalMatrix[{1,Sin[θ[]]^2}]

MetricInBasis[gamma,-spherical,metricvalues]

MetricInBasis[gamma,spherical,Inverse[metricvalues]]

TensorValues[gamma]

We enlarge the metric:

ComponentValue[{gamma[{1,-spherical},{1,-spherical}],gamma[{1,-spherical},{2,-spherical}],gamma[{1,-spherical},{3,-spherical}],gamma[{2,-spherical},{1,-spherical}],gamma[{3,-spherical},{1,-spherical}]},{0,0,0,0,0}]

The determinant of the metric:

Sqrt[Det[metricvalues]]//PowerExpand

ComponentValue[Determinant[gamma[-a,-b],spherical],%]

Remove[metricvalues];

And the epsilon tensor:

$epsilonSign

ComponentArray[epsilongamma,{-spherical,-spherical}]

$epsilonSign Detgammaspherical[]{{0,1},{-1,0}}//ToValues

ComponentValue[ComponentArray[epsilongamma,{-spherical,-spherical}],$epsilonSign Detgammaspherical[]{{0,1},{-1,0}}//ToValues]

Now we compute the values of epsilongamma in all other index configurations:

$BMUseValues=All;
$TUseValues=All;
$CVSimplify=Simplify;
$CCSimplify=Simplify;

ChangeComponents[epsilongamma[{a,spherical},{b,spherical}],epsilongamma[{-a,-spherical},{-b,-spherical}]]

TensorValues[epsilongamma]

Compute Christoffels:

AllComponentValues[Christoffel[cd,PDspherical][{-a,-spherical},{-b,-spherical},{-c,-spherical}]]

tmp=gamma[{-a,-spherical},{-b,-spherical}]

tmp=SeparateBasis[AIndex][tmp]-tmp

tmp=cd[-e][tmp]

tmp=tmp//ContractMetric//ContractBasis

tmp=tmp//ToBasis[spherical]

tmp=ComponentArray[tmp]//ToValues//Flatten

tmp=Solve[Map[#==0&,tmp],First/@Last[TensorValues[ChristoffelcdPDspherical]]]

tmp/.Rule->ComponentValue

ChangeComponents[ChristoffelcdPDspherical[{a,spherical},{-b,-spherical},{-c,-spherical}],ChristoffelcdPDspherical[{-a,-spherical},{-b,-spherical},{-c,-spherical}]]

Note that the spherical basis is not normalized (because it is a coordinate chart):

gamma[-a,-b]Basis[a,{3,-spherical}]Basis[b,{3,-spherical}]

%//ContractBasis//ToValues

Transformation between spherical and Cartesian charts. Instead of defining a new chart, we define the Cartesian chart as the DefaultChart. Note the trick of using 3 coordinate values rather than 2:

change={x[]->r[]Sin[θ[]]Cos[ϕ[]],y[]->r[]Sin[θ[]]Sin[ϕ[]],z[]->r[]Cos[θ[]]};

spherical/:CNumbersOf[spherical,TangentS2]={1,2,3};
cartesian/:CNumbersOf[cartesian,TangentS2]={1,2,3};

Compute relevant matrices for the change of basis:

MatrixForm[changeL=ComponentArray[Basis,{cartesian,-spherical}]]

MatrixForm[changeLinv=ComponentArray[Basis,{spherical,-cartesian}]]

sphericalcoordinates=Table[Coordinate[i,spherical][],{i,3}]
cartesiancoordinates=Table[Coordinate[i,cartesian][],{i,3}]/.change

MatrixForm[changeR=Outer[D,cartesiancoordinates,sphericalcoordinates]]
MatrixForm[changeRinv=Simplify[Inverse[changeR]]]

Equate them:

ComponentValue[changeL,changeR]

ComponentValue[changeLinv,changeRinv]

r[]:=1

Identity in the spherical basis:

SetAttributes[identity,Orderless]
identity[a_Symbol,-b_Symbol]:=Basis[a,{1,-spherical}]Basis[-b,{1,spherical}]+Basis[a,{2,-spherical}]Basis[-b,{2,spherical}]+Basis[a,{3,-spherical}]Basis[-b,{3,spherical}]

identity[-b,a]

3.3. Complex basis vectors

Then we have the complex basis vectors

DefTensor[m[a],S2]

DefTensor[mbar[a],S2,PrintAs->"(OverBar[m])"]

Orthogonality properties:

m/:m[a_]m[-a_]:=0
mbar/:mbar[a_]mbar[-a_]:=0
mbar/:m[a_]mbar[-a_]:=1
m/:m[-a_]mbar[a_]:=1

Components in the coordinated spherical basis:

ComponentValue[ComponentArray[m[{a,spherical}]],{0,(1/(Sqrt[2])),((I Csc[θ[]])/(Sqrt[2]))}]

ChangeComponents[m[{-a,-spherical}],m[{a,spherical}]]

ComponentValue[ComponentArray[mbar[{a,spherical}]],{0,(1/(Sqrt[2])),-((I Csc[θ[]])/(Sqrt[2]))}]

ChangeComponents[mbar[{-a,-spherical}],mbar[{a,spherical}]]

m/:m[a_]Basis[-a_,{b_,spherical}]:=m[{b,spherical}]
m/:m[-a_]Basis[a_,{b_,spherical}]:=m[{b,spherical}]
mbar/:mbar[a_]Basis[-a_,{b_,spherical}]:=mbar[{b,spherical}]
mbar/:mbar[-a_]Basis[a_,{b_,spherical}]:=mbar[{b,spherical}]

ChangeComponents[m[{a,cartesian}],m[{a,spherical}]]

ChangeComponents[m[{-a,-cartesian}],m[{-a,-spherical}]]

cd[-a][TraceBasisDummy[m[{c,spherical}]Basis[{-c,-spherical},b]]]//ToValues

Basis[a,{3,-spherical}]Basis[-b,{2,spherical}]%//ContractBasis//ToValues

3.4. Products of m, mbar

mmbarRule=mbar[a_]m[b_]->(gamma[a,b]+I epsilongamma[a,b])/2

mbar[a]mbar[c]mbar[e]mbar[f1]

mbar[a]m[c]mbar[e]m[f1]/.mmbarRule

mbar[a]m[c]mbar[e]m[f1]//.mmbarRule

Expand[%]

3.5. Orbital basis vectors

Again, we define the basis tensor with an arbitrary number of indices:

DefTensor[Orbitalt[LI[m],AnyIndices[TangentS2]],S2,PrintAs->"t"]

??Orbitalt

They are STF tensors:

SymmetryGroupOfTensor[Orbitalt[_LI,indices__]]^:=Symmetric[1+Range@Length@{indices}]
Orbitalt[_LI,___,a_,___,-a_,___]:=0
Orbitalt[_LI,___,-a_,___,a_,___]:=0

Define Cartesian components:

AllComponentValues[Orbitalt[LI[1],{a,cartesian}],{-1,-I,0}/Sqrt[2]];
AllComponentValues[Orbitalt[LI[1],{-a,-cartesian}],{-1,-I,0}/Sqrt[2]];
AllComponentValues[Orbitalt[LI[0],{a,cartesian}],{0,0,1}];
AllComponentValues[Orbitalt[LI[0],{-a,-cartesian}],{0,0,1}];
AllComponentValues[Orbitalt[LI[-1],{a,cartesian}],{1,-I,0}/Sqrt[2]];
AllComponentValues[Orbitalt[LI[-1],{-a,-cartesian}],{1,-I,0}/Sqrt[2]];

Components of the vectors in the spherical chart:

AllComponentValues[Orbitalt[LI[1],{a,spherical}],-{Sin[θ[]],Cos[θ[]],I Csc[θ[]]}Exp[I ϕ[]]/Sqrt[2]];
AllComponentValues[Orbitalt[LI[1],{-a,-spherical}],-{Sin[θ[]],Cos[θ[]],I Sin[θ[]]}Exp[I ϕ[]]/Sqrt[2]];
AllComponentValues[Orbitalt[LI[0],{a,spherical}],{Cos[θ[]],-Sin[θ[]],0}];
AllComponentValues[Orbitalt[LI[0],{-a,-spherical}],{Cos[θ[]],-Sin[θ[]],0}];
AllComponentValues[Orbitalt[LI[-1],{a,spherical}],{Sin[θ[]],Cos[θ[]],-I Csc[θ[]]}Exp[-I ϕ[]]/Sqrt[2]];
AllComponentValues[Orbitalt[LI[-1],{-a,-spherical}],{Sin[θ[]],Cos[θ[]],-I Sin[θ[]]}Exp[-I ϕ[]]/Sqrt[2]];

Arbitrary number of indices, in any chart:

Orbitalt[LI[m_Integer],first:{_Integer,basis_},seq:{_Integer,basis_}..]:=Module[{s=Length[{first,seq}],perms=Permutations[{first,seq}]},Sum[ClebschGordan[{1,mp},{s-1,m-mp},{s,m}]Plus@@((Orbitalt[LI[mp],First[#]]Orbitalt[LI[m-mp],Sequence@@Rest[#]])&/@perms)/Length[perms],{mp,-1,1}]]

Off[ClebschGordan::phy]

Orbitalt[LI[2],-a,-b]Basis[b,{2,-spherical}]//ContractBasis

%Basis[a,{3,-spherical}]//ContractBasis

%//ToValues

All components:

Table[Orbitalt[LI[mm],{2,-spherical},{2,-spherical},{2,-spherical}],{mm,-3,3}]

%//ToValues//Simplify

At the poles we have:

northpole={θ[]->0,ϕ[]->0};
southpole={θ[]->Pi,ϕ[]->0};

ToValues[Orbitalt[LI[1],{2,-spherical}]==-m[{2,-spherical}]]/.northpole
ToValues[Orbitalt[LI[1],{3,-spherical}]==-m[{3,-spherical}]]/.northpole

ToValues[Orbitalt[LI[-1],{2,-spherical}]==mbar[{2,-spherical}]]/.northpole
ToValues[Orbitalt[LI[-1],{3,-spherical}]==mbar[{3,-spherical}]]/.northpole

ToValues[Orbitalt[LI[-1],{2,-spherical}]==-m[{2,-spherical}]]/.southpole
ToValues[Orbitalt[LI[-1],{3,-spherical}]==-m[{3,-spherical}]]/.southpole

ToValues[Orbitalt[LI[1],{2,-spherical}]==mbar[{2,-spherical}]]/.southpole
ToValues[Orbitalt[LI[1],{3,-spherical}]==mbar[{3,-spherical}]]/.southpole

These three must coincide. Note the importance of the (-1)^s factor:

{ToValues[Orbitalt[LI[3],{2,-spherical},{2,-spherical},{2,-spherical}]]/.northpole,ToValues[Orbitalt[LI[-3],{2,-spherical},{2,-spherical},{2,-spherical}]]/.southpole,ToValues[(-1)^3m[{2,-spherical}]m[{2,-spherical}]m[{2,-spherical}]]}

We can compute it also using the Cartesian components:

Orbitalt[LI[3],{-a,-cartesian},{-b,-cartesian},{-c,-cartesian}]Basis[{a,cartesian},{2,-spherical}]Basis[{b,cartesian},{2,-spherical}]Basis[{c,cartesian},{2,-spherical}]

TraceBasisDummy[%]

%//ToValues

%//Simplify

%/.northpole

Orbitalt[LI[3],{1,-cartesian},{1,-cartesian},{1,-cartesian}]

%//ToValues

Cut and paste:

((Subscript[t,  1]^(1 ))^3)//InputForm

4. Spherical harmonics

STF Spherical harmonics

Harmonic label:

DefConstantSymbol[l];

Now we define the STF tensors.

DefTensor[Z[LI[l],LI[m],AnyIndices[TangentS2]],S2]

DefTensor[X[LI[l],LI[m],AnyIndices[TangentS2]],S2]

Common notation with a sign. Note the stupid exception for scalars:

DefTensor[calZ[LI[sign],LI[l],LI[m],AnyIndices[TangentS2]],S2,PrintAs->"ℨ"]

calZ[LI[1],LI[l_],LI[m_]]:=2Z[LI[l],LI[m]]
calZ[LI[1],LI[l_],LI[m_],inds___]:=Z[LI[l],LI[m],inds]
calZ[LI[-1],LI[l_],LI[m_],inds___]:=-I X[LI[l],LI[m],inds]

calZ[LI[c_.(-1)^(a_.+2b_)],LI[l_],LI[m_],inds___]:=calZ[LI[c(-1)^a],LI[l],LI[m],inds]

calZ[LI[-(-1)^(x+2y)],LI[l],LI[m],-a,-b]//InputForm

0) We need |m|≤l and s≤l:

Z/:Z[LI[l_],LI[m_],inds___]:=0/;Length[{inds}]>l;
X/:X[LI[l_],LI[m_],inds___]:=0/;Length[{inds}]>l;
calZ/:calZ[LI[sign_],LI[l_],LI[m_],inds___]:=0/;Length[{inds}]>l;

1) We introduce their symmetry properties:

SymmetryGroupOfTensor[Z[_LI,_LI,indices___]]^:=Symmetric[2+Range@Length@{indices},Cycles]

SymmetryGroupOfTensor[X[_LI,_LI,indices___]]^:=Symmetric[2+Range@Length@{indices},Cycles]

SymmetryGroupOfTensor[calZ[_LI,_LI,_LI,indices___]]^:=Symmetric[3+Range@Length@{indices},Cycles]

2) They are all traceless:

Z[_LI,_LI,___,a_,___,-a_,___]:=0
Z[_LI,_LI,___,-a_,___,a_,___]:=0
X[_LI,_LI,___,a_,___,-a_,___]:=0
X[_LI,_LI,___,-a_,___,a_,___]:=0
calZ[_LI,_LI,_LI,___,a_,___,-a_,___]:=0
calZ[_LI,_LI,_LI,___,-a_,___,a_,___]:=0

3) There are no pseudoscalars:

X[_LI,_LI]:=0

4) Orthogonality relations:

X/: Z[LI[l_],LI[m_],inds1__?AIndexQ]X[LI[l_],LI[m_],inds2__?AIndexQ]:=0/;Sort[xAct`xTensor`Private`TakePairs[IndexList[inds1,inds2]]]===Sort[IndexList[inds1,inds2]]

5) Products with epsilongamma:

ZToX=Z[l_LI,m_LI,a_,b___]:>Module[{c},-Symmetrize[epsilongamma[a,c]X[l,m,-c,b],{a,b}]];
XToZ=X[l_LI,m_LI,a_,b___]:>Module[{c},Symmetrize[epsilongamma[a,c]Z[l,m,-c,b],{a,b}]];
calZTocalZ=calZ[LI[sign_],l_LI,m_LI,a_,b___]:>Module[{c},-I Symmetrize[epsilongamma[a,c]calZ[LI[-sign],l,m,-c,b],{a,b}]];

Dropepsilon=term:(_. epsilongamma[a_,b_](Z|X|calZ)[__]):>ReplaceAll[term,{ZToX,XToZ,calZTocalZ}]/;IsIndexOf[term,-a]||IsIndexOf[term,-b];

An alternative coding (not used now):

Dropepsilon=term:(_. epsilongamma[___,a_,___](Z|X|calZ)[___,-a_,___]|_. epsilongamma[___,-a_,___](Z|X)[___,a_,___]):>ReplaceAll[term,{ZToX,XToZ,calZTocalZ}];

6) We define the relations between STF tensors with different numbers of indices. Note that the only place where epsilongamma appears is in the derivative of X_a:

cd[a_][Z[LI[l_],LI[m_]]]:=Z[LI[l],LI[m],a]

cd[b_][Z[LI[l_],LI[m_],a_?AIndexQ]]:=Z[LI[l],LI[m],a,b]-l(l+1)/2gamma[a,b]Z[LI[l],LI[m]]

cd[b_][X[LI[l_],LI[m_],a_?AIndexQ]]:=X[LI[l],LI[m],a,b]-l(l+1)/2epsilongamma[a,b]Z[LI[l],LI[m]]

cd[b_][Z[LI[l_],LI[m_],a1_?AIndexQ,a2_?AIndexQ,a___?AIndexQ]]:=Module[{b$},With[{s=1+1+Length[{a}]},Z[LI[l],LI[m],a1,a2,a,b$]+(l+s)(l-s+1)/2Symmetrize[1/2gamma[a1,a2]Z[LI[l],LI[m],a,b$]-gamma[b$,a1]Z[LI[l],LI[m],a2,a],{a1,a2,a}]]/.b$->b]

cd[b_][X[LI[l_],LI[m_],a1_?AIndexQ,a2_?AIndexQ,a___?AIndexQ]]:=Module[{b$},With[{s=1+1+Length[{a}]},X[LI[l],LI[m],a1,a2,a,b$]+(l+s)(l-s+1)/2Symmetrize[1/2gamma[a1,a2]X[LI[l],LI[m],a,b$]-gamma[b$,a1]X[LI[l],LI[m],a2,a],{a1,a2,a}]]/.b$->b]

Pure-spin harmonics

Normalization factor:

coeffK[l_,s_]:=Infinity/;s>l
coeffK[l_,s_]:=Sqrt[(2l+1)/(4Pi)/2^s (l+s)!/(l-s)!]

Asymptotics for large l+lp:

Series[coeffK[(L-Lp+lpp)/2,s]coeffK[(L+Lp+lpp)/2,sp],{L,Infinity,1}]//FullSimplify

PowerExpand[%]//FullSimplify

Normal[%]//FullSimplify

The pure-spin harmonics are:

PureSpin[l_,1,mlabel_][inds___]:=With[{s=Length[{inds}]},coeffK[l,s]calD[LI@l,LI@s,LI@mlabel](-1)^s Times@@(m/@{inds})]
PureSpin[l_,-1,mlabel_][inds___]:=With[{s=Length[{inds}]},coeffK[l,s]calD[LI@l,LI@-s,LI@mlabel]Times@@(mbar/@{inds})]

Note that both PureSpin[-1] and PureSpin[+1] with no indices are equal to the scalar spherical harmonics.

With[{j=5,m=2},{PureSpin[j,1,m][],PureSpin[j,-1,m][]}-SphericalHarmonicY[j,m,θ[],ϕ[]]/.calDRule]//Simplify

With[{j=4,m=-1},{PureSpin[j,1,m][],PureSpin[j,-1,m][]}-SphericalHarmonicY[j,m,θ[],ϕ[]]/.calDRule]//Simplify

Relations among them

ZXToY:={
Z[LI[l_],LI[m_]]:>PureSpin[l,1,m][],
Z[LI[l_],LI[m_],inds__]:>PureSpin[l,1,m][inds]+PureSpin[l,-1,m][inds],
X[LI[l_],LI[m_],inds__]:>I(PureSpin[l,1,m][inds]-PureSpin[l,-1,m][inds])}

paper1=Z[LI[l],LI[m],-a,-b,-c]-(l^2+l+2)/4gamma[-a,-b]cd[-c][Z[LI[l],LI[m]]]-(l+2)(l-1)/2 Symmetrize[gamma[-c,-a]cd[-b][Z[LI[l],LI[m]]],{-a,-b}]

Simplification[cd[-c]@cd[-b]@cd[-a]@Z[LI[l],LI[m]]-paper1]

paper2=Z[LI[l],LI[m],-a,-b,-c,-d]-(l^2+l+6)/6gamma[-a,-b]cd[-d]@cd[-c]@Z[LI[l],LI[m]]-Symmetrize[
(l+3)(l-2)/6gamma[-c,-d]cd[-b]@cd[-a]@Z[LI[l],LI[m]]+l(l+1)/3gamma[-c,-a]cd[-b]@cd[-d]@Z[LI[l],LI[m]]+ (l+3)(l-2)/3gamma[-d,-a]cd[-c]@cd[-b]@Z[LI[l],LI[m]]+(l+3)(l+1)l(l-2)/12(gamma[-c,-a]gamma[-b,-d]+gamma[-a,-b]gamma[-c,-d]/2)Z[LI[l],LI[m]],{-a,-b}];

Simplification[cd[-d]@cd[-c]@cd[-b]@cd[-a]@Z[LI[l],LI[m]]-paper2]

More checks (this are full expansions, brute force):

STFPart[cd[-c]@cd[-b]@cd[-a]@Z[LI[l],LI[m]],gamma]//ContractMetric

%//Simplification

STFPart[cd[-d]@cd[-c]@cd[-b]@cd[-a]@Z[LI[l],LI[m]],gamma]//ContractMetric

%//Simplification

Explicit expressions

spherical/:CNumbersOf[spherical,TangentS2]={2,3};
cartesian/:CNumbersOf[cartesian,TangentS2]={2,3};

HarmonicComponent[expr_]:=Simplify[ToValues[expr/.ZXToY/.calDRule]]

Components of tensors can be directly obtained:

{Z[LI[4],LI[1],{2,-spherical}],Z[LI[4],LI[1],{3,-spherical}]}

%//HarmonicComponent

{X[LI[4],LI[1],{2,-spherical}],X[LI[4],LI[1],{3,-spherical}]}

%//HarmonicComponent

However, derivatives must be computed with care. For example:

cd[{3,-spherical}][cd[{3,-spherical}][Z[LI[4],LI[1]]]]

%//HarmonicComponent

That is not a component of the second derivative of a scalar harmonic:

ContractBasis[Basis[{3,-spherical},a]Basis[{3,-spherical},b]cd[-a][cd[-b][Z[LI[4],LI[1]]]]]//ToValues

%//HarmonicComponent

The reason is the way xTensor  understands component derivatives:

cd[{3,-spherical}][cd[{3,-spherical}][Z[LI[4],LI[1]]]]

SeparateBasis[AIndex][%,IndicesOf[spherical]]

SeparateBasis[spherical][%]//Expand

%//TraceBasisDummy//ToValues//Expand

%//HarmonicComponent

Or in other words:

HarmonicComponent[Z[LI[4],LI[1]]]

D[%,{ϕ[],2}]//Simplify

ContractBasis[Basis[{2,-spherical},a]Basis[{3,-spherical},b]cd[-a][cd[-b][Z[LI[4],LI[1]]]]]//ToValues

%//HarmonicComponent

cd[-a][Basis[{2,-spherical},b]]

ContractBasis[Basis[{3,-spherical},a]Basis[{3,-spherical},b]cd[-b][X[LI[4],LI[2],-a]]]//HarmonicComponent

ContractBasis[Basis[{3,-spherical},a]Basis[{3,-spherical},b]Basis[{3,-spherical},c]cd[-a]@cd[-b][cd[-c][Z[LI[4],LI[1]]]]]//HarmonicComponent

ContractBasis[Basis[{3,-spherical},a]Basis[{3,-spherical},b]Basis[{3,-spherical},c]cd[-a]@cd[-b][X[LI[4],LI[1],-c]]]//HarmonicComponent

New example:

Z[LI[4],LI[3],a,b]X[LI[4],LI[3],-a,c]

%//SeparateBasis[spherical]

ContractBasis[%,AIndex]

%//HarmonicComponent

%//TraceBasisDummy//ToValues//Simplify

Pure-orbital harmonics

Here we show that the Cartesian components of the pure-orbital harmonics are eigenfunctions of the Laplacian, but not of the z component of angular momentum. Note that this form for the Laplacian is only valid on scalars:

LaplacianOnScalars[expr_]:=D[expr,θ[],θ[]]+Cos[θ[]]/Sin[θ[]]D[expr,θ[]]+D[expr,ϕ[],ϕ[]]/Sin[θ[]]^2

Note that we use Separate twice:

?SeparateBasis

Laplacian[expr_]:=SeparateBasis[spherical][SeparateBasis[spherical][DummyToBasis[spherical][ContractBasis[gamma[a,b]cd[-a][cd[-b][expr]]]],IndicesOf[Free]],IndicesOf[Free]]//TraceBasisDummy//ToValues//Simplify

Definition:

PureOrbital[j_,l_,m_][inds__]:=With[{s=Length[{inds}]},Sum[ClebschGordan[{l,m-mp},{s,mp},{j,m}]SphericalHarmonicY[l,m-mp,θ[],ϕ[]]Orbitalt[LI[mp],inds],{mp,-s,s}]]

PureOrbital[4,3,-2][{3,-cartesian},{3,-cartesian},{2,-cartesian}]//ToValues//Simplify

LaplacianOnScalars[%]/%//Simplify

We have these two orbital harmonics:

O1=With[{j=4,m=-2},Sqrt[j(j+1)]Sqrt[j/(2j+1)]PureOrbital[j,j+1,m][c]identity[d,-c]//ContractBasis//ToValues//Simplify]

O2=With[{j=4,m=-2},Sqrt[j(j+1)]Sqrt[(j+1)/(2j+1)]PureOrbital[j,j-1,m][c]identity[d,-c]//ContractBasis//ToValues//Simplify]

which add up to a spin harmonic (note that the radial part cancels)

newZ=O1+O2//Simplify

newZ-Z[LI[4],LI[-2],c]identity[d,-c]//ContractBasis//HarmonicComponent//Simplify

which has the correct geometric Laplacian:

gamma[a,b]cd[-a]@cd[-b]@Z[LI[4],LI[-2],c]//ContractMetric

gamma[a,b]cd[-a]@cd[-b]@Z[LI[l],LI[m],c]//ContractMetric//Simplify

gamma[a,b]cd[-a]@cd[-b]@Z[LI[l],LI[m],c,d]//ContractMetric//Simplification//Simplify

gamma[a,b]cd[-a]@cd[-b]@Z[LI[l],LI[m],c,d,e]//ContractMetric//Simplification//Simplify

gamma[a,b]cd[-a]@cd[-b]@Z[LI[l],LI[m],c,d,e,f]//ContractMetric//Simplification//Simplify

Laplacian[newZ]

%+19 newZ//Simplify

Let us study the tangential part:

O1tangential=O1/.Basis[_,{1,-spherical}]->0//Simplify

O2tangential=O2/.Basis[_,{1,-spherical}]->0//Simplify

The following expressions have changed and now are far more complicated. In particular the sums should yield zero. It must be related to the issue of whether the cnumbers of the bases are {1,2,3} or {2,3}. Removing the Christoffel terms containing a component 1 and removing the result in the direction 1 we get zero again:

Laplacian[O1tangential]

%+19 O1tangential//Simplify

%/.ChristoffelcdPDspherical[___,{1,_},___]->0

Laplacian[O2tangential]

%+19 O2tangential//Simplify

%/.ChristoffelcdPDspherical[___,{1,_},___]->0

Factor[l(l+1)-s^2]

Now we study the Laplacian acting on the scalar components:

basis[spherical]:={Basis[-d,{1,spherical}], Basis[-d,{2,spherical}],Basis[-d,{3,spherical}]};
basis[cartesian]:={Basis[-d,{1,cartesian}], Basis[-d,{2,cartesian}],Basis[-d,{3,cartesian}]}

O1 basis[spherical]//ContractBasis//Simplify

LaplacianOnScalars[%]/%//Simplify

O2 basis[spherical]//ContractBasis//Simplify

LaplacianOnScalars[%]/%//Simplify

newZ basis[spherical]//ContractBasis//Simplify

LaplacianOnScalars[%]/%//Simplify

Cartesian components:

O1 basis[cartesian]//ContractBasis//ToValues//Simplify

LaplacianOnScalars[%]/%//Simplify

O2 basis[cartesian]//ContractBasis//ToValues//Simplify

LaplacianOnScalars[%]/%//Simplify

newZ basis[cartesian]//ContractBasis//ToValues//Simplify

LaplacianOnScalars[%]/%//Simplify

5. Products of harmonics

Products of calDs

This checks that the formula in the paper is correct:

Off[ClebschGordan::phy]

With[{j1=4,m1p=1,m1=-1,j2=6,m2p=-2,m2=4},calDEdmonds[j1,m1p,m1][alpha,beta,gamma]calDEdmonds[j2,m2p,m2][alpha,beta,gamma]-Sum[ClebschGordan[{j1,m1},{j2,m2},{j,m1+m2}]ClebschGordan[{j1,m1p},{j2,m2p},{j,m1p+m2p}]calDEdmonds[j,m1p+m2p,m1+m2][alpha,beta,gamma],{j,Abs[j1-j2],j1+j2}]]//Simplify

Products of PureSpins

Product coefficients. We do not want to define tensors, but need to state how the main functions act on them:

Format[coeffE[{s_,l_,m_},{sp_,lp_,mp_},lpp_]]:=Subsuperscript["E",StringJoin[ToString[s],",",ToString[l],",",ToString[m],";",ToString[lpp]],StringJoin[ToString[sp],",",ToString[lp],",",ToString[mp]]]

coeffE/:FindIndices[coeffE[___]]=IndexList[];

_?FirstDerQ[coeffE[__]]^=0;

IsIndexOf[coeffE[___],__]^:=False;

Substitution rule:

coeffERule=coeffE[{s_,l_,m_},{sp_,lp_,mp_},lpp_]:>coeffK[lp,Abs[sp]]coeffK[l,Abs[s]]/coeffK[lpp,Abs[s+sp]]ClebschGordan[{lp,mp},{l,m},{lpp,m+mp}]ClebschGordan[{lp,sp},{l,s},{lpp,s+sp}]

Scalar case:

coeffE[{0,l,m},{0,lp,mp},lpp]/(ClebschGordan[{l,0},{lp,0},{lpp,0}]ClebschGordan[{lp,mp},{l,m},{lpp,m+mp}])/.coeffERule//Simplify

Mathematica knows how to expand factorials near Infinity (Stirling expansion)

Series[Gamma[x+1],{x,Infinity,1}]//Simplify

Expansion of the prefactor of the E coefficients for large L:

Series[coeffK[(L+Lp+lpp)/2,sp]coeffK[(L-Lp+lpp)/2,s],{L,Infinity,1}]//Simplify

PowerExpand[%]//Simplify

Checks of the formulas in the paper:

With[{lp=3,mp=2,l=3,mlabel=0},PureSpin[lp,1,mp][a,b,c]PureSpin[l,1,mlabel][d,e,f]-Sum[coeffE[{3,lp,mp},{3,l,mlabel},j]PureSpin[j,1,mp+mlabel][a,b,c,d,e,f],{j,Max[Abs[lp-l],6],lp+l}]]//Simplify

%/.coeffERule

%/.calDRule

With[{lp=3,mp=-2,l=3,mlabel=-1},PureSpin[lp,-1,mp][a,b,c]PureSpin[l,-1,mlabel][d,e,f]-Sum[coeffE[{-3,lp,mp},{-3,l,mlabel},j]PureSpin[j,-1,mp+mlabel][a,b,c,d,e,f],{j,Max[Abs[lp-l],6],lp+l}]]//Simplify

%/.coeffERule

%/.calDRule

Collect[With[{lp=3,mp=-2,l=5,mlabel=-3},PureSpin[lp,-1,mp][a,b]PureSpin[l,1,mlabel][c,d,e,f]-Sum[coeffE[{-2,lp,mp},{4,l,mlabel},j]PureSpin[j,1,mp+mlabel][e,f](-mbar[a]m[c])(-mbar[b]m[d]),{j,Max[Abs[lp-l],0],lp+l}]]//Expand,{epsilongamma[__]},Simplify]

%/.coeffERule//Simplify

%/.calDRule//Simplify

Collect[With[{lp=3,mp=-2,l=5,mlabel=-3},PureSpin[lp,1,mp][a,b]PureSpin[l,-1,mlabel][c,d,e,f]-Sum[coeffE[{2,lp,mp},{-4,l,mlabel},j]PureSpin[j,-1,mp+mlabel][e,f](-m[a]mbar[c])(-m[b]mbar[d]),{j,Max[Abs[lp-l],0],lp+l}]]//Expand,{epsilongamma[__]},Simplify]

%/.coeffERule//Simplify

%/.calDRule//Simplify

Collect[With[{lp=5,mp=-2,l=3,mlabel=-3},PureSpin[lp,-1,mp][a,b,e,f]PureSpin[l,1,mlabel][c,d]-Sum[coeffE[{-4,lp,mp},{2,l,mlabel},j]PureSpin[j,-1,mp+mlabel][e,f](-mbar[a]m[c])(-mbar[b]m[d]),{j,Max[Abs[lp-l],0],lp+l}]]//Expand,{epsilongamma[__]},Simplify]

%/.coeffERule//Simplify

%/.calDRule//Simplify

Products of X and Z

Now we implement the final product formula. We first define the TTensor can calTTensor tensors with even number of indices (note that the a,b indices are reorganized in the paper):

TTensor[1][ainds_List,binds_List]:=(-1)^Length[ainds]Times@@(mbar/@ainds)Times@@(m/@binds);
TTensor[-1][ainds_List,binds_List]:=(-1)^Length[ainds]Times@@(m/@ainds)Times@@(mbar/@binds);

calTTensor[sign_][ainds_List,binds_List]:=(1/2)(TTensor[-1][ainds,binds]+sign  TTensor[1][ainds,binds]);

calTTensor[sign_][ainds_List,binds_List]:=(-1)^Length[ainds](sign Times@@(mbar/@ainds)Times@@(m/@binds)+Times@@(m/@ainds)Times@@(mbar/@binds))/2

calTTensor[-1][{a,c,e},{b,d,f}]

and the rule is:

ProductRule:={rest_. (tp:(Z|X))[LI[lp_],LI[mp_],a___](t:(Z|X))[LI[l_],LI[mlabel_],b___]:>Module[{sp=Length[{a}],s=Length[{b}],ab,aa,sigma,sigmap,coeff,coeffp},
{coeffp,sigmap}=Switch[tp,Z,{1,1},X,{I,-1}];
{coeff,sigma}=Switch[t,Z,{1,1},X,{I,-1}];
ab=Sequence@@Drop[{a},s];
aa=Drop[{a},-sp+s];
rest  coeffp coeff If[s===0,1/2,1]If[sp===0,1/2,1]Sum[
coeffE[{sp,lp,mp},{s,l,mlabel},lpp]calZ[LI[(-1)^(l+lp-lpp)sigma sigmap],LI[lpp],LI[mp+mlabel],a,b]+sigma coeffE[{sp,lp,mp},{-s,l,mlabel},lpp](calTTensor[+1][aa,{b}]calZ[LI[(-1)^(l+lp-lpp)sigma sigmap],LI[lpp],LI[mp+mlabel],ab]+calTTensor[-1][aa,{b}]calZ[LI[-(-1)^(l+lp-lpp)sigma sigmap],LI[lpp],LI[mp+mlabel],ab]),{lpp,Abs[lp-l],lp+l}]]/;Length[{a}]≥Length[{b}],rest_. Z[LI[l_],LI[mlabel_]]^2:>
rest Sum[coeffE[{0,l,mlabel},{0,l,mlabel},lpp]Z[LI[lpp],LI[2mlabel]],{lpp,0,2l}]}

This is the formal (non-expanding) form. It is exactly the previous rule, without the Sum command. Note the special case for a square of a harmonic. Note also that this rule multiplies Z's and X's, but not calZ's.

FormalProductRule[lc_]:={rest_. (ta:(Z|X))[LI[la_],LI[ma_],a___](tb:(Z|X))[LI[lb_],LI[mb_],b___]:>Module[{sa=Length[{a}],sb=Length[{b}],sigmaa,sigmab,coeffa,coeffb,indicesA,indicesZ},
{coeffa,sigmaa}=Switch[ta,Z,{1,1},X,{I,-1}];
{coeffb,sigmab}=Switch[tb,Z,{1,1},X,{I,-1}];
indicesZ=Sequence@@(Drop[{a},sb]);
indicesA=Drop[{a},sb-sa];rest  coeffa coeffb  If[sa===0,1/2,1]If[sb===0,1/2,1](coeffE[{sa,la,ma},{sb,lb,mb},lc]calZ[LI[(-1)^(lb+la-lc)sigmaa sigmab],LI[lc],LI[ma+mb],a,b]+sigmab coeffE[{sa,la,ma},{-sb,lb,mb},lc](calTTensor[+1][indicesA,{b}]calZ[LI[(-1)^(lb+la-lc)sigmaa sigmab],LI[lc],LI[ma+mb],indicesZ]+calTTensor[-1][indicesA,{b}]calZ[LI[-(-1)^(lb+la-lc)sigmaa sigmab],LI[lc],LI[ma+mb],indicesZ]))]/;Length[{a}]≥Length[{b}],rest_. Z[LI[la_],LI[ma_]]^2:>
rest  coeffE[{0,la,ma},{0,la,ma},lc]calZ[LI[(-1)^(-lc)],LI[lc],LI[2ma]]/2}

It is important to note that when both harmonics have the same labels la and ma then the E-coefficients vanish for odd lc, which ranges from 0 to 2la through even integers. In order to make that point clearer we have used calZ / 2 instead of Z in the last expression, because calZ vanishes for odd l, but Z doesn't.

Examples

Example1:

expression=X[LI[4],LI[-1],-a,-b,-c]X[LI[3],LI[2],-d,-e]

comp[expr_]:=ContractBasis[Basis[a,{2,-spherical}]Basis[b,{2,-spherical}]Basis[c,{2,-spherical}]Basis[d,{2,-spherical}]Basis[e,{2,-spherical}]expr]

r1=expression//comp//HarmonicComponent

expression/.ProductRule

%//comp//ToValues

%/.coeffERule

r2=%//HarmonicComponent

Simplify[r1-r2]

Compare total timings:

Timing[HarmonicComponent[comp[expression]]]

Timing[HarmonicComponent[ToValues@comp[expression/.ProductRule]/.coeffERule]]

Example2. There is a contracted index

expression=X[LI[4],LI[-1],-c,-a,-b]X[LI[3],LI[2],b,-d]

comp[expr_]:=ContractBasis[Basis[a,{2,-spherical}]Basis[c,{2,-spherical}]Basis[d,{2,-spherical}]expr]

r1=expression//comp//HarmonicComponent

expression/.ProductRule

%//comp//ToValues

%/.coeffERule

r2=%//HarmonicComponent

Simplify[r1-r2]

Example 3. Equal numbers of indices:

expression=X[LI[2],LI[-1],-a,-b]X[LI[2],LI[2],-d,-e]

comp[expr_]:=ContractBasis[Basis[a,{2,-spherical}]Basis[b,{2,-spherical}]Basis[d,{2,-spherical}]Basis[e,{2,-spherical}]expr]

r1=comp[expression]//HarmonicComponent

expression/.ProductRule

%//comp//ToValues

%/.coeffERule

r2=%//HarmonicComponent

Simplify[r1-r2]

Example 4. Equal numbers of indices and contracted indices:

expression=X[LI[2],LI[-1],-a,-b]X[LI[2],LI[2],a,b]

r1=expression//HarmonicComponent

expression/.ProductRule

%/.coeffERule

r2=%//HarmonicComponent

Simplify[r1-r2]

Example 5, with scalar harmonics:

Z[LI[3],LI[2]]Z[LI[3],LI[2]]/.ProductRule

%/.coeffERule

%//HarmonicComponent

SphericalHarmonicY[3,2,θ[],ϕ[]]^2

Simplify[%-%%]

Example 6, with an epsilon term:

Z[LI[4],LI[3],a,b]X[LI[4],LI[3],c,d]

%/.ProductRule

%/.coeffERule//Simplification

%//.mmbarRule//Simplification

r1=gamma[-b,-d]%//ContractMetric//Simplification

Z[LI[4],LI[3],a,b]X[LI[4],LI[3],c,-b]

%/.ProductRule

%/.coeffERule//Simplification

r2=%//.mmbarRule//Simplification

Simplify[r1-r2]

Example 7:

expression=Z[LI[l1],LI[m1],a]Z[LI[l2],LI[m2],-a]

tmp=expression/.FormalProductRule[l]

r1=%/.{l1->l2,m1->m2}

expression/.{l1->l2,m1->m2}

r2=%/.FormalProductRule[l]

Simplify[r1-r2]

Particular case:

TraceBasisDummy[Z[LI[2],LI[1],{a,spherical}]Z[LI[3],LI[-2],{-a,-spherical}]]

HarmonicComponent[%]

Sum[tmp/.{l1->2,m1->1,l2->3,m2->-2},{l,1,5}]

%/.coeffERule//HarmonicComponent

That is correct. However, here r1 and r2 are different:

expression=X[LI[l1],LI[m1],a]Z[LI[l2],LI[m2],-a]

tmp=expression/.FormalProductRule[l]

r1=tmp/.{l1->l2,m1->m2}

r2=expression/.{l1->l2,m1->m2}

Particular case:

tmp/.{l1->l2,m1->m2}/.{l2->2,m2->1}

Sum[%,{l,0,4}]

%/.coeffERule

Table[coeffE[{1,3,-1},{1,3,3},l],{l,0,6}]/.coeffERule

This is due to the fact that coeffE[{s1,l1,m1},{s2,l2,m2},l] is zero for odd l when l1=l2 and either s1=s2 or m1=m2. This comes from the fact that CG[{l,m},{l,m},{lp,2m}] is zero for odd lp.

Example 8:

expression=Z[LI[l1],LI[m1]]Z[LI[l2],LI[m2]]

tmp=%/.FormalProductRule[l]

r1=tmp/.{l1->l2,m1->m2}

expression/.{l1->l2,m1->m2}

r2=%/.FormalProductRule[l]

Simplify[r1-r2]

Particular case:

Z[LI[2],LI[1]]Z[LI[3],LI[-2]]

HarmonicComponent[%]

Sum[tmp/.{l1->2,m1->1,l2->3,m2->-2},{l,1,5}]

%/.coeffERule//HarmonicComponent

That is correct. Another particular case:

Z[LI[2],LI[1]]Z[LI[2],LI[1]]

HarmonicComponent[%]

Sum[tmp/.{l1->2,m1->1,l2->2,m2->1},{l,0,4}]

%/.coeffERule//HarmonicComponent

That is correct too.

Dump

End[]

EndPackage[]

TimeUsed[]

MemoryInUse[]/1024/1000.

MaxMemoryUsed[]/1024/1000.

SetDirectory["~/projects/secondpt/mathematica/"]

DumpSave["xPert_0.mx"]

Exit[]

Get["xPert_0.mx"]

MemoryInUse[]/1024/1000.