# Analytical form of 2d integrals relevant to graphene

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.

• Why do you have two different cases xiN marked the same? It is necessary to separate identifiers, for example,xiN and xiN1. Undefined {Δxi[k,1,c,1],Δxi[k,-1,c,1],Δxi[k,1,c,1]+Δxi[k,-1,c,1]} therefore not reproduced Plot Commented Sep 8, 2019 at 17:27
• @AlexTrounev I've corrected the definitions. I was simplifying the things and forgot to remove kappa and to add the definition of Δxi. As for the xiN, they are not the same. The 1st line defines the function for s=+1 or -1. The second line with the 3rd parameter equal to 0 defines the sum of two cases, +1 and -1. This is obtain by simplifying the sum of integrands. Commented Sep 8, 2019 at 17:55
• Is h2 supposed to have a big discontinuity at x=1? This must be the source of the discontinuities at k=30... Commented Sep 8, 2019 at 23:10
• There is a problem with the analytical solution. Therefore, we need to get this solution. Commented Sep 9, 2019 at 4:50
• @MelaGo If we assume that the sum of Δxi[k,1,c,1]+Δxi[k,-1,c,1] should be zero, which my calculations seem to show (at variance with the paper), then there should be no f term in the expression for Δxi, h2 should be modified accordingly and it should be continuous. The reason for the separation of xe and xi into f and h is to have s-dependent and s-independent parts. Commented Sep 9, 2019 at 7:34

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}}$$

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.

• What you've used here are basically the quadratic AGM transformations for the complete elliptic integrals. Commented Nov 13, 2019 at 8:52
• @J.M. Do you know what AGM means? Also, people now are forgetting much about the elliptic integrals. Can you recommend some nice book on the topic? Commented Nov 13, 2019 at 8:56
• By "AGM" I meant the classical arithmetic-geometric mean of Gauss. The most useful book that readily comes to mind was the one I linked to in my earlier answer, which is Byrd and Friedman's; if I think of others, I'll send them your way. Commented Nov 13, 2019 at 9:01
• I also should hasten to add something about your remark, "Notice that for $k>1$ an imaginary part appears. Thus, we will be needing only the real part." You can use the formulae here so that you don't need to get rid of imaginary parts, which yields the replacement rule EllipticE[(4 λ_)/(1 + λ_)^2] :> 2 λ EllipticE[1/λ^2]/(1 + λ) + (1/λ - 1) EllipticK[1/λ^2]. Commented Nov 13, 2019 at 9:03
• @J.M. I see, it's very useful. Commented Nov 13, 2019 at 9:07