Analytical form of 2d integrals relevant to graphene

Question

This question is continuation of my previous post. Alex Trounev was very helpful in fixing a crucial typo in the analytic solution known from the article "Density Dependent Exchange Contribution to ∂𝜇/∂𝑛 and Compressibility in Graphene". Here, the question is about obtaining the analytic solution in fully automatic way with MA. It consists of 2 parts:

1. I know the analytic solution for xe (but do not know how to obtain it) and I can verify it numerically.

2. I have numerical solution for xi and a fixed analytical solution that works in some parameter range.

   I am asking because I feel that there are other typos that I am not able to fix myself.

---

We will be needing some auxiliary functions [Eqs.(5,6) in the paper] that I define beforehands

(*i[1]=Integrate[EllipticK[y^2]-EllipticE[y^2],{y,0,1/x},Assumptions->x>1]*)
i[1]=1/(2 x) π (-HypergeometricPFQ[{-(1/2),1/2,1/2},{1,3/2},1/x^2]+HypergeometricPFQ[{1/2,1/2,1/2},{1,3/2},1/x^2]);

(*i[2]=Integrate[(EllipticK[y^2]-EllipticE[y^2]-π/4 y^2)/y^3,{y,0,x},Assumptions->0≤x≤1]*)
i[2]=3/256 π x^2 (HypergeometricPFQ[{1,1,3/2,5/2},{2,3,3},x^2]+3 HypergeometricPFQ[{1,1,5/2,5/2},{2,3,3},x^2]);

f[x_] =Piecewise[{{EllipticE[x^2],x<=1},{x EllipticE[1/x^2]-(x-1/x)EllipticK[1/x^2],True}}];
h[x_] =Piecewise[{{x(π/4Log[4/x]-π/8)-x i[2],x<=1},{x i[1],x>1}}];
h2[x_]=Piecewise[{{x(π/4Log[4/x]+π/8)-x i[2],x<=1},{x i[1],x>1}}];

These expressions I want to reproduce with MA

Calculation of xe

With the definitions

p[k_,q_,θ_]:=Sqrt[k^2+q^2-2 k q Cos[θ]]
Fs[k_,q_,θ_,si_,sf_]:=1/2(1+si sf (k-q Cos[θ] )/p[k,q,θ] )
n[k_,s_]:=UnitStep[1- k s]
n0[s_]:=1/2(1-s)
δn[k_,s_]:=n[k,s]-n0[s]

we get the first 2d integral that I want to know how to compute analytically

xeN[k_?NumericQ,s_]:=-(1/(2π))NIntegrate[
Fs[k,q,θ,s,1]δn[p[k,q,θ],1]+Fs[k,q,θ,s,-1]δn[p[k,q,θ],-1],
{q,0,∞},{θ,0,2π},PrecisionGoal->2]

xeN[k_?NumericQ,0]:=-(1/(2π))NIntegrate[UnitStep[1- p[k,q,θ]],
{q,0,∞},{θ,0,2π},PrecisionGoal->4]  (* sum of the 1 and -1 cases*)

Here, xe is the known analytical expression for xeN but how to get this with MA?:

xe[k_,s_]:=1/π (-f[k]-s h[k])

Verification

is time consuming

eval[ 1]=Table[{k,xeN[k, 1]},{k,0.1,2,0.2}]
eval[-1]=Table[{k,xeN[k,-1]},{k,0.1,2,0.2}]
eval[ 0]=Table[{k,xeN[k, 0]},{k,0.1,2,0.2}];   (*the sum of the 1 and -1 cases*)

and plot

Plot[{xe[k,1],xe[k,-1],xe[k,1]+xe[k,-1]},{k,0,2},
PlotRange->{-1.05,0},PlotTheme->{"FrameGrid","BoldColor"},FrameStyle->14,FrameLabel->{"k","Subscript[x, e]"},Epilog->{PointSize[Medium],Blue,Point[eval[1]],Orange,Point[eval[-1]],Red,Point[eval[0]]}]

