I have two integrals that I am calculating using nested NIntegrate. One is for the singlet helium atom and another for triplet. (here, I am calculating variationally the approximate energies and eigenstates of 1s2p state of neutral helium).

First Integral (Singlet):

f1Final[R2_?NumericQ,Î¸2_?NumericQ,Î±_?NumericQ,Î²_?NumericQ]:=2*Pi*NIntegrate[

((Î±^3 Î²^5 (E^(-R2 Î±-R1 Î²) R1 Cos[Î¸1]+E^(-R1 Î±-R2 Î²) R2 Cos[Î¸2])^2) R1^2 Sin[Î¸1])/(2 Ï€^2 Sqrt[R1^2+R2^2-2 R1 R2 Cos[Î¸2]]),

{R1,0,âˆž},{Î¸1,0,Ï€},

Method->{“GlobalAdaptive”,”MaxErrorIncreases”->10000,Method->”GaussKronrodRule”},AccuracyGoal->5,PrecisionGoal->5]

fFinal[{Î±_?NumericQ,Î²_?NumericQ}]:=2*Pi*NIntegrate[

f1Final[R2,Î¸2,Î±,Î²] R2^2 Sin[Î¸2],{R2,0,âˆž},{Î¸2,0,Ï€},

Method->{“GlobalAdaptive”,”MaxErrorIncreases”->10000,Method->”GaussKronrodRule”},AccuracyGoal->5,PrecisionGoal->5]

Second Integral (Triplet):

f1FinalMin[R2_?NumericQ,Î¸2_?NumericQ,Î±_?NumericQ,Î²_?NumericQ]:=2*Pi*NIntegrate[

((Î±^3 Î²^5 (-E^(-R2 Î±-R1 Î²) R1 Cos[Î¸1]+E^(-R1 Î±-R2 Î²) R2 Cos[Î¸2])^2) R1^2 Sin[Î¸1])/(2 Ï€^2 Sqrt[R1^2+R2^2-2 R1 R2 Cos[Î¸2]]),

{R1,0,âˆž},{Î¸1,0,Ï€},

Method->{“GlobalAdaptive”,”MaxErrorIncreases”->10000,Method->”GaussKronrodRule”},AccuracyGoal->5,PrecisionGoal->5]

fFinalMin[{Î±_?NumericQ,Î²_?NumericQ}]:=2*Pi*NIntegrate[

f1FinalMin[R2,Î¸2,Î±,Î²]*R2^2*Sin[Î¸2],

{R2,0,âˆž},{Î¸2,0,Ï€},

Method->{“GlobalAdaptive”,”MaxErrorIncreases”->10000,Method->”GaussKronrodRule”},AccuracyGoal->5,PrecisionGoal->5]

The first integral is for the singlet state. The second for triplet state. The wave function goes to zero, for the singlet state, when R1=R2R1=R2 and θ1\theta1 = θ2\theta2. For some reason, the two integrals give the same values for any input of α\alpha and β\beta. Does anyone know why this might be the case? The two integrals should differ, because the triplet state (and therefore the integrand above for f1FinalMin) goes to zero for R1=R2R1=R2 and θ1\theta1 = θ2\theta2. NIntegrate doesn’t seem to be working properly.

Also, if anyone has an idea for a better way to perform these integrals, that would be great.

Thanks!

=================

=================

1 Answer

1

=================

For problems like these, I like to take small steps and check each one, so please bear with me. First, we write the two exponential functions, which I am calling f1 and f2. Then write the next two more complicated functions, called gSing for singlet and gTrip for triplet. Then write an expression for the denominator and an expression for the volume elements.

f1 = Exp[-Î±*r1 – Î²*r2];

f2 = Exp[-Î±*r2 – Î²*r1];

gSing = f2*r1*Cos[Î¸1] + f1*r2*Cos[Î¸2];

gTrip = -f2*r1*Cos[Î¸1] + f1*r2*Cos[Î¸2];

d = Sqrt[r1^2 + r2^2 – 2*r1*r2*Cos[Î¸2]];

vol = r1^2*Sin[Î¸1]*r2^2*Sin[Î¸2];

The above expressions do not have some of your constant factors in them, but that is okay for now. Next we write expressions for the singlet and triplet integrands, and for the difference.

hSing = (gSing^2/d)*vol;

hTrip = (gTrip^2/d)*vol;

hDiff = FullSimplify[hSing – hTrip]

The result looks like this r31r32sin(2θ1)sin(2θ2)e(α+β)(−(r1+r2))√r21+r22−2r1r2cos(θ2)

\frac{r_1^3 r_2^3 \sin(2\theta_1) \sin(2\theta_2) e^{(\alpha +\beta)(-(r_1+r_2))}}{\sqrt{r_1^2+r_2^2-2r_1 r_2 \cos (\theta_2)}}

