Casimir energy of disk-plane system

All contents of these notebooks are copyright &copy Noah Graham, 2016, all rights reserved.  These programs may be used and modified for noncommercial research purposes, provided that citation to

N. Graham, "Exact Electromagnetic Casimir Energy of a Disk Opposite a Plane," arXiv:1606.01090, Phys. Rev. A 94, 032509 (2016)

is included in all publications or other products in which the program or any programs derived from it were used.  This program and all accompanying materials are provided "AS IS," WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.

(* δ=0 electric δ=1 magnetic *)
<<Spheroidal6`
ParallelTable[<<Spheroidal6`,{i,1,$KernelCount}]

nmδToIndex[n_,m_,δ_]:=2((n^2)+n +m)+δ+1
indexTon[i_]:=Floor[Sqrt[((i-1)/2)]]
indexTom[i_]:=((i-nmδToIndex[indexTon[i],0,indexToδ[i]])/2)
indexToδ[i_]:=Mod[i-1,2]
minnSphere=1
minnSpheroid=0
maxn = 4
prec=2
minindexSphere=2(minnSphere^2)+1
minindexSpheroid=2(minnSpheroid^2)+1

maxindex=2(maxn^2)
ncoeff[n_,m_]:=(((2n+1)(n-m)!)/((n+m)!))
zero=(10^(-9))(-I)

integraloptions:=Sequence[MaxRecursion->16,AccuracyGoal->MachinePrecision-2,Method->DoubleExponential]
integraloptions4 :=Sequence[PrecisionGoal->2prec,AccuracyGoal->4prec,Method->GlobalAdaptive,MaxRecursion->24]

sumoptions:=Sequence[Method->AlternatingSigns,VerifyConvergence->False]
s3Ratio[n_,m_,zeroγ]:=s3Ratio[n,m,zeroγ]=((-I Derivative[0,0,0,1][LegendreQ][n,m,3,zero])/(LegendreQ[n,m,3,zero]))
s3Ratio[n_,m_,γ_]:=s3Ratio[n,m,γ]=((-I SpheroidalS3Prime[n,m,zero,I γ])/(SpheroidalS3[n,m,zero,I γ]))
s1Ratio[n_,m_,γ_]:=s1Ratio[n,m,γ]=((SpheroidalS1[n,m,zero,I γ])/(SpheroidalS3[n,m,zero,I γ]))
s1PrimeRatio[n_,m_,γ_]:=s1PrimeRatio[n,m,γ]=((SpheroidalS1Prime[n,m,zero,I γ])/(SpheroidalS3Prime[n,m,zero,I γ]))
a1ScaledZero[n_?NumberQ,m_]:=If[m<0,(((-1))^m)a1ScaledZero[n,-m](Sqrt[((ncoeff[n,m])/(ncoeff[n,-m]))]),If[OddQ[n-m],0,a1ScaledZero[n,m]=ncoeff[n,m]s3Ratio[n,m,zeroγ](((2^(m-1)) π(((-1))^((m-n)/2))Gamma[((n+1)/2)])/( Gamma[(1/2)-(m/2)-(n/2)] Gamma[(n/2)+1]((n-m)/2)!))]]
a1Scaled[n_?NumberQ,m_,γ_]:=If[m<0,(((-1))^m)a1Scaled[n,-m,γ](Sqrt[((ncoeff[n,m])/(ncoeff[n,-m]))]),If[OddQ[n-m],0,a1Scaled[n,m,γ]=ncoeff[n,m]s3Ratio[n,m,γ]
If [m==0&&γ≠0,γ SpheroidalS1[n,0,0,I γ]SpheroidalS2[n,0,0,I γ]SpheroidalRadialFactor[n,0,I γ],NIntegrate[SpheroidalPS[n,m,η,I γ]((Cos[γ(Sqrt[1-(η^2)])])/(Sqrt[1-(η^2)])),{η,0,1},Evaluate[integraloptions]]]]]
b1Scaled[n_?NumberQ,m_,γ_]:=If[m<0,(((-1))^m)b1Scaled[n,-m,γ](Sqrt[((ncoeff[n,m])/(ncoeff[n,-m]))]),If[OddQ[n-m],0,b1Scaled[n,m,γ]=ncoeff[n,m]s3Ratio[n,m,γ]
If[m==0,γ ((SpheroidalS1[n,0,0,I γ])^2)SpheroidalRadialFactor[n,0,I γ],NIntegrate[SpheroidalPS[n,m,η,I γ]((Sin[γ(Sqrt[1-(η^2)])])/(Sqrt[1-(η^2)])),{η,0,1},Evaluate[integraloptions]]]]]
b1ScaledLeadingOrder[n_?NumberQ,m_]:=If[m<0,(((-1))^m)b1ScaledLeadingOrder[n,-m](Sqrt[((ncoeff[n,m])/(ncoeff[n,-m]))]),If[OddQ[n-m],0,b1ScaledLeadingOrder[n,m]=ncoeff[n,m]s3Ratio[n,m,zeroγ]If[m==0,KroneckerDelta[n],(((((-1))^n)(2^(m-2))m Gamma[(n/2)]Gamma[((n+m+1)/2)])/(Gamma[((n+3)/2)]((n-m)/2)!))]]]