This is actually the Fig.1 from the paper.

Calculation of xi

The second integral is simpler from the numerical point of view

xiN[k_?NumericQ,kc_?NumericQ,s_]:=-(1/(2π))NIntegrate[
Fs[k,q,θ,s,1]n0[1]+Fs[k,q,θ,s,-1]n0[-1],
{q,0,kc},{θ,0,2π},PrecisionGoal->4]+kc/2

xiN[k_?NumericQ,kc_?NumericQ,0]:=-(1/(2π))NIntegrate[1,
{q,0,kc},{θ,0,2π},PrecisionGoal->4]+kc   (* sum of the 1 and -1 cases*)

but not so simple analytically...

Verification

c=30;
ival[1] =Table[{k,xiN[k,c,1]}, {k,0.1,40,2}];
ival[-1]=Table[{k,xiN[k,c,-1]},{k,0.1,40,2}];
ival[0] =Table[{k,xiN[k,c,0]}, {k,0.1,40,2}];   (*the sum of the 1 and -1 cases*)

against analytical solution

xi[k_,s_,kc_]:=kc/π (-f[k/kc]+s h2[k/kc])
Δxi[k_,s_,kc_]:=xi[k,s,kc]+kc/2

and plot

Plot[{Δxi[k,1,c],Δxi[k,-1,c],Δxi[k,1,c]+Δxi[k,-1,c]},{k,0,40},
PlotRange->{-20,20},PlotTheme->{"FrameGrid","BoldColor"},FrameStyle->14,FrameLabel->{"k","Subscript[x, i]"},Epilog->{PointSize[Medium],Blue,Point[ival[1]],Orange,Point[ival[-1]],Red,Point[ival[0]]}]

Thanks to the Alex Trounev answer, h was replaced with h2 (sign correction) thereby fixing the small-$k$ behavior, however, there are still problems at high $k$ as can be seen from the plot.

Answer 1

I managed to find an answer to my question and would like to share my experience with the community. It turns out that MA can do a lot with elliptic integrals, however, it does not know some identities, in particular Eq.8.126 from Table Of Integrals, Series And Products by Gradshteyn and Ryzhik. I retype them below in original notations $$K\!\left(\frac{2\sqrt{k}}{1+k}\right)=(1+k)K(k),$$ and $$E\!\left(\frac{2\sqrt{k}}{1+k}\right)=\frac1{1+k}\left(2E(k)-(1-k^2)K(k)\right).$$ One should remember that they are valid for $0\le k\le 1$ and one should beware of slightly different MA conventions. We will be using the following rules

rules={EllipticK[(4 λ_)/(1+λ_)^2]-> (1+λ)EllipticK[λ^2],
       EllipticE[(4 λ_)/(1+λ_)^2]-> 1/(1+λ) (2EllipticE[λ^2]-(1-λ^2)EllipticK[λ^2])};

for simplification. Notice that for $k>1$ an imaginary part appears. Thus, we will be needing only the real part.

Numerical verification

g[1]=Plot[ EllipticK[(4 k )/(k+1)^2],{k,0,3},
          PlotRange->All,PlotTheme->{"Monochrome","Frame"},
Mesh->21,MeshStyle->{Red,PointSize[Medium]},PlotStyle->None];
g[2]=Plot[Re[(1+k)EllipticK[k^2]],{k,0,3},
          PlotRange->All,PlotTheme->{"Monochrome","Frame"}];
g[5]=Show[g[2],g[1],FrameLabel->{"k","K(k^2)"},LabelStyle->14];
g[3]=Plot[ EllipticE[(4 k )/(k+1)^2],{k,0,3},
           PlotRange->All,PlotTheme->{"Monochrome","Frame"},
Mesh->21,MeshStyle->{Red,PointSize[Medium]},PlotStyle->None];
g[4]=Plot[Re[1/(1+k) (2EllipticE[k^2]-(1-k^2)EllipticK[k^2])],{k,0,3},
           PlotRange->All,PlotTheme->{"Monochrome","Frame"}];
