Tackle a general elliptic integral?

Questions
1

I’m currently trying to reproduce some of the calculations in Phys Rev B, 18 4509 (1978), doi:10.1103/PhysRevB.18.4509. You shouldn’t need to read the paper to understand my problem though.

Specifically, I’m interested in evaluating the following integral (from Appendix B in the paper):

G_1=\frac{1}{2\Delta p}\left[1-\frac{i}{\pi}\int_{-1}^1 \sqrt{\frac{1-\left(x+\frac{\epsilon }{2}\right)^2}{1-x^2}}\frac{dx}{\Omega +x}\right]

\epsilon and \Omega are both real but are not otherwise restricted in their values except they are interdependent and come from the physical model parameters. \Delta p is another related parameter but isn’t involved in the integral itself. Obviously, depending on their values, the integrand itself will become complex at certain points in the interval, and/or the integral will no longer converge. In the latter case the integral is imaginary but is physically meaningful, since the imaginary part of the integral is what contributes to the electronic density of states in this model (the main thing I’m after). The authors used Byrd’s Handbook of Elliptic Integrals to arrive at the following expression for G_1 (Mathematica’s Integrate was not able to express the integral in terms of fundamental elliptic integrals):

G_1= \left\{\begin{array}{lr} \left(\frac{1}{2\Delta p}\right)\left\{1+\left(\frac{2}{\pi}\right)\left(1+\xi\right)\left[\alpha\Pi\left(-\beta\xi, \left|\xi\right|\right)-\Pi\left(-\xi,\left|\xi\right|\right)\right]\right\}, & \left|\epsilon\right|>4\\ \left(\frac{1+\eta}{\pi\Delta p}\right)\text{sgn}\left(\epsilon\right)\left[\alpha\Pi\left(-\beta\eta, \left|\eta\right|\right)-\Pi\left(-\eta,\left|\eta\right|\right)\right]+\text{sgn}\left(\gamma\right)G_2-\left(\frac{i\eta}{\pi\Delta p}\right)\left[\beta\Pi\left(-\alpha\left(1+\eta\right), \left|\eta'\right|\right)-\Pi\left(1+\eta,\left|\eta'\right|\right)\right], & \left|\epsilon\right|<4\\ \end{array}\right.

Where

\alpha=(\Delta p + 1)^2/(\Omega+1)
\beta=(\Delta p - 1)^2/(\Omega+1)
\gamma=(\Delta p^2 - 1)/\Delta p
\eta=\xi^{-1}=\frac{1}{4}\epsilon
\eta'=1-\eta^2
G_2=\frac{\Delta p^2 - 1}{4\Delta p^2\sqrt{\Omega^2-1}}

I’ve tried evaluating the above expression and it is not successfully reproducing the results in the paper. Since I’m not super-familiar with elliptic integrals I was hoping to use Mathematica in a way to verify the above expression and find the cause of the problem, although I am concurrently trying to use the handbook to arrive at the expression by hand. As I understand it there are algorithms for expressing general elliptic integrals in terms of the fundamental ones but I’m not sure how/if they’re implemented in the Integrate method within Mathematica or if there’s a different, special approach I could use.

Thanks!

2

Welcome to Wolfram Community! Please make sure you know the rules: https://wolfr.am/READ-1ST Please post your Wolfram Language code.

3

Here are some things I’ve tried:

Integrate[(1/(x + \[CapitalOmega])) ((1 - (x + \[Epsilon]/2)^2)/(
   1 - x^2))^(1/2), {x, -1, 1}]

Which doesn’t go anywhere. By doing some manipulations by hand, getting the radical completely in the denominator, then doing polynomial division on the term outside, I can get three separate integrals that are listed in a table of integrals. The simplest of these is the following:

Integrate[1/((1 - x) (-1 - x) (1 - \[Epsilon]/2 - x) (-1 - \[Epsilon]/
     2 - x))^(1/2), {x, -1, 1}]

Which also doesn’t result in an answer. However, replacing the variables with explicit values gives an actual answer:

In[1]:= Integrate[1/((1-x)(-1-x)(3-x)(-3-x))^(1/2),{x,-1,1}]
Out[1]= 2/3 EllipticK[1/9]

As soon as the roots aren’t from difference of squares terms though, it doesn’t work again:

Integrate[1/((1 - x) (-1 - x) (2 - x) (-3 - x))^(1/2), {x, -1, 1}]

I think I’m just going to have to simplify the entire thing by hand to get it in terms of EllipticK and EllipticPi, then evaluate it at that point.

4

Assumptions on epsilon and omega will generally help: Are they real? positive? negative? et cetera… By default Mathematica will assume everything can be complex which might be ‘overkill’ and does allow for simple answers.

As a first step I would verify the results the author found: try some ‘random’ values and check if the results are in agreement (exact and numerical integration).

5

Thanks for your reply. Omega and epsilon are real. Due to the complexity of the integral adding in the assumptions didn’t change anything (the kernel just crashes after a while).

I’ve already done some attempts at verifying the authors’ results. In fact, the lack of agreement was why I’m trying to re-do their solving of this integral to begin with. Their formula for G_1 only works in one situation. To explain more I need to write some more definitions:

\epsilon=(\omega-\epsilon_t)(\omega-\epsilon_p)-4
\Delta p=1-(\omega-\epsilon_p)\Delta\epsilon_t
\Omega=\frac{\Delta p^2 + 1}{2\Delta p} + \frac{\epsilon}{2}

Here, the true fundamental parameters are \epsilon_t, the energy of a certain orbital (diagonal Hamiltonian matrix element), \epsilon_p is another orbital energy, and \Delta\epsilon_t is a perturbation of a surface ion’s energy due to the change in atomic environment. The formula for G_1 only works when \Delta\epsilon_t=0, and only gives the correct value for the imaginary part (contributing to the density of states). The real part disagrees with what is plotted in the paper. Similarly though, solving the integral for G1 numerically is pretty unstable and I’m unable to get agreement with the paper except in the \Delta\epsilon_t=0 case for the imaginary part.

Here’s some code showing some more playing around with this system:

\[Epsilon]p = -9.0;
\[Epsilon]t = -6.0;
\[CapitalDelta]\[Epsilon]t = 0.0;
\[Epsilon][w_] := (w - \[Epsilon]p) (w - \[Epsilon]t) - 4
\[CapitalDelta]p[w_] := 
 1 - (w + \[Epsilon]p) \[CapitalDelta]\[Epsilon]t
\[CapitalOmega][w_] := ((\[CapitalDelta]p[w])^2 + 1)/(
  2 \[CapitalDelta]p[w]) + \[Epsilon][w]/2
Plot[\[Epsilon][w], {w, -12, -3}]
Plot[\[CapitalDelta]p[w], {w, -12, -3}]
Plot[\[CapitalOmega][w], {w, -12, -3}]

Manipulate[
 Plot[{Re[(1/(x + \[CapitalOmega][
          w])) ((1 - (x - \[Epsilon][w]/2)^2)/(1 - x^2))^(1/2)], 
   Im[(1/(x + \[CapitalOmega][w])) ((1 - (x - \[Epsilon][w]/2)^2)/(1 -
           x^2))^(1/2)]}, {x, -1, 1}], {w, -12.5, -3}]

\[Alpha][w_] := (\[CapitalDelta]p[w] + 1)^2/(\[CapitalOmega][w] + 1)
\[Beta][w_] := (\[CapitalDelta]p[w] - 1)^2/(\[CapitalOmega][w] + 1)
\[Gamma][w_] := ((\[CapitalDelta]p[w])^2 - 1)/\[CapitalDelta]p[w]
\[Eta][w_] := 1/4 \[Epsilon][w]
\[Eta]1[w_] := 1 - (\[Eta][w])^2
\[Xi][w_] := 1/\[Eta][w]
G2[w_] := ((\[CapitalDelta]p[w])^2 - 1)/(
 4 (\[CapitalDelta]p[w])^2 Sqrt[(\[CapitalOmega][w])^2 - 1])

t = -4.7; N[
 Piecewise[{{1/(
     2 \[CapitalDelta]p[t]) (1 + 
       2/Pi (1 + \[Xi][
           t]) (\[Alpha][
            t] EllipticPi[-\[Beta][t] \[Xi][t], (Abs[\[Xi][t]])^2] - 
          EllipticPi[-\[Xi][t], (Abs[\[Xi][t]])^2])), 
    Abs[\[Epsilon][t]] > 
     4}, {(1 + \[Eta][t])/(Pi \[CapitalDelta]p[t])
       Sign[\[Epsilon][
        t]] (\[Alpha][
          t] EllipticPi[-\[Beta][t] \[Eta][t], (Abs[\[Eta][t]])^2] - 
        EllipticPi[-\[Eta][t], (Abs[\[Eta][t]])^2]) + 
     Sign[\[Gamma][t]] G2[t] - ((I \[Eta][t])/(
       Pi \[CapitalDelta]p[
         t])) (\[Beta][
          t] EllipticPi[-\[Alpha][t] (1 + \[Eta][t]), (Abs[\[Eta]1[
             t]])^2] - 
        EllipticPi[(1 + \[Eta][t]), (Abs[\[Eta]1[t]])^2]), 
    Abs[\[Epsilon][t]] < 4}}]]
1/(2 \[CapitalDelta]p[t]) (1 - 
   I/Pi NIntegrate[(1/(x + \[CapitalOmega][
           t])) ((1 - (x - \[Epsilon][t]/2)^2)/(1 - x^2))^(1/
         2), {x, -1, 1}])

Plot[-2/Pi Abs[w + 9] Re[
   1/(2 \[CapitalDelta]p[w]) (1 - 
      I/Pi NIntegrate[(1/(x + \[CapitalOmega][
              w])) ((1 - (x - \[Epsilon][w]/2)^2)/(1 - x^2))^(1/
            2), {x, -1, 1}])], {w, -12.5, -3}, 
 PlotRange -> {{-3, -12}, {-2, 2}}]
Plot[-1/Pi (w + 9) Re[
   Piecewise[{{1/(
       2 \[CapitalDelta]p[w]) (1 + 
         2/Pi (1 + \[Xi][
             w]) (\[Alpha][
              w] EllipticPi[-\[Beta][w] \[Xi][w], (Abs[\[Xi][w]])^2] -
             EllipticPi[-\[Xi][w], (Abs[\[Xi][w]])^2])), 
      Abs[\[Epsilon][w]] > 
       4}, {(1 + \[Eta][w])/(Pi \[CapitalDelta]p[w])
         Sign[\[Epsilon][
          w]] (\[Alpha][
            w] EllipticPi[-\[Beta][w] \[Eta][w], (Abs[\[Eta][w]])^2] -
           EllipticPi[-\[Eta][w], (Abs[\[Eta][w]])^2]) + 
       Sign[\[Gamma][w]] G2[w] - ((I \[Eta][w])/(
         Pi \[CapitalDelta]p[
           w])) (\[Beta][
            w] EllipticPi[-\[Alpha][w] (1 + \[Eta][w]), (Abs[\[Eta]1[
               w]])^2] - 
          EllipticPi[(1 + \[Eta][w]), (Abs[\[Eta]1[w]])^2]), 
      Abs[\[Epsilon][w]] < 4}}]], {w, -12.5, -3}, 
 PlotRange -> {{-3, -12.5}, {-2, 2}}]
Plot[-2/Pi Abs[w + 9] Im[
   1/(2 \[CapitalDelta]p[w]) (1 - 
      I/Pi NIntegrate[(1/(x + \[CapitalOmega][
              w])) ((1 - (x - \[Epsilon][w]/2)^2)/(1 - x^2))^(1/
            2), {x, -1, 1}])], {w, -12.5, -3}, 
 PlotRange -> {{-3, -12}, {0, 3}}]
Plot[-2/Pi Abs[w + 9] Im[
   Piecewise[{{1/(
       2 \[CapitalDelta]p[w]) (1 + 
         2/Pi (1 + \[Xi][
             w]) (\[Alpha][
              w] EllipticPi[-\[Beta][w] \[Xi][w], (Abs[\[Xi][w]])^2] -
             EllipticPi[-\[Xi][w], (Abs[\[Xi][w]])^2])), 
      Abs[\[Epsilon][w]] > 
       4}, {(1 + \[Eta][w])/(Pi \[CapitalDelta]p[w])
         Sign[\[Epsilon][
          w]] (\[Alpha][
            w] EllipticPi[-\[Beta][w] \[Eta][w], (Abs[\[Eta][w]])^2] -
           EllipticPi[-\[Eta][w], (Abs[\[Eta][w]])^2]) + 
       Sign[\[Gamma][w]] G2[w] - ((I \[Eta][w])/(
         Pi \[CapitalDelta]p[
           w])) (\[Beta][
            w] EllipticPi[-\[Alpha][w] (1 + \[Eta][w]), (Abs[\[Eta]1[
               w]])^2] - 
          EllipticPi[(1 + \[Eta][w]), (Abs[\[Eta]1[w]])^2]), 
      Abs[\[Epsilon][w]] < 4}}]], {w, -12.5, -3}, 
 PlotRange -> {{-3, -12.5}, {0, 3}}]
6

Hi Kevin,

just a banal observation: Have you seen that the indefinite integral does give an answer?!

EllipticInt.png
EllipticInt.png652×382 21.9 KB

Regards – Henrik
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=EllipticInt.png&userId=32203