a2Scaled[n_?NumberQ,m_,γ_]:=If[m<0,(((-1))^(m+1))a2Scaled[n,-m,γ](Sqrt[((ncoeff[n,m])/(ncoeff[n,-m]))]),If[EvenQ[n-m] || m==0,0,a2Scaled[n,m,γ]=((m ncoeff[n,m])/(s3Ratio[n,m,γ]))NIntegrate[SpheroidalPS[n,m,(Sqrt[1-(u^2)]),I γ]((Sin[γ u])/u),{u,0,1},Evaluate[integraloptions]]]]
a2ScaledLeadingOrder[n_?NumberQ,m_]:=If[m<0,(((-1))^(m+1))a2ScaledLeadingOrder[n,-m](Sqrt[((ncoeff[n,m])/(ncoeff[n,-m]))]),If[EvenQ[n-m] || m==0,0,a2ScaledLeadingOrder[n,m]=((m ncoeff[n,m])/(s3Ratio[n,m,zeroγ]))(((((-1))^((-1+m-n)/2)) (2^(m-1)) π Gamma[(n/2)])/(Gamma[-(m/2)-(n/2)] Gamma[((3+n)/2)] Gamma[(1/2) (1-m+n)]))]]
b2Scaled[n_?NumberQ,m_,γ_]:=If[m<0,(((-1))^(m+1))b2Scaled[n,-m,γ](Sqrt[((ncoeff[n,m])/(ncoeff[n,-m]))]),If[EvenQ[n-m],0,If[m==0,-(((2n+1)SpheroidalPS[n,0,1,I γ])/(s3Ratio[n,0,γ])),b2Scaled[n,m,γ]=-((m ncoeff[n,m])/(s3Ratio[n,m,γ]))NIntegrate[SpheroidalPS[n,m,(Sqrt[1-(u^2)]),I γ]((Cos[γ u])/u),{u,0,1},Evaluate[integraloptions]]]]]
b2ScaledZero[n_?NumberQ,m_]:=If[m<0,(((-1))^(m+1))b2ScaledZero[n,-m](Sqrt[((ncoeff[n,m])/(ncoeff[n,-m]))]),If[EvenQ[n-m],0,b2ScaledZero[n,m]=-((ncoeff[n,m])/(s3Ratio[n,m,zeroγ]))(((((-1))^((-1+m-n)/2)) (2^m) π Gamma[((1+n)/2)])/( Gamma[-(m/2)-(n/2)] Gamma[1+(n/2)] Gamma[(1/2) (1-m+n)]))]]

