GV/PT/DT/D4-D2-D0 invariants

Accompanying notebook for the articles “Quantum Geometry, Stability and Modularity” [arxiv:2301.08066]
by Sergei Alexandrov, Soheyla Feyzbakhsh, Albrecht Klemm, Boris Pioline and Thorsten Schimannek
and ``Quantum Geometry and Mock Modularity” [arxiv:2312.nnnnn] by Sergei Alexandrov, Soheyla Feyzbakhsh, Albrecht Klemm, Boris Pioline

ClearAll[GVTab,DTTab,PTTab,PTTabGv,GVTabExt,GVTabExt2,h1,h1fun,h1int];
SetDirectory[NotebookDirectory[]];

(* Download from http://www.th.physik.uni-bonn.de/Groups/Klemm/data.php, and save same directory as this notebook *)
<<GVPTDTlib.m

(* Download the datafile X*.m for your favorite hypergeometric CY3 from http://www.th.physik.uni-bonn.de/Groups/Klemm/data.php *)
<<X10.m

(* maximal genus in the table of GV invariants *)
gmaxGV=Length[GVTab]-1;
(* highest genus attainable with standard boundary conditions *)
gInt=gMaxInteg[kappa,CYreg];
(* maximal degree in the table of GV invariants *)
QmaxGV=Length[Transpose[GVTab]];
(* extract Fourier coefficients in generating series of D4D2D0 invariants *)
JTab=CoefficientList[h1int,q];
(* restore leading powers of q in vector valued generating series *)
h1=Table[h1int[[mu+1]] q^(-(kappa+c2 )/24+FractionalPart[mu^2/(2kappa )+mu  /2]),{mu,0,Length[h1int]-1}];

(* CYname, Euler number, H^3, c2(TY).H, effective regulator rho, maximal attainable genus, maximal available genus *)
{CY,chi,kappa,c2,CYreg,QmaxGV,gInt,gmaxGV}

(* check CY condition and degree kappa *)
{Plus@@CYwts-Plus@@CYmod,kappa-Times@@CYmod/Times@@CYwts}

(* GV invariants saturating the Castelnuovo constraint, for degrees multiple of kappa *)
Table[GVSaturating[kappa,c2,kappa k],{k,1,10}]

(* Table of GV[g,d] row=genus, column= degree; data from direct integration *)
GVTab//MatrixForm

(* Table of DT[d,n] row=degree, 1st column=nmin *)
DTTab//MatrixForm

(* Table of PT[d,n] row=degree, 1st column=nmin *)
PTTab//MatrixForm

(* Table of PT[d,n] row=degree, 1st column=nmin; computed from GVTab + unknown Gv[g,d] *)
PTTabGv//MatrixForm

(* Table of GV[g,d] row=genus, column= degree; data from GVTab + modular constraints *)
GVTabExt//MatrixForm

(*  Same as previous table, but different format; row= degree, 1st column= gMax[d]; data from direct integration + modularity *)
GVTabExt2//MatrixForm

(* check GVTabExt2 against known GV invariants *)
Table[Join[{m,gMax[kappa,m]},Table[If[NumericQ[GVTabExt2[[m,delta+2]]],GVTabExt2[[m,delta+2]]-Gv[gMax[kappa,m]-delta,m]/.repGV,GVTabExt2[[m,delta+2]]],{delta,0,5}]],{m,1,Length[GVTabExt2]}]/.repGV//MatrixForm

(* GVTabExt2 with 1st column=degree, and ? replaced by known GV invariants; if result depends on unknown PT, leave ? *)
Table[Join[{m,gMax[kappa,m]},Table[If[GVTabExt2[[m,delta+2]]==="?"||GVTabExt2[[m,delta+2]]==="",Gv[gMax[kappa,m]-delta,m]/.repGV,If[NumericQ[GVTabExt2[[m,delta+2]]],GVTabExt2[[m,delta+2]],"?"]],{delta,0,5}]],{m,1,Length[GVTabExt2]}]/.repGV//MatrixForm

(* compare polar part with naive ansatz (4.11) *)
Series[Take[h1-ZPolarNaive[kappa,c2,1,Table[0,{mu,0,kappa-1}]]/.repDT,Floor[kappa/2+1]],{q,0,0}]

(* predictions from naive ansatz (4.11) for non-polar terms *)
With[{dn=2},Take[ZPolarNaiveInt[kappa,c2,1,Table[dn,{mu,0,kappa-1}]]/.repDT,Floor[kappa/2+1]]]

(* Extra constraints for the direct integration method; increase dn until unknown coefficients appear *)
With[{dn=5},GVPrediction[kappa,CYreg,gInt+25]]

Determination of PT and DT invariants from GV invariants

Qmax=Min[QmaxGV,Floor[-kappa/2+1/2Sqrt[kappa(8(gmaxGV-1)+kappa)]]]

NMax=100;
(* generating function of PT invariants *)
ZPTnew[QMax_,q_,y_]:=Exp[Sum[Sum[(((-1))^(g+1))GVTab[[g+1,Q]]Sum[(1/k)(q^((1-g)k)) (y^(Q k))(((1-(q^k)))^(2g-2)),{k,1,QMax}],{g,0,gMax[kappa,Q]}],{Q,1,QMax}]]/.q->-q;
(* y-expand the generating function of PT invariants *)
ZPTexpy=Simplify[Series[ZPTnew[Qmax,q,y],{y,0,Qmax}]];

(* q-expand each coefficient in y, after multiplying by the power of q corresponding to the Castelnuovo bound *)
ZPTexp=Table[Series[(q^(-MMin[kappa,i-1]))ZPTexpy[[3,i]],{q,0,NMax}],{i,1,Length[ZPTexpy[[3]]]}];
(* extract the coefficients of the expansion *)
PTTabsh=CoefficientList[ZPTexp,q];
(* add zeros to the first row to complete to a square table *)
PTTabsh[[1]]=Join[PTTabsh[[1]],Array[0 &,NMax]];
(* adds the first column encoding the shifts of the D0-charge *)
PTTab=Transpose[Join[{Table[-MMin[kappa,mu],{mu,0,Qmax}]},Transpose[PTTabsh]]];
MatrixForm[PTTab]

(* same for DT invariants *)
ZDTexp=Table[Series[MacMahon[-q,NMax]^chi (q^(-MMin[kappa,i-1]))ZPTexpy[[3,i]],{q,0,NMax}],{i,1,Length[ZPTexpy[[3]]]}];
(* extract the coefficients of the expansion *)
DTTabsh=CoefficientList[ZDTexp,q];
(* adds the first column encoding the shifts of the D0-charge *)
DTTab=Transpose[Join[{Table[-MMin[kappa,mu],{mu,0,Qmax}]},Transpose[DTTabsh]]];
MatrixForm[DTTab]

(* extending the table of GV invariants by unknown entries *)
Qmax=(Min[QmaxGV,Floor[-kappa/2+1/2Sqrt[kappa(8(gmaxGV-1)+kappa)]]])+4;
GVtemp=Array[Gv[#1-1,#2] &, {Max[Length[GVTab],gMax[kappa,Qmax]]+10,Length[Transpose[GVTab]]+10}];
GVtemp[[1;;Length[GVTab],1;;Length[Transpose[GVTab]]]]=GVTab;
NMax=60;
(* generating function of PT invariants *)
ZPTnew[QMax_,q_,y_]:=Exp[Sum[Sum[(((-1))^(g+1))GVtemp[[g+1,Q]]Sum[(1/k)(q^((1-g)k)) (y^(Q k))(((1-(q^k)))^(2g-2)),{k,1,QMax}],{g,0,gMax[kappa,Q]}],{Q,1,QMax}]]/.q->-q;

(* y-expand the generating function of PT invariants *)
ZPTexpy=Simplify[Series[ZPTnew[Qmax,q,y],{y,0,Qmax}]];

(* q-expand each coefficient in y, after multiplying by the power of q corresponding to the Castelnuovo bound *)
ZPTexp=Table[Series[(q^(-MMin[kappa,i-1]))ZPTexpy[[3,i]],{q,0,NMax}],{i,1,Length[ZPTexpy[[3]]]}];
(* extract the coefficients of the expansion *)
PTTabsh=Simplify[CoefficientList[ZPTexp,q]];
(* add zeros to the first row to complete to a square table *)
PTTabsh[[1]]=Join[PTTabsh[[1]],Array[0 &,NMax]];
(* adds the first column encoding the shifts of the D0-charge *)
PTTabGv=Transpose[Join[{Table[-MMin[kappa,mu],{mu,0,Qmax}]},Transpose[PTTabsh]]];
MatrixForm[PTTabGv]

Determination of J invariants from modular generating functions at rank 1

(* the generating function is h1fun in X.m *)
h1fun

NMax=100;

With[{r=1},
(* number of polar coefficients in elliptic genus: *)npolar=Plus@@Table[Ceiling[(kappa r^3+c2 r)/24-FractionalPart[mu^2/(2kappa r)+mu r /2]],{mu,0,r kappa/2}];
(* total number of polar coefficients in Ansatz: *)
npolart=Plus@@Ceiling[Table[-kappa r (1-r/2+k/(kappa r))^2/2+(4kappa r^3+c2 r)/24,{k,0,r kappa/2}]];
lmax=Catch[Do[
(* count monomials in e4 and e6 with total weight 2kappa r^3+r c2/2-2-Mod[kappa r,2]-2l *)
ListPowers=Table[Map[{Count[#,4],Count[#,6]}&,IntegerPartitions[2kappa r^3+r c2/2-2-Mod[kappa r,2]-2l,10000 (* some arbitrary large number !*),{4,6}]],{l,0,lmax}];
Vars=Flatten[Table[Table[s[ListPowers[[l+1,i,1]],ListPowers[[l+1,i,2]]],{i,Length[ListPowers[[l+1]]]}],{l,0,lmax}]];
If[Length[Vars]>=npolart,Throw[lmax]] ;,{lmax,1,10}]];]

DThTab={Table[Series[If[EvenQ[kappa],Th[1,kappa,k,q,NMax],dTh[1,kappa,k,q,NMax]],{q,0,NMax}],{k,0, kappa-1}]};
wTheta=If[EvenQ[kappa],1/2,3/2];
Do[
AppendTo[DThTab,Series[q D[Last[DThTab],q]-wTheta/12 E2[q]Last[DThTab],{q,0,NMax}]];
wTheta+=2,{l,lmax}];

h1ser=Assuming[q>0,Series[ReplaceAll[h1fun,{Eta->eta[q,NMax],e4->E4[q,NMax],e6->E6[q,NMax],DdTh[p_]:>DThTab[[p+1]],DTh[p_]:>DThTab[[p+1]]}],{q,0,NMax}]];

(* removing fractional part of the powers of q *)
h1int=If[kappa==1,{q^((kappa +c2 )/24) h1ser},Array[q^((kappa +c2 )/24-FractionalPart[(#-1)^2/(2kappa )+(#-1) /2]) h1ser[[#]] &,kappa]]

(* check eta-product formula for mu=0, ell=0 *)
Assuming[q>0,Series[DThTab[[1,1]]/If[Mod[kappa,2]==0,2(-1)^(kappa/2)eta[q^(2kappa),NMax]^2/eta[q^kappa,NMax],(-1)^((kappa-1)/2)kappa eta[q^kappa,NMax]^3],{q,0,NMax}]]

Compute D4D2D0 invariants with r=1 from Soheyla’s thm 1.1

(* Minimal values of the spectral flow parameter n *)
NTab=Array[nMin [kappa,#1-1,#2-1+MMin[kappa,#1-1]]&,{IntegerPart[kappa/2]+1,10}];
(* adds the first column encoding the shifts of the D0-charge *)
 NTab=Transpose[Join[{Table[-MMin[kappa,mu],{mu,0,IntegerPart[kappa/2]}]},Transpose[NTab]]];
MatrixForm[NTab]

(* Abelian D4-D2-D0 indices for minimal legal value of spectral flow *)
Print["Rigorous prediction from rank 0 wall-crossing:"];
DTs=Table[Table[S2[kappa,c2,mu,M,nMin[kappa,mu,M]],{M,MMin[kappa,mu],12+MMin[kappa,mu]}],{mu,0,kappa/2}]/.repPT;
DTs//MatrixForm
DTs/.repJ//MatrixForm

(* try if n-1 units of spectral flow already produces the correct result *)
DTs1=Table[Table[S2Lemma[kappa,c2,mu,M,nMin[kappa,mu,M]-1],{M,MMin[kappa,mu],12+MMin[kappa,mu]}],{mu,0,kappa/2}]/.repPT;
DTs1//MatrixForm
Simplify[DTs1/.repJ]//MatrixForm

(* table of (x,α) values for n-1 units of spectral flow *)
Charges1=Table[Table[mumToxa[kappa,{mu,M},nMin[kappa,mu,M]-1],{M,MMin[kappa,mu],12+MMin[kappa,mu]}],{mu,0,kappa/2}];
Charges1//MatrixForm

(* check if coefficients stabilize with increasing value of spectral flow*)
Table[Table[Table[S1[kappa,c2,mu,M,n],{n,1,nMin[kappa,mu,M]+1}],{M,MMin[kappa,mu],3+MMin[kappa,mu]}],{mu,0,kappa/2}]/.repPTGv

Compare modular prediction with Eq 4.27 in paper I, using n-dn units of spectral flow, where n is the minimal legal value; construct scatter plot

Rpointx={};Wpointx={};

(* result for n-dn *)
dn=1;mmax=15;
DTs=Table[Table[S2Lemma[kappa,c2,mu,M,nMin[kappa,mu,M]-dn],{M,MMin[kappa,mu],mmax+MMin[kappa,mu]}],{mu,0,kappa/2}]/.repPT;
JDT=DTs/.repJ//MatrixForm
JTab[[1;;IntegerPart[kappa/2]+1,1;;mmax]]//MatrixForm

Rpoint={};
Wpoint={};
For[mu=0,mu<=kappa/2,For[M=MMin[kappa,mu],M<=mmax+MMin[kappa,mu],(*Print[{mu,M}];*)
If[NumericQ[JDT[[1,mu+1,M-MMin[kappa,mu]+1]]],If[JDT[[1,mu+1,M-MMin[kappa,mu]+1]]==JTab[[mu+1,M-MMin[kappa,mu]+1]],
AppendTo[Rpoint,mumToxa[kappa,{mun[kappa,mu,nMin[kappa,mu,M]-dn],mn[kappa,mu,M,nMin[kappa,mu,M]-dn]}]],AppendTo[Wpoint,mumToxa[kappa,{mun[kappa,mu,nMin[kappa,mu,M]-dn],mn[kappa,mu,M,nMin[kappa,mu,M]-dn]}]]
]];M++];mu++];

Rpointx=Union[Flatten[Append[{Rpointx},Rpoint],1]];
Wpointx=Union[Flatten[Append[{Wpointx},Wpoint],1]];

{Length[Rpointx],Length[Wpointx]}

(* plot in (x,α) plane; BMT is the curve where the BMT final wall is tangent to parabola  *)
(* blue dots are points where Eq 4.27 is consistent with modular prediction, red otherwise *)
Pl=Plot[{fcon[x],x/2+1,3/4(x+1),Sqrt[2x]},{x,0,5},Epilog->{PointSize[0.02],Blue,Point[Rpointx],Red,Point[Wpointx]},PlotRange->{1,4},AxesLabel->{Style[x,14],Style[α,14]},Ticks->{{1,2,3,4},{1,2,3,4}},PlotLegends->{"f(x)","x/2+1","3/4(x+1)","BMT"}];
Show[Pl,Table[ParametricPlot[xan[x,0,n],{n,0,4},PlotStyle->Dotted],{x,0,4,1/2}]]

(* plot in (x,δ) plane; BMT is the curve where the BMT final wall is tangent to parabola  *)
(* blue dots are points where Eq 4.27 is consistent with modular prediction, red otherwise *)
Pl=Plot[{1/2x(x+1)-2/3x fcon[x],1/6x(x-1),0,-2/3 x(Sqrt[2x]-3/4(x+1))},{x,0,4.5},Epilog->{PointSize[0.02],Blue,Point[xaToxdelta/@Rpointx],Red,Point[xaToxdelta/@Wpointx]},PlotRange->{0,4},AxesLabel->{Style[x,14],Style[δ,14]},Ticks->{{1,2,3,4},{1,2,3,4}},PlotLegends->{"g(x)","1/6x(x-1)","Cast","BMT"}]

GV and PT invariants from Abelian D4-D2-D0 indices

(* approximately the minimal genus that can be treated by our method *)
gMin[kappa_,Q_]:=Floor[((Q^2)/(3 kappa))+((2Q)/3)+1];
(* the maximal degree which gives rise to the genus satisfying the Castelnuovo bound *)
QmaxGVg=Qmin[kappa,gmaxGV]

QMax=QmaxGVg+7;
dM=2;
UpBoundM[kappa_,Q_,Qm_]:=Max[1-gMin[kappa,Q],((((Qm^2)-6 Q Qm+3(Q^2)))/(6 kappa))-(Qm/6)-(Q/2)];

ZGV[QMax_,x_,y_]:=Sum[((MoebiusMu[k])/k)Log[1+Sum[Sum[(PT[Q,M]/.repPTfromJ/.repJ/.repPT) (((-1))^M)(((-x))^(k M)) (y^(k Q)),{M,MMin[kappa,Q],Ceiling[UpBoundM[kappa,Q,QMax]]+2dM}],{Q,1,QMax}]],{k,1,QMax}]
(* y-expand the generating function of PT invariants *)
ZGVexpy=Simplify[Series[ ZGV[QMax,x,y],{y,0,QMax}]];

GVTabExt=Array[0 &, {gMax[kappa,QMax]+1,QMax}];
GVTabExt2=Array["" &, {QMax,gMax[kappa,QMax]+2}];
Do[ZGVexp=ZGVexpy[[3,Q]];
GVTabExt2[[Q,1]]=gMax[kappa,Q];
Do[If[g>=Max[gMin[kappa,Q]-dM,0],
GVTabExt[[g+1,Q]]=Limit[(x^(g-1)) ZGVexp,x->0];
ZGVexp=Series[ZGVexp-GVTabExt[[g+1,Q]](x^(1-g))(((1+x))^(2(g-1))),{x,0,gMax[kappa,QMax]}],
GVTabExt[[g+1,Q]]=Gv[g,Q]];GVTabExt2[[Q,gMax[kappa,Q]-g+2]]=GVTabExt[[g+1,Q]],{g,gMax[kappa,Q],0,-1}],{Q,1,QMax}];
MatrixForm[GVTabExt]
MatrixForm[GVTabExt2]

(* combining known and predicted GV invariants *)
GVTabExtFull=Array[Gv[#1-1,#2] &, {Length[GVTabExt],Length[Transpose[GVTab]]}];
GVTabExtFull[[1;;Length[GVTab]]]=GVTab;
GVTabExtFull[[Length[GVTab]+1;;Length[GVTabExt],1;;Length[Transpose[GVTabExt]]]]=GVTabExt[[Length[GVTab]+1;;Length[GVTabExt]]];

(* using the combined table of GV invariants *)
Qmax=(Min[QmaxGV,Floor[-kappa/2+1/2Sqrt[kappa(8(gmaxGV-1)+kappa)]]])+4;
NMax=50;
(* generating function of PT invariants *)
ZPTnew[QMax_,q_,y_]:=Exp[Sum[Sum[(((-1))^(g+1))GVTabExtFull[[g+1,Q]]Sum[(1/k)(q^((1-g)k)) (y^(Q k))(((1-(q^k)))^(2g-2)),{k,1,QMax}],{g,0,gMax[kappa,Q]}],{Q,1,QMax}]]/.q->-q;
(* y-expand the generating function of PT invariants *)
ZPTexpy=Simplify[Series[ZPTnew[Qmax,q,y],{y,0,Qmax}]];

(* q-expand each coefficient in y, after multiplying by the power of q corresponding to the Castelnuovo bound *)
ZPTexp=Table[Series[(q^(-MMin[kappa,i-1]))ZPTexpy[[3,i]],{q,0,NMax}],{i,1,Length[ZPTexpy[[3]]]}];
(* extract the coefficients of the expansion *)
PTTabsh=Simplify[CoefficientList[ZPTexp,q]];
(* add zeros to the first row to complete to a square table *)
PTTabsh[[1]]=Join[PTTabsh[[1]],Array[0 &,NMax]];
(* adds the first column encoding the shifts of the D0-charge *)
PTTabExt=Transpose[Join[{Table[-MMin[kappa,mu],{mu,0,Qmax}]},Transpose[PTTabsh]]];
MatrixForm[PTTabExt]

Charge r=2 D4-D2-D0 invariants from rank 0 wall-crossing (unpublished)

(* maximal genus in the table of GV invariants *)
gmaxGV=Length[GVTab]-1;
(* highest genus attainable with standard boundary conditions *)
gInt=gMaxInteg[kappa,CYreg];
(* maximal degree in the table of GV invariants *)
QmaxGV=Length[Transpose[GVTab]];
(* extract Fourier coefficients in generating series of D4D2D0 invariants *)
JTab=CoefficientList[h1int,q];
(* restore leading powers of q in vector valued generating series *)
h1=Table[h1int[[mu+1]] q^(-(kappa+c2 )/24+FractionalPart[mu^2/(2kappa )+mu  /2]),{mu,0,Length[h1int]-1}];

Jr2=Table[Table[If[CondMP[kappa,mu,M],SPTDT[kappa,c2,mu,M,Ceiling[bb2[kappa,mu,M]]],J2[mu,M]],{M,M2Min[kappa,mu],4+M2Min[kappa,mu]}],{mu,0,kappa}]/.repPT/.repDT;
Print["Rigorous prediction from rank 0 wall-crossing:"];
Jr2//MatrixForm
Print["Overall powers of q:"];
Table[q^(-((kappa/3)+(c2/12)))q^(FractionalPart[((mu^2)/(4 kappa))+mu]),{mu,0,kappa}]
Print["Discriminant of BMT inequality:"];
Table[Table[BMTDiscriminant[ChernFromD4D2D0[kappa,{2,mu,M}]],{M,M2Min[kappa,mu],4+M2Min[kappa,mu]}],{mu,0,kappa}]//MatrixForm
Print["Naive ansatz from Eq. (5.20)in [arXiv:2204.02207] :"];
Jr2an=Table[Table[(((-1))^(M+ChiO[kappa,c2,2]+1))(ChiO[kappa,c2,2]-2 mu-M)DT[mu,M],{M,M2Min[kappa,mu],4+M2Min[kappa,mu]}],{mu,0,kappa}]/.repDT;
Jr2an//MatrixForm

Charge r=2 D4-D2-D0 invariants from rank -1 wall-crossing (Thm 1 in paper II)

(* Minimal values of the spectral flow parameter n *)
Nr2Tab=Array[nr2Min [kappa,#1-1,#2-1+Mr2Min[kappa,#1-1]]&,{kappa+1,5}];
(* adds the first column encoding the shifts of the D0-charge *)
 Nr2Tab=Transpose[Join[{Table[-Mr2Min[kappa,mu],{mu,0,kappa}]},Transpose[Nr2Tab]]];
MatrixForm[Nr2Tab]

(* D4D2D0 indices for minimal legal value of n *)
DTr2s=Table[Table[Sr2[kappa,c2,mu,M,nr2Min[kappa,mu,M]],{M,Mr2Min[kappa,mu],5+Mr2Min[kappa,mu]}],{mu,0,kappa}]/.repPT;
Simplify[DTr2s]//MatrixForm
DTr2=Simplify[DTr2s/.repJ/.repJ2];
DTr2//MatrixForm
Simplify[DTr2/.repPTfromJ/.repPT/.repJ]//MatrixForm
DTR2=Simplify[DTr2/.repPTExt];
DTR2//MatrixForm

Which[(* specific for X10 *)
CY == "X10",
repJ2S={J2[0,1]->DTR2[[1,2]],J2[0,2]->DTR2[[1,3]],J2[1,-1]->DTR2[[2,1]],J2[1,0]->DTR2[[2,2]],J2[1,1]->DTR2[[2,3]]},
(* specific for X8 *)
CY == "X8",
repJ2S={J2[0,1]->DTR2[[1,2]],J2[0,2]->DTR2[[1,3]],J2[1,-1]->DTR2[[2,1]],J2[1,0]->DTR2[[2,2]],J2[1,1]->DTR2[[2,3]],J2[1,2]->DTR2[[2,4]],J2[2,-2]->DTR2[[3,1]]},
(* specific for X66 *)
CY == "X66",
repJ2S={J2[1,-1]->-1/2,a[22,7,0]->(638909706639409329426483218315175219485418151277413253009/23098972503443511196547086232111387166720000000)},
(* other CYs *)
True,repJ2S={}];

(* try if n-1 units of spectral flow already produces the correct result *)
DTr2s1m=Table[Table[Sr2[kappa,c2,mu,M,nr2Min[kappa,mu,M]-1],{M,Mr2Min[kappa,mu],6+Mr2Min[kappa,mu]}],{mu,0,kappa}]/.repPT;
Simplify[DTr2s1m]//MatrixForm
DTr21m=Simplify[DTr2s1m/.repJ/.repJ2/.repJ2S];
DTr21m//MatrixForm
DTR21m=Simplify[DTr21m/.repPTExt/.repJ2S];
DTR21m//MatrixForm

(* table of (x,\alpha) values for n-1 units of spectral flow *)
Charges2=Table[Table[mumToxar2[kappa,{mu,M},nr2Min[kappa,mu,M]-1],{M,Mr2Min[kappa,mu],8+Mr2Min[kappa,mu]}],{mu,0,kappa}];
Charges2//MatrixForm

(* try if n-2 units of spectral flow already produces the correct result *)
DTr2s2m=Table[Table[Sr2[kappa,c2,mu,M,nr2Min[kappa,mu,M]-2],{M,Mr2Min[kappa,mu],6+Mr2Min[kappa,mu]}],{mu,0,kappa}]/.repPT;
Simplify[DTr2s2m]//MatrixForm
DTr22m=Simplify[DTr2s2m/.repJ/.repJ2/.repJ2S];
DTr22m//MatrixForm
DTR22m=Simplify[DTr22m/.repPTExt/.repJ2S];
DTR22m//MatrixForm

(* table of (x,\alpha) values for n-2 units of spectral flow *)
Charges2=Table[Table[mumToxar2[kappa,{mu,M},nr2Min[kappa,mu,M]-2],{M,Mr2Min[kappa,mu],8+Mr2Min[kappa,mu]}],{mu,0,kappa}];
Charges2//MatrixForm

(* Try with n+1 units of spectral flow *)
DTr2s1=Table[Table[Sr2[kappa,c2,mu,M,nr2Min[kappa,mu,M]+1],{M,Mr2Min[kappa,mu],5+Mr2Min[kappa,mu]}],{mu,0,kappa}]/.repPT;
Simplify[DTr2s1]//MatrixForm
DTr21=Simplify[DTr2s1/.repJ/.repJ2/.repJ2S];
DTr21//MatrixForm
DTR21=Simplify[DTr21/.repPTExt/.repJ2S];
DTR21//MatrixForm

(* try with n+2 units of spectral flow *)
DTr2s2=Table[Table[Sr2[kappa,c2,mu,M,nr2Min[kappa,mu,M]+2],{M,Mr2Min[kappa,mu],5+Mr2Min[kappa,mu]}],{mu,0,kappa}]/.repPT;
Simplify[DTr2s2]//MatrixForm
DTr22=Simplify[DTr2s2/.repJ/.repJ2/.repJ2S];
DTr22//MatrixForm
DTR22=Simplify[DTr22/.repPTExt];
DTR22//MatrixForm

(* constraints on GV invariants for X10 *)
PT1365=Solve[DTR22[[2,2]]==DTR2[[2,2]],PT[13,-65]][[1,1]]
DTR21/.PT1365//MatrixForm
DTR2/.PT1365//MatrixForm
DTR21m/.PT1365//MatrixForm

(* additional check for X8 *)
J2Tab=CoefficientList[h2int,q];
PT18=Solve[{DTR2[[3,2]]==J2Tab[[3,2]],DTR2[[3,3]]==J2Tab[[3,3]]+(1/4)JTab[[1,2]],DTR2[[3,4]]==J2Tab[[3,4]]},{PT[18,-66],PT[18,-65],PT[18,-64]}][[1]]
DTR21[[3,1]]/.PT18

(* finding unkown coefficient for X66 *)
sa=Solve[DTR21[[1,1]]==2,PT[8,-24]][[1,1]]
DTR2[[1,2]]/.sa

Chamber structure for rank -1 charges used to compute r=2 D4-D2-D0 invariants

(* table of most polar terms in generating series of charge 2 D4-D2-D0, using rank 0 wall-crossing *)
Jr2=Table[Table[If[CondMP[kappa,mu,M],SPTDT[kappa,c2,mu,M,Ceiling[bb2[kappa,mu,M]]],J2[mu,M]],{M,M2Min[kappa,mu],4+M2Min[kappa,mu]}],{mu,0,kappa}]/.repPT/.repDT

Rankm1Chamber[ Q_, m_] :=Module[{Li,Ind,Gr},Li={};
	Gr={Dotted,BMT[{-kappa,0,Q,-m}],Dashing[None],Flatten[{
(* first P1 *)
Table[
Ind=(((-1))^(Chi1[kappa, c2,Q, m,Qp,mp]))Chi1[kappa, c2,Q, m,Qp,mp]PTw[kappa,c2,Qp,mp,((2Q-3Qp)/kappa)-1] J[SFmu[kappa,Q-Qp],SFM[kappa,Q-Qp,m-mp+Qp]]/.repPT/.repJ;
If[Ind===0,{},AppendTo[Li,{"(Subscript[P, 1]):",Ind,Qp,mp}];
{Hue[.8-.4(Qp)/Floor[UL1Qp[kappa,Q,m]]],WMS[{-kappa,0,Q,-m},{0,kappa,1/2kappa+Q-Qp,1/6kappa+mp-m-Qp},"(Subscript[P, 1]):"<>ToString[Ind]]}],{Qp,0,Floor[UL1Qp[kappa,Q,m]]},
{mp,Ceiling[LL1mp[kappa,Q,Qp]],Floor[UL1mp[kappa,Q,m,Qp]]}],
(* now P2 *)
th+=2dth;Table[
Ind=(1/2)(((-1))^(Chi2[kappa,c2, Q, mup,Mp]+Chi2[kappa, c2,Q,mup,Mpp]+1))
Chi2[kappa, c2,Q, mup,Mp]Chi2[kappa,c2, Q, mup,Mpp]
PT[Q-2mup+3kappa,m-Mp-Mpp+2Q-6mup+5kappa] J[SFmu[kappa,mup],SFM[kappa,mup,Mp]]J[SFmu[kappa,mup],SFM[kappa,mup,Mpp]]/.repPT/.repJ;
If[Ind===0,{},AppendTo[Li,{"(Subscript[P, 2]):",Ind,mup,Mp,Mpp}];
{Hue[.9-(mup)/10/Floor[UL1Qp[kappa,Q,m]]],
(* to be checked :*)
WMS[{-kappa,0,Q,-m},{0,kappa,1/2kappa+mup,1/6kappa-Mp},"(Subscript[P, 2]):"<>ToString[Ind]]}],
{mup,Ceiling[(kappa/2)(1-((3m)/Q))],Floor[(1/2)(Q+3kappa)]},{Mp,Ceiling[LL2Mp[kappa,mup]],Floor[UL2sumMp[kappa,Q,m,mup]]-Ceiling[LL2Mp[kappa,mup]]},{Mpp,Ceiling[LL2Mp[kappa,mup]],Floor[UL2sumMp[kappa,Q,m,mup]]-Mp}],
(* now P3 *)
th=Pi/3;Table[
Ind=(((-1))^(Chi3[kappa, c2,Q, m,Qp,mp]))Chi3[kappa,c2, Q, m,Qp,mp]PT[Qp,mp] J2[SFmur2[kappa,Q-Qp],SFMr2[kappa,Q-Qp,m-mp+2Qp]]/.repPT/.repJ/.repJ2/.repJ22;
If[Ind===0,{},AppendTo[Li,{"(Subscript[P, 3]):",Ind,Qp,mp}];
{Hue[1/20-(Qp)/20/Floor[UL1Qp[kappa,Q,m]]],
(* need to fix this:*)
WMS[{-kappa,0,Q,-m},{0,2kappa,2kappa+Q-Qp,4/3kappa+mp-m-2Qp},"(Subscript[P, 3]):"<>ToString[Ind]]}]
,{Qp,0,Floor[UL3Qp[kappa,Q,m]]},{mp,Ceiling[LL3mp[kappa,Qp]],Floor[UL3mp[kappa,Q,m,Qp]]}]
},2]};
Print[Li,"->",{PT[Q,m],Sum[Li[[i,2]],{i,Length[Li]}]}/.repPT];Gr];

(* customize this routine to beautify the plots below: *)
 HalfCircle[b0_,b1_,t_]:=Module[{},th+=dth;If[(th>.8Pi/2&&dth>0),th=Pi-.1;dth=-1.3dth];
(*th=N[ArcTan[1/20Log[Abs[ToExpression[StringSplit[t,":"][[2]]]]]]]; *){Circle[{Simplify[(b0+b1)/2],0},Abs[Simplify[b1-b0]]+RandomReal[{0,0.05}],{0,Pi}],Text[Style[t,12],{(b0+b1)/2,0}+Abs[Simplify[b1-b0]]{Cos[th],Sin[th]}]}];

repJ22={}; (* additional rule for evaluating J2[]*)

Chamber structure for leading terms in X10:

(* X10 *) repJ22={J2[1,-1]->-21/4};

(* initial angle for the labels: *) th=-Pi/30+.05;
(* increment angle between consecutive walls: *) dth=Pi/30;
With[{mu=0,M=0,d=0},
Print["Minimal value of M:",Mr2Min[kappa,mu]];
Print["Minimal legal value of n:",nr2Min[kappa,mu,M]];
Mun=munr2[kappa,mu,nr2Min[kappa,mu,M]-d];
Mn=mnr2[kappa,mu,M,nr2Min[kappa,mu,M]-d];
Print["Chamber structure for γ=["<>ToString[-kappa]<>","<>ToString[Mun]<>",0,"<>ToString[-Mn]<>"]:"];
Show[Graphics[Rankm1Chamber[Mun,Mn]],Axes->{True,False}]]

PT[6,-15]/.repPT

(* initial angle for the labels: *) th=-Pi/30+.05;
(* increment angle between consecutive walls: *) dth=Pi/30;
With[{mu=1,M=-1,d=0},
Print["Minimal value of M:",Mr2Min[kappa,mu]];
Print["Minimal legal value of n:",nr2Min[kappa,mu,M]];
Mun=munr2[kappa,mu,nr2Min[kappa,mu,M]-d];
Mn=mnr2[kappa,mu,M,nr2Min[kappa,mu,M]-d];
Print["Chamber structure for γ=["<>ToString[-kappa]<>","<>ToString[Mun]<>",0,"<>ToString[-Mn]<>"]:"];
Show[Graphics[Rankm1Chamber[Mun,Mn]],Axes->{True,False}]]

PT[7,-19]/.repPT

Chamber structure for leading terms in X8:

(*X8 *) repJ22={J2[2,-2]->-17};

(* initial angle for the labels: *) th=-Pi/30+.05;
(* increment angle between consecutive walls: *) dth=Pi/30;
With[{mu=0,M=0,d=0},
Print["Minimal value of M:",Mr2Min[kappa,mu]];
Print["Minimal legal value of n:",nr2Min[kappa,mu,M]];
Mun=munr2[kappa,mu,nr2Min[kappa,mu,M]-d];
Mn=mnr2[kappa,mu,M,nr2Min[kappa,mu,M]-d];
Print["Chamber structure for γ=["<>ToString[-kappa]<>","<>ToString[Mun]<>",0,"<>ToString[-Mn]<>"]:"];
Show[Graphics[Rankm1Chamber[Mun,Mn]],Axes->{True,False}]]

(* initial angle for the labels: *) th=-Pi/30+.05;
(* increment angle between consecutive walls: *) dth=Pi/30;
With[{mu=1,M=-1,d=0},
Print["Minimal value of M:",Mr2Min[kappa,mu]];
Print["Minimal legal value of n:",nr2Min[kappa,mu,M]];
Mun=munr2[kappa,mu,nr2Min[kappa,mu,M]-d];
Mn=mnr2[kappa,mu,M,nr2Min[kappa,mu,M]-d];
Print["Chamber structure for γ=["<>ToString[-kappa]<>","<>ToString[Mun]<>",0,"<>ToString[-Mn]<>"]:"];
Show[Graphics[Rankm1Chamber[Mun,Mn]],Axes->{True,False}]]

(* initial angle for the labels: *) th=-Pi/30+.05;
(* increment angle between consecutive walls: *) dth=Pi/30;
With[{mu=2,M=-2,d=0},
Print["Minimal value of M:",Mr2Min[kappa,mu]];
Print["Minimal legal value of n:",nr2Min[kappa,mu,M]];
Mun=munr2[kappa,mu,nr2Min[kappa,mu,M]-d];
Mn=mnr2[kappa,mu,M,nr2Min[kappa,mu,M]-d];
Print["Chamber structure for γ=["<>ToString[-kappa]<>","<>ToString[Mun]<>",0,"<>ToString[-Mn]<>"]:"];
Show[Graphics[Rankm1Chamber[Mun,Mn]],Axes->{True,False}]]

SJ[kappa,c2,18,-66]/.repPT/.repJ

SJ[kappa,c2,18,-65]/.repPT/.repJ

SJ[kappa,c2,18,-64]/.repPT/.repJ

mumToxa[kappa,{18,-64}]

r=2 D4-D2-D0 invariants from mock modular generating functions

<<X10.m

(* restore leading powers of q in vector valued generating series *)
h1=Table[h1int[[mu+1]] q^(-(kappa+c2 )/24+FractionalPart[mu^2/(2kappa )+mu  /2]),{mu,0,Length[h1int]-1}];

(* the generating function of rational DT invariants is h2fun+(-1)^(mu+kappa)G[mu+1]h1fun^2 in CY.m,
where h2fun is a weakly holomorphic vector valued modular form and G[mu] a mock modular form with specific shadow*)
h2fun

NMax=25;

(* Generating series of first 100 Hurwitz class numbers, using data from https://oeis.org/A058306 *)
ClassNum={-1,0,0,1,1,0,0,1,1,0,0,1,4,0,0,2,3,0,0,1,2,0,0,3,2,0,0,4,2,0,0,3,3,0,0,2,5,0,0,4,2,0,0,1,4,0,0,5,10,0,0,2,2,0,0,4,4,0,0,3,4,0,0,5,7,0,0,1,4,0,0,7,3,0,0,7,4,0,0,5,6,0,0,3,4,0,0,6,2,0,0,2,6,0,0,8,6,0,0,3}/{12,1,1,3,2,1,1,1,1,1,1,1,3,1,1,1,2,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
ClassF[kappa_,mu_,q_]:=If[mu==0,-1/12,0]+Sum[ClassNum[[4n+1-mu]] q^(n-mu/4),{n,1,Min[(Length[ClassNum]-1)/4,kappa NMax]}];
(* the same as a function of tau, to avoid branch cuts when computing action of Hecke operator *)
ClassFt[kappa_,mu_,tau_]:=If[mu==0,-1/12,0]+Sum[ClassNum[[4n+1-mu]] Exp[2Pi I tau (n-mu/4)],{n,1,Min[(Length[ClassNum]-1)/4,kappa NMax]}];

(* compute image under of generating series of Hurwitz class number under Hecke_kappa ASSUMING kappa is prime *)
HeckeClassFprime=FullSimplify[Assuming[{q>0,q<1},Series[Table[kappa DivisorSum[kappa,If[IntegerQ[mu #/kappa],1/#^2 Sum[Exp[I Pi mu^2 b #/2/kappa^2]ClassFt[kappa,Mod[mu #/kappa,2],kappa/#^2 tau+b/#],{b,0,#-1}],0]&],{mu,0,2kappa-1}]/.tau->Log[q]/(2Pi I),{q,0,NMax}]]]

Which[ (* prescription in B.16 [arXiv:2204.02207] for solution to mock-modular anomaly *)
kappa==1,HeckeClassF=HeckeClassFprime,
PrimeQ[kappa],HeckeClassF=1/2HeckeClassFprime,
kappa==4,HeckeClassF=1/2HeckeClassFprime-1/8Table[If[Mod[mu,2]==0,HeckeClassFprime[[mu+1]]+HeckeClassFprime[[Mod[mu+4,2kappa]+1]],0],{mu,0,2kappa-1}],
kappa==8,
HeckeClassF=1/2HeckeClassFprime-1/6Table[If[Mod[mu,2]==0,HeckeClassFprime[[mu+1]]+HeckeClassFprime[[Mod[mu+8,2kappa]+1]],0],{mu,0,2kappa-1}],
kappa==9,
HeckeClassF=
1/2HeckeClassFprime-1/10Table[If[Mod[mu,3]==0,HeckeClassFprime[[mu+1]]+HeckeClassFprime[[Mod[mu+6,2kappa]+1]]+HeckeClassFprime[[Mod[mu+12,2kappa]+1]],0],{mu,0,2kappa-1}], 
kappa==16,1/2HeckeClassFprime-Table[1/6If[Mod[mu,2]==0,HeckeClassFprime[[mu+1]]+HeckeClassFprime[[Mod[mu+16,2kappa]+1]]
+HeckeClassFprime[[Mod[mu+12,2kappa]+1]],0]+1/60 If[Mod[mu,4]==0,HeckeClassFprime[[mu+1]]+HeckeClassFprime[[Mod[mu+8,2kappa]+1]]+HeckeClassFprime[[Mod[mu+16,2kappa]+1]]
+HeckeClassFprime[[Mod[mu+24,2kappa]+1]],0],{mu,0,2kappa-1}], 
True,Print["I don't know to evaluate the modified Hecke operator !"]];

(* determine number of Serre derivatives needed to saturate the expected dimension of space of weakly holomorphic modular forms *)
With[{r=2},npolar=Plus@@Table[Ceiling[(kappa r^3+c2 r)/24-FractionalPart[mu^2/(2kappa r)+mu r /2]],{mu,0,r kappa/2}];
(* total number of polar coefficients in Ansatz: *)
npolart=Plus@@Ceiling[Table[-kappa r (1-r/2+k/(kappa r))^2/2+(4kappa r^3+c2 r)/24,{k,0,r kappa/2}]];
lmax=Catch[Do[
(* count monomials in e4 and e6 with total weight 2kappa r^3+r c2/2-2-Mod[kappa r,2]-2l *)
ListPowers=Table[Map[{Count[#,4],Count[#,6]}&,IntegerPartitions[2kappa r^3+r c2/2-2-Mod[kappa r,2]-2l,10000 (* some arbitrary large number !*),{4,6}]],{l,0,lmax}];
Vars=Flatten[Table[Table[s[ListPowers[[l+1,i,1]],ListPowers[[l+1,i,2]]],{i,Length[ListPowers[[l+1]]]}],{l,0,lmax}]];
If[Length[Vars]>=npolart,Throw[lmax]] ;,{lmax,1,10}]];
]

(* construct table of iterated Serre derivatives of theta series up to order lmax *)
DThTab={Table[Series[Th[2,kappa,k,q,NMax],{q,0,NMax}],{k,0, kappa}]};
wTheta=1/2;
Do[
AppendTo[DThTab,Series[q D[Last[DThTab],q]-wTheta/12 E2[q]Last[DThTab],{q,0,NMax}]];
wTheta+=2,{l,lmax}];

(* explicit q-expansions of Eisenstein series, to speed up computation *)
E2[q_]:=1-24 q-72 (q^2)-96 (q^3)-168 (q^4)-144 (q^5)-288 (q^6)-192 (q^7)-360 (q^8)-312 (q^9)-432 (q^10)-288 (q^11)-672 (q^12)-336 (q^13)-576 (q^14)-576 (q^15)-744 (q^16)-432 (q^17)-936 (q^18)-480 (q^19)-1008 (q^20)-768 (q^21)-864 (q^22)-576 (q^23)-1440 (q^24)-744 (q^25)-1008 (q^26)-960 (q^27)-1344 (q^28)-720 (q^29)-1728 (q^30)-768 (q^31)-1512 (q^32)-1152 (q^33)-1296 (q^34)-1152 (q^35)-2184 (q^36)-912 (q^37)-1440 (q^38)-1344 (q^39)-2160 (q^40)-1008 (q^41)-2304 (q^42)-1056 (q^43)-2016 (q^44)-1872 (q^45)-1728 (q^46)-1152 (q^47)-2976 (q^48)-1368 (q^49)-2232 (q^50);
E4[q_]:=1+240 q+2160 (q^2)+6720 (q^3)+17520 (q^4)+30240 (q^5)+60480 (q^6)+82560 (q^7)+140400 (q^8)+181680 (q^9)+272160 (q^10)+319680 (q^11)+490560 (q^12)+527520 (q^13)+743040 (q^14)+846720 (q^15)+1123440 (q^16)+1179360 (q^17)+1635120 (q^18)+1646400 (q^19)+2207520 (q^20)+2311680 (q^21)+2877120 (q^22)+2920320 (q^23)+3931200 (q^24)+3780240 (q^25)+4747680 (q^26)+4905600 (q^27)+6026880 (q^28)+5853600 (q^29)+7620480 (q^30)+7150080 (q^31)+8987760 (q^32)+8951040 (q^33)+10614240 (q^34)+10402560 (q^35)+13262640 (q^36)+12156960 (q^37)+14817600 (q^38)+14770560 (q^39)+17690400 (q^40)+16541280 (q^41)+20805120 (q^42)+19081920 (q^43)+23336640 (q^44)+22891680 (q^45)+26282880 (q^46)+24917760 (q^47)+31456320 (q^48)+28318320 (q^49)+34022160 (q^50);
E6[q_]:=1-504 q-16632 (q^2)-122976 (q^3)-532728 (q^4)-1575504 (q^5)-4058208 (q^6)-8471232 (q^7)-17047800 (q^8)-29883672 (q^9)-51991632 (q^10)-81170208 (q^11)-129985632 (q^12)-187132176 (q^13)-279550656 (q^14)-384422976 (q^15)-545530104 (q^16)-715608432 (q^17)-986161176 (q^18)-1247954400 (q^19)-1665307728 (q^20)-2066980608 (q^21)-2678616864 (q^22)-3243917376 (q^23)-4159663200 (q^24)-4923450504 (q^25)-6175361808 (q^26)-7261732800 (q^27)-8954092224 (q^28)-10337619600 (q^29)-12685958208 (q^30)-14429092608 (q^31)-17456963832 (q^32)-19805530752 (q^33)-23615078256 (q^34)-26481071232 (q^35)-31587041304 (q^36)-34949354832 (q^37)-41182495200 (q^38)-45660250944 (q^39)-53291422800 (q^40)-58391525808 (q^41)-68210360064 (q^42)-74092255776 (q^43)-85796909856 (q^44)-93416358672 (q^45)-107049273408 (q^46)-115589884032 (q^47)-133109345376 (q^48)-142375996728 (q^49)-162473866632 (q^50);

(* expansion of weakly holomorphic h_{2,mu}^{(0)} *)
ZWeak2=Table[Series[h2fun/.{e4->E4[q],e6->E6[q],Eta->eta[q],DTh[l_]:>DThTab[[l+1,k]]},{q,0,NMax}],{k,Length[h2fun]}]

(* mock-modular contribution *) ZMock2=With[{r=2},Table[Series[-Sum[(-1)^(mu-2mu1+kappa)HeckeClassF[[Mod[mu-2mu1+kappa,2kappa]+1]]h1[[Mod[mu1,kappa]+1]]h1[[Mod[mu-mu1,kappa]+1]],{mu1,0,kappa-1}],{q,0,NMax}],{mu,0,r kappa/2}]]

(* contribution to charge 2 rational D4-D2-D0 from Abelian D4-D2-D0 invariants *) 
ZAbelian2=With[{r=2},Table[Series[If[Mod[mu+kappa,r]==0,1/r^2Simplify[h1[[Mod[(mu+kappa)/r,kappa]+1]]/.q->q^r,q>0],0],{q,0,NMax}],{mu,0,kappa}]]

Print["Modularity predicts the following rank 2 rational DT invariants:"];
ZRat2=Simplify[ZWeak2-ZMock2]

Print["Modularity predicts the following rank 2 integer DT invariants:"];
ZInt2=Simplify[ZWeak2-ZAbelian2-ZMock2]

(* remove leading powers of q in vector valued generating series *)
h2rat=With[{r=2},Table[ZRat2[[mu+1]] q^((kappa r^3+c2 r )/24-FractionalPart[mu^2/(2kappa r )+mu r /2]),{mu,0,Length[ZInt2]-1}]]

(* remove leading powers of q in vector valued generating series *)
h2int=With[{r=2},Table[ZInt2[[mu+1]] q^((kappa r^3+c2 r )/24-FractionalPart[mu^2/(2kappa r )+mu r /2]),{mu,0,Length[ZInt2]-1}]]