g[6]=Show[g[4],g[3],FrameLabel->{"k","E(k^2)"},LabelStyle->14];
g[7]=GraphicsRow[{g[5],g[6]},ImageSize->800]

Now let us consider one particular integral that was problematic before, namely xi. Other integrals can be computed in a similar way.

Computation of xi

Taking into account the definitions above, the integrand reads

ui=Simplify[Fs[k,q,θ,s,1]n0[1]+Fs[k,q,θ,s,-1]n0[-1]]
Δui=ui-1/2//Simplify

$$\frac12 \left[ 1-\frac{s (k-q \cos\theta)}{\sqrt{k^2+q^2-2 k q \cos\theta}}\right]$$

The constant term is trivial to integrate, we turn to the second one

$$\delta x_i=\frac1{4\pi}\int_0^{k_c}\mathrm{d}q\int_0^{2\pi}\mathrm{d}\theta\frac{s (k-q \cos\theta)}{\sqrt{k^2+q^2-2 k q \cos\theta}}$$

The angular integral reads

ri[1]=2 Integrate[Δui,{θ,0,π},Assumptions->k>0&&q>0&&k!=q]
(*-((s ((k+q) EllipticE[(4 k q)/(k+q)^2]+(k-q) EllipticK[(4 k q)/(k+q)^2]))/k)*)

now using the replacement rules

ri[2]=(((ri[1]/.q->x k)//Simplify)/.rules)//Simplify
(* -2 s EllipticE[x^2] *)

the integral takes a very simples form. We do the second integration discarding the imaginary part in accordance with the analysis above

ri[3]=Integrate[ri[2]/.{x-> q/k},{q,0,kc},Assumptions->kc>k>0]
(* -(1/32) k π s (-((k^2 HypergeometricPFQ[{1,1,3/2,3/2},{2,2,3},k^2/kc^2])/kc^2)+8 (1+Log[16]-2 Log[k]+2 Log[kc]-I MeijerG[{{},{1/2,1/2,1}},{{-1,0,0},{}},k^2/kc^2])) *)

ri[4]=ri[3]/.MeijerG[_,_,_]->0
(* -(1/32) k π s (-((k^2 HypergeometricPFQ[{1,1,3/2,3/2},{2,2,3},k^2/kc^2])/kc^2)+8 (1+Log[16]-2 Log[k]+2 Log[kc])) *)

ri[5]=Integrate[ri[2]/.{x-> q/k},{q,0,kc},Assumptions->k/kc>1&&kc>0]
(* -kc π s HypergeometricPFQ[{-(1/2),1/2,1/2},{1,3/2},kc^2/k^2]*)

Numerical verification

xiN[k_?NumericQ,kc_?NumericQ,s_]:=-(1/(2π))NIntegrate[Fs[k,q,θ,s,1]n0[1]+Fs[k,q,θ,s,-1]n0[-1],{q,0,kc},{θ,0,2π},PrecisionGoal->4]+kc/2
Δxi[k_,kc_,s_]:=Piecewise[{{(k s)/4 (1/2 -(k/(4 kc))^2 HypergeometricPFQ[{1,1,3/2,3/2},{2,2,3},k^2/kc^2]+Log[(4 kc)/k]),k<=kc},{(kc s )/2 HypergeometricPFQ[{-(1/2),1/2,1/2},{1,3/2},kc^2/k^2],k>kc}}]

c=30;
ival[1]=Table[{k,xiN[k,c,1]},{k,0.1,40,2}];
ival[-1]=Table[{k,xiN[k,c,-1]},{k,0.1,40,2}];

Plot[{Δxi[k,c,1],Δxi[k,c,-1]},{k,0,40},PlotRange->{-20,20},
PlotTheme->{"FrameGrid","Monochrome","BoldColor"},
FrameStyle->14,FrameLabel->{"k","Subscript[x, i]"},Epilog->{PointSize[Medium],Blue,Point[ival[1]],Orange,Point[ival[-1]]}]

Thus, analytic expression for Δxi is the main result here.