sa1[m_,γ_]:=If[m<0,sa1[-m,γ],sa1[m,γ]=-(((m-1)!!)/((m-2)!!))If[EvenQ[m],1,(2/π)]+NSum[a1Scaled[n,m,γ]SpheroidalPS[n,m,0,I γ]-a1ScaledZero[n,m]LegendreP[n,m,0],{n,m,∞,2},Evaluate[sumoptions]]]
sb1[m_,γ_]:= If[m<0, sb1[-m,γ],sb1[m,γ]=-γ ((m!!)/((m-1)!!))If[OddQ[m],1,(2/π)]+(1/I)NSum[I(b1Scaled[n,m,γ]SpheroidalPS[n,m,0,I γ]-γ b1ScaledLeadingOrder[n,m]LegendreP[n,m,0]),{n,m,∞,2},Evaluate[sumoptions]]]
sa2[m_,γ_]:=If[m<0,-sa2[-m,γ],sa2[m,γ]= -γ (((m-1)!!)/((m-2)!!))If[EvenQ[m],1,(2/π)]+(1/I)NSum[I(a2Scaled[n,m,γ]SpheroidalPSPrime[n,m,0,I γ]-γ a2ScaledLeadingOrder[n,m]Derivative[0,0,1][LegendreP][n,m,0]),{n,m+1,∞,2},Evaluate[sumoptions]]]
sb2[m_,γ_]:=If[m<0,-sb2[-m,γ],sb2[m,γ]=((m!!)/((m-1)!!))If[OddQ[m],1,(2/π)]+NSum[b2Scaled[n,m,γ]SpheroidalPSPrime[n,m,0,I γ]-b2ScaledZero[n,m]Derivative[0,0,1][LegendreP][n,m,0],{n,m+1,∞,2},Evaluate[sumoptions]]]
q1[m_,γ_]:=((sa2[m,γ])/(sb2[m,γ]))
q2[m_,γ_]:=((sa1[m,γ])/(sb1[m,γ]))
α1Scaled[n_,m_,γ_]:=((SpheroidalPS[n,m,0,I γ]s1Ratio[n,m,γ])/(sa1[m,γ]-q1[m,γ]sb1[m,γ]))
β1Scaled[n_,γ_]:=((SpheroidalPS[n,0,0,I γ]s1Ratio[n,0,γ])/(sb1[0,γ]))

α2Scaled[n_,m_,γ_]:=((SpheroidalPSPrime[n,m,0,I γ]s1PrimeRatio[n,m,γ])/(sa2[m,γ]-q2[m,γ]sb2[m,γ] ))
β2Scaled[n_,γ_]:=((SpheroidalPSPrime[n,0,0,I γ]s1PrimeRatio[n,0,γ])/(sb2[0,γ]))

𝒯matrixElementEE[n_,np_,m_,γ_]:=If[OddQ[n-m]||OddQ[np-m],0,
If[n≠np,0,-s1Ratio[n,m,γ]]+If[m==0,β1Scaled[np,γ]b1Scaled[n,0,γ],α1Scaled[np,m,γ](a1Scaled[n,m,γ]-q1[m,γ]b1Scaled[n,m,γ])]((SpheroidalS3Prime[np,m,zero,I γ])/(SpheroidalS3Prime[n,m,zero,I γ]))]
𝒯matrixElementME[n_,np_,m_,γ_]:=If[EvenQ[n-m]||OddQ[np-m],0,If[m==0,0,α1Scaled[np,m,γ](a2Scaled[n,m,γ]-q1[m,γ]b2Scaled[n,m,γ])]((-I SpheroidalS3Prime[np,m,zero,I γ])/(SpheroidalS3[n,m,zero,I γ]))]
𝒯matrixElementEM[n_,np_,m_,γ_]:=If[OddQ[n-m]||EvenQ[np-m],0,If[m==0,0,α2Scaled[np,m,γ](a1Scaled[n,m,γ]-q2[m,γ]b1Scaled[n,m,γ])]((SpheroidalS3[np,m,zero,I γ])/(-I SpheroidalS3Prime[n,m,zero,I γ]))]
𝒯matrixElementMM[n_,np_,m_,γ_]:=If[EvenQ[n-m]||EvenQ[np-m],0,If[n≠np,0,-s1PrimeRatio[n,m,γ]]+If[m==0,β2Scaled[np,γ]b2Scaled[n,0,γ],α2Scaled[np,m,γ](a2Scaled[n,m,γ]-q2[m,γ]b2Scaled[n,m,γ])]((SpheroidalS3[np,m,zero,I γ])/(SpheroidalS3[n,m,zero,I γ]))]