The difference in the two integrands contains a factor of sin(2θ1)\sin(2\theta_1) which integrates to zero. Therefore, the difference in the two integrals must be zero, even though the integrands are different.

If we go back and look at where that sin2θ1sin2\theta_1 term is coming from, it looks right. So, I’m going to suggest maybe the formula is wrong. Since you mentioned the 1S2P state, I would expect one electron to have the ground state (1S) wavefunction and the other to have the excited state (2P) wavefunction. So, that might be something to check / verify. Whatever the problem is, as long as that sin2θ1sin2\theta_1 term is there, Mathematica is right to give the same answer for the two integrals.

Another thing that might help is to look at Problem 5.11 in Griffiths. He gives a hint about a similar integral for the (1S)2(1S)^2 state. He says we should integrate w.r.t. d3r2d^3r_2 first. He also says to integrate r2r_2 from zero to r1r_1 and then from r1r_1 to infinity. Introduction to Quantum Mechanics by David J. Griffiths, second edition.

Since we have taken electron 1 to be on the polar axis with θ1=0\theta_1=0, maybe we should not be integrating w.r.t. θ1\theta_1 at all. If we set θ1=0\theta_1=0 to begin with, and remove the sin(θ1)\sin(\theta_1) term from the volume element expression, then we get a new expression for the difference in integrands. Using the hint in Griffiths, here is how this new expression could be integrated symbolically.

vol2 = r1^2*r2^2*Sin[Î¸2];

hSing2 = (gSing^2/d)*vol2 /. Î¸1 -> 0;

hTrip2 = (gTrip^2/d)*vol2 /. Î¸1 -> 0;

hDiff2 = FullSimplify[hSing2 – hTrip2]

Now to integrate hDiff2, we want to integrate w.r.t. θ2\theta_2 first, per Griffiths. We divide the integral into an outer integral with r2>r1r_2>r_1 and and inner integral with r2

inner1 = Integrate[hDiff2, {Î¸2, 0, Pi},

Assumptions -> {Î± > 0, Î²] > 0, r2 > 0, r2 < r1}]
Next, integrate w.r.t. r2r_2 from zero to r1r_1 and from r1r_1 to infinity. And finally, combined the results and integrate w.r.t. r1r_1 from zero to infinity. outer2 = Integrate[outer1, {r2, r1, Infinity},
Assumptions -> {Î± > 0, Î² > 0, r1 > 0}]

inner2 = Integrate[inner1, {r2, 0, r1},

Assumptions -> {Î± > 0, Î² > 0, r1 > 0}]

diffInt = Integrate[inner2 + outer2, {r1, 0, Infinity},

Assumptions -> {Î± > 0, Î² > 0}]

So, that is the basic concept. The justification is that if θ1\theta_1 was allowed to vary, then the distance between electrons would depend on θ1\theta_1, but it does not. Don’t forget all the constant factors I have left out, including a factor of 2 on that second volume element.

(+1) for taking the effort – this confirms that it’s not a Mathematica issue.

– Jens

Jan 25 at 21:58

did this work for you @Louis? I agree with you about θ1\theta_1 (although Griffiths doesn’t do what you did, which is confusing me more) but if I do the integral as you suggest i do not get the right result for singlet and triplet states; i took care of the constants too. the singlet state fshould give (assuming that the integral result is singInt) singInt+ 1/2 ($\alpha$^2 + $\beta$^2) – 2 $\alpha$ – $\beta$ = -2.12161. what is bothering me is that i get the right result for the singlet state but the triplet state is the same which makes sense given your analysis but shouldn’t be the case.

– Pawan Dhakal

Jan 26 at 5:39

@PawanDhakal. Yes, it worked for me. That is, for the Coulomb interaction energy of the electrons in the 1S2P state, I get 0.548 for the singlet and 0.445 for the triplet in units of q2Z/a0q^2 Z/a_0. I get positive values for repulsion, which seems to agree with the energy level diagram for He. Are we calculating the same thing? What is the right answer? Also, what textbook(s) are you looking at? I can show my full calculation, if you want. I will check it a few more times.

– LouisB

Jan 26 at 22:11

@Louis Yes. the repulsion should give positive values. The total energies for singlet state should be -57.76 eV (about -2.12 atomic units) and triplet should be lower by 0.22 eV. If you have the chance to look at P.J.E Peebles (Quantum Mechanics, Chapter 6), he does the calculation for 1s2p state and derives the results. But, when I do what you suggest, the energy values (after minimization) are -67.74eV for singlet and -67.79eV for triplet. The book I was mentioning was the Griffiths one, which doesn’t do θ1=0\theta_1 = 0 during the calculation of the ground state energy.

– Pawan Dhakal

Jan 27 at 4:38