(* First, I compute the tt-component of the 2D stress-enetgy tensor numerically. (This is the method provided in the paper)*)
I1[us_]:=(1/4)*(I/(8*(Pi^3)))*NIntegrate[ Exp[-2*Pi*y]*Gamma[1+4*I*y]/((4*y)^(4*I*y))NIntegrate[Exp[-2*Pi*x]*Gamma[1-4*I*x]*((4*x)^(4*I*x))*Exp[I*us*(y-x)]/(x-y), {x,0,y,Infinity},Method->{"PrincipalValue","SingularPointIntegrationRadius"->1/10}], {y,0,Infinity}]

IDirac=(1/4)*(-1/(8*Pi^2))*NIntegrate[ Exp[-4*Pi*y]*Gamma[1+4*I*y]*Gamma[1-4*I*y],{y,0,Infinity}]

I2[us_]:=(1/4)*(I/(8*(Pi^3)))*NIntegrate[ Exp[-2*Pi*y]*Gamma[1-4*I*y]*((4*y)^(4*I*y))*Exp[-2*Pi*x]*Gamma[1-4*I*x]*((4*x)^(4*I*x))*Exp[-I*us*(y+x)]/(x+y), {x,0,Infinity}, {y,0,Infinity}];

(* This is equation 6.29 of the our paper. *)

Itot[us_]:=I1[us]+I2[us]+(-0.00020723299881758508`-1.2502718732455542`*^-24 I)

p1table=Table[{j,Re[Itot[j]]},{j,0,20}]

{{0,-0.00018986468681910176`},{1,-0.0001626524701290107`},{2,-0.00013567939266059052`},{3,-0.00010998055254791695`},{4,-0.00008651421116921893`},{5,-0.00006601250048808712`},{6,-0.0000488801564419963`},{7,-0.00003517123013101301`},{8,-0.000024642624825468423`},{9,-0.00001685576544310216`},{10,-0.000011287818461655399`},{11,-7.422263558573984`*^-6},{12,-4.8055741631620075`*^-6},{13,-3.0715522453837523`*^-6},{14,-1.9425344294973226`*^-6},{15,-1.2179724230397407`*^-6},{16,-7.583842863747843`*^-7},{17,-4.6959911756416523`*^-7},{18,-2.894980926196908`*^-7},{19,-1.7784597715285433`*^-7},{20,-1.0895425155569628`*^-7}}

p1tablep={{0,-0.00018986468681910176`},{1,-0.0001626524701290107`},{2,-0.00013567939266059052`},{3,-0.00010998055254791695`},{4,-0.00008651421116921893`},{5,-0.00006601250048808712`},{6,-0.0000488801564419963`},{7,-0.00003517123013101301`},{8,-0.000024642624825468423`},{9,-0.00001685576544310216`},{10,-0.000011287818461655399`},{11,-7.422263558573984`*^-6},{12,-4.8055741631620075`*^-6},{13,-3.0715522453837523`*^-6},{14,-1.9425344294973226`*^-6},{15,-1.2179724230397407`*^-6},{16,-7.583842863747843`*^-7},{17,-4.6959911756416523`*^-7},{18,-2.894980926196908`*^-7},{19,-1.7784597715285433`*^-7},{20,-1.0895425155569628`*^-7}}

p1=ListPlot[p1tablep, PlotStyle->{Black,Thickness[0.001]}]

p2=Plot[1/(24*Pi)((-8/(-4-4*ProductLog[Exp[(-x)/4]])^3)-24/((-4-4*ProductLog[Exp[(-x)/4]])^4))-1/(768*Pi), {x,0,20},PlotStyle->{Black,Thickness[0.001]}]

Show[p1,p2,Frame->True,FrameStyle->Directive[Black,Thickness[0.003]],LabelStyle->{FontFamily->"Times" ,Black},Axes->False,FrameLabel->{"(v/M)",""}]