𝒯matrixElement[n_,δ_,np_,δp_,m_,γ_]:=
(I^(np-n))(Sqrt[ ((ncoeff[np,m])/(ncoeff[n,m]))])If[δ==0,If[δp==0,𝒯matrixElementEE[n,np,m,γ],𝒯matrixElementEM[n,np,m,γ]],If[δp==0,𝒯matrixElementME[n,np,m,γ],𝒯matrixElementMM[n,np,m,γ]]]

fresnel[κ_,δplanar_]:=(((-1))^(δplanar)) 

translationMatrixElement[κ_,kperp_?NumberQ,δplanar_,n_,m_,δ_,θ_,ϕk_]:=
Module[{Cosθq=((I kperp Cos[ϕk])/κ),SinθqSinϕq=(kperp/(I κ)) Sin[ϕk],SinθqCosϕq=((Sqrt[(kperp^2)+(κ^2)])/κ)},
(Sqrt[ ((4 π ncoeff[n,m])/(n(n+1)))])((((( SinθqCosϕq+I SinθqSinϕq)/(SinθqCosϕq-I SinθqSinϕq))))^(m/2)) (κ/kperp) ((Cosθq)/(1-((Cosθq)^2)))If[δplanar==δ,
(Cosθq(n+1) SinθqCosϕq+I m ((SinθqSinϕq)/(Cosθq))) LegendreP[n,m,3,Cosθq]-(n+1-m) SinθqCosϕq LegendreP[1+n,m,3,Cosθq],
(((-1))^δ) (((n+1)SinθqSinϕq-I m SinθqCosϕq) LegendreP[n,m,3,Cosθq]-(n+1-m) ((SinθqSinϕq)/(Cosθq)) LegendreP[1+n,m,3,Cosθq])]]/;θ==(π/2)

translationMatrixElement[κ_,kperp_?NumberQ,δplanar_,n_,m_,δ_,θ_,ϕk_]:=
Module[{Cosθq=((I(Sqrt[(kperp^2)+(κ^2)])Cos[θ]-kperp Cos[ϕk]Sin[θ])/(I κ)),SinθqSinϕq=(kperp/(I κ)) Sin[ϕk],SinθqCosϕq= Cos[θ](kperp/(I κ))Cos[ϕk]+Sin[θ]((Sqrt[(kperp^2)+(κ^2)])/κ)},
(Sqrt[ ((4 π ncoeff[n,m])/(n(n+1)))])((((( SinθqCosϕq+I SinθqSinϕq)/(SinθqCosϕq-I SinθqSinϕq))))^(m/2)) (κ/kperp) If[δplanar==δ,
I m LegendreP[n,m,3,Cosθq]Sin[θ] ((SinθqSinϕq)/(1-((Cosθq)^2)))+(Sin[θ]SinθqCosϕq Cosθq-Cos[θ] (1-((Cosθq)^2))) Derivative[0,0,0,1][LegendreP][n,m,3,Cosθq],
(((-1))^δ) (I m LegendreP[n,m,3,Cosθq] (Cos[θ]-((SinθqCosϕq)/(1-((Cosθq)^2))) Cosθq Sin[θ])+Sin[θ] SinθqSinϕq Derivative[0,0,0,1][LegendreP][n,m,3,Cosθq])]]

