problem 2

This is exponent for S(ω)(e^(-i ω t)). We are supposed to expand it around (Subscript[ω, 0]), and the final result shows time advance.

I a/(ω0-ω-I G)- I ω t
Normal@Series[%,{ω,ω0,1}]
Collect[%,ω,FullSimplify]

This is a plot showing (Subscript[f, in])(t), where for simplicity we divide by oscillating part to stablize it

(E^(-((t^2)/(2 (σ^2)))-I t ω0))/Exp[-I t ω0]/.{σ->1,G->15,a->1125}/.ω0->1;
Plot[Re[%],{t,-6,6}]

Take some example parameters, we are below computing the output signals (Subscript[f, out])(t), where G is Γ in the problem set. To make it more precise, we set a pretty high precision and working precision. To be consistent with input, we also divide it by Exp[-i (Subscript[ω, 0]) t]

Exp[-I ω t]Exp[I a/(ω0-ω-I G)]Exp[-(ω-ω0)^2 σ^2/2]/Sqrt[2Pi]σ /.{σ->1,G->15,a->1125}
%/Exp[-I ω0 t]/.ω0->1
ddo[t_]:=Re[NIntegrate[((E^(I t+((1125 I)/((1-15 I)-ω))-(1/2) (((-1+ω))^2)-I t ω))/(Sqrt[2 π]))//SetPrecision[#,100]&,{ω,-Infinity,Infinity},WorkingPrecision->60]]

You can now plot output signal! The way I suggest is to do ListPlot. It may take some time, if you have quicker way to realize it, please tell me! You will notice that, for the example parameters, the Gaussian is not distorted too much with a time advance around 5σ (the peak now is around -5 rather than 0). However, the magnitude of the signal is highly suppressed!

tab=Monitor[Table[{t,ddo[t]},{t,-8,1,0.01}],t];//AbsoluteTiming

ListPlot[tab]

It doesn’t violate causality even though it has time advance. Because its norm is inconsequential compared to input signal, so you shall find it’s triggered by negative-time-tail of input signal.

problem 3

real scalar

Evaluation of low-energy part of sum rule

Mlow=-g[0]+g[2](s^2+t^2+u^2)+g[3] s t u+g[4](s^2+t^2+u^2)^2/.{u->-s-t};

evalMlow[k_(*this is the number of subtraction, even number*)][expr_]:=Module[{todo=expr/(s(s(s+t))^(k/2))},Residue[todo,{s,0}]+Residue[todo,{s,-t}]]

light=Table[evalMlow[i][Mlow]//Expand,{i,2,4,2}]

Evaluation at high-energy: you should use partial-wave expansion. Let’s first consider real scalar.

Pj[j_,x_,d_,nders_]:=(Hypergeometric2F1[-j+nders,j+d-3+nders,((d-2)/2)+nders,((1-x)/2)] (Pochhammer[-j,nders] Pochhammer[j+d-3,nders]))/((((-2))^nders) Pochhammer[((d-2)/2),nders])(*this is the Legendre function in general dimensions with any derivatives*)

(Table[(2s+t)/((s(s+t))^(k/2+1))pj[J,1+2t/s,4,0]/.pj->Pj,{k,2,4,2}]/.s->m^2)m^2(*this additonal m^2 is due to the definition of average*);
heavy=Normal@Series[%,{t,0,2}]

Let’s define the ave operator to make eq (1) of problem set more manifest. Although, this is not necessary, you can eyeball it easily!

aveofheavy=ave/@heavy//.{ave[a_+b_]:>ave[a]+ave[b],ave[a_ t^(pow_:1)]:>ave[a]t^pow}

Flatten@CoefficientList[aveofheavy-light,t]

(*to get null constraint, we can take two g[4]*)

Flatten@CoefficientList[aveofheavy-light,t];
{%[[3]],-2%[[4]]}//Expand
Total@%/.ave[a_]:>a//FullSimplify

The null constraint is actually from crossing symmetry, to make it manifest, we can consider full 2-subtracted sum rule instead of considering 2-subtraction and 4-subtraction.

evalMlowfull[k_(*this is the number of subtraction, even number*)][expr_]:=Module[{todo=expr/((s-sp)(s(s+t))^(k/2))},Residue[todo,{s,0}]+Residue[todo,{s,-t}]+Residue[todo,{s,sp}]]

lightfull=Collect[Normal@Series[evalMlowfull[2][Mlow],{sp,0,2},{t,0,2}],{sp,t}]

((2s+t)/((s-sp)(s+t+sp)(s(s+t))^(k/2))pj[J,1+2t/s,4,0]/.pj->Pj/.k->2/.s->m^2)m^2(*this additonal m^2 is due to the definition of average*);
heavyfull=Collect[Normal@Series[%,{sp,0,2},{t,0,2}],{sp,t},Simplify]

(*it's all equivalent to above*)

complex scalar

Mlowcomp[s_,t_]:=-g[0]+g[1]t+g[2,1]t^2+g[2,2]s u+g[3,1]s t u+g[3,2]t^3+g[4,1]t^4+g[4,2]t^2 s u+g[4,3]s^2 u^2/.{u->-s-t};

light1=Normal@Series[evalMlowfull[2][Mlowcomp[s,t]],{sp,0,2},{t,0,2}]
light2=Normal@Series[evalMlowfull[2][Mlowcomp[s,-s-t]+Mlowcomp[-s-t,s]],{sp,0,2},{t,0,2}]

((2s+t)/((s-sp)(s+t+sp)(s(s+t))^(k/2))pj[J,1+2t/s,4,0]/.pj->Pj/.k->2/.s->m^2)m^2(*this additonal m^2 is due to the definition of average*);
heavy1=Collect[Normal@Series[Normal@Series[% cpb ,{sp,0,4},{t,0,4}]/.{sp->sp ep,t->t ep},{ep,0,2}]/.ep->1,{sp,t},FullSimplify];
heavy2=Collect[Normal@Series[Normal@Series[%%(cpp+cpb),{sp,0,4},{t,0,4}]/.{sp->sp ep,t->t ep},{ep,0,2}]/.ep->1,{sp,t},FullSimplify];

evalMlowfullodd[k_(*this is the number of subtraction, even number*)][expr_]:=Module[{todo=expr/((s-sp)(s+t)^(k))},Residue[todo,{s,0}]+Residue[todo,{s,-t}]+Residue[todo,{s,sp}]]

(*if we include single subtraction, change 3 to 1*)

light3=Collect[Normal@Series[evalMlowfullodd[1][Mlowcomp[s,-s-t]-Mlowcomp[-s-t,s]],{sp,0,3},{t,0,3}],{sp,t}]

listwhat[expr_,what_]:=DeleteDuplicates[Cases[expr,what,{0,Infinity}]](*auxiliary code, to list what you want*)

Let’s then eyeball rep1, rep2 and rep3 below to count how many null constraints we have

((1/((s-sp)(s+t)^k)+1/((s+t+sp)s^k))pj[J,1+2t/s,4,0]/.pj->Pj/.{s->m^2})m^2/.k->1//Together;
heavy3=Collect[Normal@Series[Normal@Series[%(cpp-cpb),{sp,0,4},{t,0,4}]/.{sp->sp ep,t->t ep},{ep,0,3}]/.ep->1,{sp,t},FullSimplify];

(CoefficientList[light1-heavy1,{sp,t}]//Flatten)//.{a___,0,b___}:>{a,b}
rep1=Table[Solve[%[[i]]==0,Take[listwhat[%[[i]],i_g],1]]//Flatten,{i,1,Length@%}]//Simplify//DeleteDuplicates

(CoefficientList[light2-heavy2,{sp,t}]//Flatten)//.{a___,0,b___}:>{a,b}
rep2=Table[Solve[%[[i]]==0,Take[listwhat[%[[i]],i_g],1]]//Flatten,{i,1,Length@%}]//Simplify//DeleteDuplicates

(CoefficientList[light3-heavy3,{sp,t}]//Flatten)//.{a___,0,b___}:>{a,b}
rep3=Table[Solve[%[[i]]==0,Take[listwhat[%[[i]],i_g],1]]//Flatten,{i,1,Length@%}]//Simplify//DeleteDuplicates

rep1
rep2

rep3

A smarter way to count null-constraints. Note: my code is not general and compact enough, make your own better codes!

Flatten@CoefficientList[light1,{sp,t}]~Join~Flatten@CoefficientList[light2,{sp,t}]~Join~Flatten@CoefficientList[light3,{sp,t}];
pos0=Position[%,0];
Delete[%%,pos0]
mattonulls=Reverse@NullSpace[Coefficient[%,#]&/@listwhat[%,i_g]];

Flatten@CoefficientList[heavy1,{sp,t}]~Join~Flatten@CoefficientList[heavy2,{sp,t}]~Join~Flatten@CoefficientList[heavy3,{sp,t}];
Delete[%,pos0];
tocount=(mattonulls.%//Simplify)//.{a___,0,b___}:>{a,b}//DeleteDuplicates
tocountm4=Coefficient[tocount,1/m^4]//.{a___,0,b___}:>{a,b}//DeleteDuplicates
tocountm6=Coefficient[tocount,1/m^6]//.{a___,0,b___}:>{a,b}//DeleteDuplicates
tocountm8=Coefficient[tocount,1/m^8]//.{a___,0,b___}:>{a,b}//DeleteDuplicates

tocountm6/(Coefficient[tocountm6/.J->0,cpb])//DeleteDuplicates
%//Length

tocountm8//Simplify
%//Length

(*check 3 are independent*)

Table[c[i],{i,1,3}].tocountm8;
Flatten@CoefficientList[%,{cpb,cpp,J}];
Solve[%==0]//Flatten

problem 4

eyeball work

lightpart={g[2],g[3],g[4],0};
heavypart={1/m^4,(3-2(J(J+1)))/m^6,1/(2m^8),((J (1+J) (-8+J+(J^2)))/( (m^8)))};

Let’s take M=1 for simplicity. It can be easily restored.

prove (Subscript[g, 4])<=(Subscript[g, 2])/(2M^4)

lightpart[[1]]/2-lightpart[[3]]
heavypart[[1]]/2-heavypart[[3]]//Simplify
Plot[% m^8,{m,1,100}]

prove (Subscript[g, 3])<=3(Subscript[g, 2])/M^2

3lightpart[[1]]-lightpart[[2]]
3heavypart[[1]]-heavypart[[2]]//Simplify
Table[{J,m,% m^6},{J,0,40,2},{m,1,8,0.01}]//Flatten[#,1]&;
ListContourPlot[%,PlotLegends->Automatic]

linear programming on Mathematica

prove (Subscript[g, 4])<=(Subscript[g, 2])/(2M^4). We take our basis to be {g[2],g[4],0} where the third one is null constraint.

norm={0,1,0};
obj={-1,0,0};
mat=Table[Table[heavypart[[i]],{J,0,40,2},{m,1,20,0.1}]//Flatten[#,1]&,{i,{1,3,4}}]//Transpose;//AbsoluteTiming

LinearOptimization[-obj,{mat,0*mat[[;;,1]]},{{norm},{1}}]
(%.{g[2],g[4],0})>=0

prove (Subscript[g, 3])<=3(Subscript[g, 2])/M^2. We take our basis to be {g[2],g[3],0}

norm={0,1,0};
obj={-1,0,0};
mat=Table[Table[heavypart[[i]],{J,0,40,2},{m,1,20,0.1}]//Flatten[#,1]&,{i,{1,2,4}}]//Transpose;//AbsoluteTiming

LinearOptimization[-obj,{mat,0*mat[[;;,1]]},{{norm},{1}}]
(%.{g[2],g[3],0})>=0

prove (Subscript[g, 3])>=-10.61249(Subscript[g, 2])/M^2. We take our basis to be {g[2],g[3],0}

norm={0,-1,0};
obj={-1,0,0};
mat=Table[Table[heavypart[[i]],{J,0,40,2},{m,1,40,0.01}]//Flatten[#,1]&,{i,{1,2,4}}]//Transpose;//AbsoluteTiming

LinearOptimization[-obj,{mat,0*mat[[;;,1]]},{{norm},{1}}]
(%.{g[2],g[3],0})>=0

add (Subscript[n, 5]) and repeate

n5=1/m^10((23d^2-12d-20)JJsq+(-21d-2)JJsq^2+4JJsq^3)/.{d->4,JJsq->J(J+1)};
heavypartadd=heavypart~Join~{n5}

norm={0,1,0,0};
obj={-1,0,0,0};
mat=Table[Table[heavypartadd[[i]],{J,0,40,2},{m,1,20,0.1}]//Flatten[#,1]&,{i,{1,3,4,5}}]//Transpose;//AbsoluteTiming

LinearOptimization[-obj,{mat,0*mat[[;;,1]]},{{norm},{1}}]
(%.{g[2],g[4],0,0})>=0

norm={0,1,0,0};
obj={-1,0,0,0};
mat=Table[Table[heavypartadd[[i]],{J,0,40,2},{m,1,20,0.1}]//Flatten[#,1]&,{i,{1,2,4,5}}]//Transpose;//AbsoluteTiming

LinearOptimization[-obj,{mat,0*mat[[;;,1]]},{{norm},{1}}]
(%.{g[2],g[3],0,0})>=0

norm={0,-1,0,0};
obj={-1,0,0,0};
mat=Table[Table[heavypartadd[[i]],{J,0,40,2},{m,1,40,0.01}]//Flatten[#,1]&,{i,{1,2,4,5}}]//Transpose;//AbsoluteTiming

LinearOptimization[-obj,{mat,0*mat[[;;,1]]},{{norm},{1}}]
(%.{g[2],g[3],0,0})>=0

using sdpb

below is SDPB wrapper, taken from Simon. You will need SDPB.m and you will also need to install docker.

For windows, you have to replace \->//

SetDirectory[NotebookDirectory[]];

<<"SDPB.m"
dockerdir=Directory[];
solvesdp[sdp_(*SDP[objective,norm,positiveMatrices]*),xmlfile_:"sdp.xml"]:=(WriteBootstrapSDP[xmlfile,sdp](*write down xml file*);RunProcess[{(*"powershell","docker"*)(*this is for Windows Users*)"/Applications/Docker.app/Contents/Resources/bin/docker"(*this is for Mac Users*),"run","-v",StringJoin[{dockerdir,":/data"}],"wlandry/sdpb:2.4.0","pvm2sdp","1024",StringJoin[{"/data/",xmlfile}],StringJoin[{"/data/",StringReplace[xmlfile,".xml"->""],"/"}]},"StandardOutput"](*convert xml file to sdp*);{RunProcess[{(*"powershell","docker"*)"/Applications/Docker.app/Contents/Resources/bin/docker","run","-v",StringJoin[{dockerdir,":/data"}],"wlandry/sdpb:2.4.0","sdpb","--precision=2048","--procsPerNode=1"(*,"--findDualFeasible"*),"-s"(*,"--maxIterations","1500"*),StringJoin[{"/data/",StringReplace[xmlfile,".xml"->""],"/"}]},"StandardOutput"](*run sdpb*),StringReplace[dockerdir~~"/"~~xmlfile,".xml"->"_out"](*return the output directory*)})
(* read the output, and return the extremal functional y. You need the 'norm' part of the sdp to reconstruct the full vector. *)(*borrow from simon*)
getFunctional[outDir_,norm_]:=Module[{α,yVector},
(* first, get the termination reason, which is the first line of out.txt *)
Clear[terminateReason];
ToExpression[ReadList[outDir~~"/out.txt",String][[1]]];
Print[terminateReason];
(* then get the extremal functional. the output file y.txt contains a newline-separated list of numbers, with e for 10^ *)
(* the first line is not useful, drop it. *)
ytmp=ToExpression["{"~~StringReplace[Import[outDir~~"/y.txt","Text"],{"e"->"*^","\n"->","}]~~"}"][[2;;]];
(* the vector has one component removed using y.norm==1. Recover it. *)
yVector = Insert[ytmp,α,maxIndexBy[norm,Abs]];
yVector/.Solve[Collect[yVector.norm,α,Expand]==1,α][[1]]]

(Subscript[g, 4])<= 1/2 (Subscript[g, 2])

norm={0,-1,0};
obj={-1,0,0};
mat=Table[{{Table[heavypart[[i]],{i,{1,3,4}}] m^8/.{m->1+x}}},{J,0,40,2}];

positiveMatrices=PositiveMatrixWithPrefactor[DampedRational[1,{},1/E,x],#]&/@mat;

(*after you run one, you have to delete xml file and checking points, otherwise it interferes with the next*)

DeleteFile[Directory[]~~"/dosdpb.xml"]
DeleteFile[Directory[]~~"/dosdpb.ck"~~"/checkpoint.json"]

outDir =solvesdp[SDP[obj,norm,positiveMatrices],"dosdpb.xml"];//AbsoluteTiming

(*this is output*)

outDir[[1]]

getFunctional[outDir[[2]],norm]
(%.{g[2],g[4],0})>=0

(Subscript[g, 3])<=3(Subscript[g, 2])/M^2

norm={0,-1,0};
obj={-1,0,0};
mat=Table[{{Table[heavypart[[i]],{i,{1,2,4}}] m^8/.{m->1+x}}},{J,0,40,2}];

positiveMatrices=PositiveMatrixWithPrefactor[DampedRational[1,{},1/E,x],#]&/@mat;

(*after you run one, you have to delete xml file and checking points, otherwise it interferes with the next*)

DeleteFile[Directory[]~~"/dosdpb.xml"]
DeleteFile[Directory[]~~"/dosdpb.ck"~~"/checkpoint.json"]

outDir =solvesdp[SDP[obj,norm,positiveMatrices],"dosdpb.xml"];//AbsoluteTiming

getFunctional[outDir[[2]],norm]
(%.{g[2],g[3],0})>=0

(Subscript[g, 3])>=-10.61249(Subscript[g, 2])/M^2

norm={0,1,0};
obj={-1,0,0};
mat=Table[{{Table[heavypart[[i]],{i,{1,2,4}}] m^8/.{m->1+x}}},{J,0,40,2}];

positiveMatrices=PositiveMatrixWithPrefactor[DampedRational[1,{},1/E,x],#]&/@mat;

(*after you run one, you have to delete xml file and checking points, otherwise it interferes with the next*)

DeleteFile[Directory[]~~"/dosdpb.xml"]
DeleteFile[Directory[]~~"/dosdpb.ck"~~"/checkpoint.json"]

outDir =solvesdp[SDP[obj,norm,positiveMatrices],"dosdpb.xml"];//AbsoluteTiming

getFunctional[outDir[[2]],norm]
(%.{g[2],g[3],0})>=0

problem 5

We consider ansatz that f is a polynomial in p: (1-p^2)^2(1+ p^2 + p^3 )+...

ruleint={(1-p^2)^2p^(a_:1)Cos[p b]:>((8 HypergeometricPFQ[{(1/2)+(a/2)},{(1/2),(7/2)+(a/2)},-((b^2)/4)])/((1+a) (3+a) (5+a))),(1-p^2)^2Cos[p b]:>(((8 HypergeometricPFQ[{(1/2)+(a/2)},{(1/2),(7/2)+(a/2)},-((b^2)/4)])/((1+a) (3+a) (5+a)))/.a->0)};

make basis of our functional. Let’s truncate to n=5.

1/Pi Table[(1-p^2)^2 p^(i),{i,{0,2,3,4,5}}]Cos[p b];
basfb=%/.ruleint//Expand;

norms=Table[Integrate[basfb[[i]],{b,-Infinity,Infinity}],{i,1,Length[basfb]}];//AbsoluteTiming

norms

width=Table[Integrate[basfb[[i]]b^2//Expand,{b,-Infinity,Infinity}],{i,1,Length[basfb]}];//AbsoluteTiming

matodo=Table[{{basfb b^10//Expand}}/.b->x,{x,-200,200,1/20}]//N;//AbsoluteTiming
matodolargeb=Table[{{Normal@Series[basfb b^3,{b,Infinity,0}]}}/.b->x,{x,{-300,300,-400,400,-1000,1000,-2000,2000}}]//N;//AbsoluteTiming
objtodo=-width;
normtodo=norms;
positiveMatrices=PositiveMatrixWithPrefactor[DampedRational[1,{},1/E,x],#]&/@(matodo~Join~matodolargeb);

DeleteFile[Directory[]~~"/prob5.xml"]
DeleteFile[Directory[]~~"/prob5.ck"~~"/checkpoint.json"]

outDir =solvesdp[SDP[objtodo,normtodo,positiveMatrices],"prob5.xml"];//AbsoluteTiming

getFunctional[outDir[[2]],normtodo];
-%.objtodo

Pi^2//N

(*we can improve it to make more and more accurate by considering more refined sample points of b, but *)

test positivity

getFunctional[outDir[[2]],normtodo];
toplot=%.basfb;

Plot[toplot,{b,-100,100},WorkingPrecision->100]