𝒯matrix[γ_?NumberQ]:=ParallelTable[Module[{n=indexTon[i],m=indexTom[i],δ=indexToδ[i],np=indexTon[ip],mp=indexTom[ip],δp=indexToδ[ip]},
If[m≠mp,0,𝒯matrixElement[n,δ,np,δp,m,γ]]],{i,minindexSpheroid,maxindex},{ip,minindexSpheroid,maxindex}]

convertMatrixElement[n_,np_,m_,δ_,κ_,radius_,pow_]:=If[Mod[n-np,2]≠0,0,Module[{maxnround=maxn+Mod[maxn,2]},
((((Sqrt[n(n+1)])))^pow)(Sqrt[((ncoeff[np,m])/(ncoeff[n,m]))])(((-1))^((n-np)/2))Spheroidal`Private`angcoeff[np,m,SetPrecision[(((κ radius))^2),24],24,(maxnround/2)][[((n-np+maxnround)/2)+1]]]]

convertMatrix[κ_?NumberQ,radius_,pow_]:=Table[Module[{n=indexTon[i],m=indexTom[i],δ=indexToδ[i],np=indexTon[ip],mp=indexTom[ip],δp=indexToδ[ip]},
If[m≠mp || δ≠δp,0,convertMatrixElement[n,np,m,δ,κ,radius,pow]]],{i,minindexSphere,maxindex},{ip,minindexSpheroid,maxindex}]


translationMatrixElement[κ_,kperp_?NumberQ,δplanar_,n_,m_,δ_]:=(Sqrt[ ((4 π ncoeff[n,m])/(n(n+1)))])If[δplanar==δ,(kperp/κ)Derivative[0,0,1][LegendreP][n,m,((Sqrt[(κ^2)+(kperp^2)])/κ)],(((-1))^δ)(κ/kperp)I m LegendreP[n,m,((Sqrt[(κ^2)+(kperp^2)])/κ)]]
𝒰matrixElement[κ_?NumberQ,separation_,θ_,ℓ_,m_,δ_,ℓp_,mp_,δp_]:=
If[m≠mp,0,(1/(4 π κ))NIntegrate[kperp((Exp[-2separation (Sqrt[(κ^2)+(kperp^2)])])/(Sqrt[(κ^2)+(kperp^2)]))Sum[translationMatrixElement[κ,kperp,δplanar,ℓ,m,δ]fresnel[κ,δplanar]
(((-1))^m)translationMatrixElement[κ,kperp,δplanar,ℓp,m,δp],{δplanar,0,1}],{kperp,0,∞},Evaluate[integraloptions2]]]/;θ==0

𝒰matrixElementIntegrand[κ_?NumberQ,separation_,θ_,ℓ_,m_,δ_,ℓp_,mp_,δp_,kperp_,ϕk_]:=(kperp/(Sqrt[(kperp^2)+(κ^2)]))Exp[-2separation (Sqrt[(κ^2)+(kperp^2)])]Sum[translationMatrixElement[κ,kperp,δplanar,ℓ,m,δ,θ,ϕk]fresnel[κ,δplanar]
(((-1))^(δp-δplanar))translationMatrixElement[κ,kperp,δplanar,ℓp,-mp,δp,-θ,ϕk],{δplanar,0,1}]

𝒰matrixElement[κ_?NumberQ,separation_,θ_,ℓ_,m_,δ_,ℓp_,mp_,δp_]:=0/;θ==(π/2)&&(OddQ[ℓ-ℓp+δ-δp]||(ℓ==Abs[m]&&ℓp==Abs[mp]&&fresnel[κ,δ]==-fresnel[κ,δp]&&Sign[m]==Sign[mp]))
𝒰matrixElement[κ_?NumberQ,separation_,θ_,ℓ_,m_,δ_,ℓp_,mp_,δp_]:=
(𝒰matrixElement[κ,separation,θ,ℓ,m,δ,ℓp,mp,δp]=
(((-1))^(δp-δ))(𝒰matrixElement[κ,separation,θ,ℓp,-mp,δp,ℓ,-m,δ]=(((-1))^(mp-m))(𝒰matrixElement[κ,separation,θ,ℓ,-m,δ,ℓp,-mp,δp]=(((-1))^(δ-δp))(𝒰matrixElement[κ,separation,θ,ℓp,mp,δp,ℓ,m,δ]=(((-1))^(m-mp))(1/(4 (π^2) κ))NIntegrate[
𝒰matrixElementIntegrand[κ,separation,θ,ℓ,m,δ,ℓp,mp,δp,kperp,ϕk]+𝒰matrixElementIntegrand[κ,separation,θ,ℓ,m,δ,ℓp,mp,δp,kperp,-ϕk],
{kperp,0,∞},{ϕk,0,(π/2)},Evaluate[integraloptions4]]))))/;θ==(π/2)

𝒰matrixElement[κ_?NumberQ,separation_,θ_,ℓ_,m_,δ_,ℓp_,mp_,δp_]:=
(𝒰matrixElement[κ,separation,θ,ℓ,m,δ,ℓp,mp,δp]=
(((-1))^(δp-δ))(𝒰matrixElement[κ,separation,θ,ℓp,-mp,δp,ℓ,-m,δ]=(((-1))^(mp-m))(𝒰matrixElement[κ,separation,θ,ℓ,-m,δ,ℓp,-mp,δp]=(((-1))^(δ-δp))(𝒰matrixElement[κ,separation,θ,ℓp,mp,δp,ℓ,m,δ]=(((-1))^(m-mp))(1/(8 (π^2) κ))NIntegrate[
𝒰matrixElementIntegrand[κ,separation,θ,ℓ,m,δ,ℓp,mp,δp,kperp,ϕk]+𝒰matrixElementIntegrand[κ,separation,θ,ℓ,m,δ,ℓp,mp,δp,kperp,π-ϕk]
+𝒰matrixElementIntegrand[κ,separation,θ,ℓ,m,δ,ℓp,mp,δp,kperp,-ϕk]+𝒰matrixElementIntegrand[κ,separation,θ,ℓ,m,δ,ℓp,mp,δp,kperp,π+ϕk],
{kperp,0,∞},{ϕk,0,(π/2)},Evaluate[integraloptions4]]))))

𝒰matrix[κ_?NumberQ,separation_,θ_]:=ParallelTable[Module[{ℓ=indexTon[i],m=indexTom[i],δ=indexToδ[i],ℓp=indexTon[ip],mp=indexTom[ip],δp=indexToδ[ip]},
If[m==-mp && fresnel[κ,δ]==-fresnel[κ,δp]&&ℓ==ℓp,0, 𝒰matrixElement[κ,separation,θ,ℓp,mp,δp,ℓ,m,δ]]],{i,minindexSphere,maxindex},{ip,minindexSphere,maxindex}]

𝒯matrixSpherical[κ_?NumberQ,radius_]:=convertMatrix[κ,radius,1].𝒯matrix[I κ radius].Transpose[convertMatrix[κ,radius,-1]]

𝒩matrix[κ_?NumberQ,separation_,radius_,θ_]:=𝒯matrixSpherical[κ,radius].𝒰matrix[κ,separation,θ]

printAndReturn[κ_,x_]:=(Print[κ," ",x," ",SessionTime[]];x)
casimirIntegrand[κ_?NumberQ,separation_,radius_,θ_]:=printAndReturn[κ,Log[Det[IdentityMatrix[maxindex-minindexSphere+1]-𝒩matrix[κ,separation, radius,θ]]]]

max=4
min=(1/128)
casimirEnergy[separation_,radius_,θ_]:=(1/(2 π))NIntegrate[casimirIntegrand[κ,separation, radius,θ],{κ,(min/(Max[radius,separation])),(max/(Min[radius,separation]))},PrecisionGoal->prec]

DistributeDefinitions["Global`"]
Module[{separation=(3/2),radius=1,θ=(π/2)},casimirEnergy[separation,radius,θ]]

NotebookSave[]