Wolfram Cloud Document

In[]:=

Ntot=100;a=1;Λ=10^-10;rVals=Range[1,Ntot-1];tolerance=10^-5;ℓmaxLimit=100;constructK3DDirichlet[N_,Λ_,ℓ_]:=Module[{K,j},K=ConstantArray[0.,{N,N}];Do[j=i;K[[i,i]]=Λ+ℓ(ℓ+1)/j^2+(j+1/2)^2/j^2+(j-1/2)^2/j^2;If[i==1,K[[i,i]]=Λ+ℓ(ℓ+1)+9/4];If[i<N,K[[i,i+1]]=-(j+1/2)^2/(j(j+1));K[[i+1,i]]=K[[i,i+1]];],{i,1,N}];K/a^2];computeEntropy[n_,Ω_]:=Module[{ΩA,ΩB,ΩC,solveA,γ,β,Ξraw,Ξ,ξvals,id},Quiet@Module[{},ΩA=Ω[[1;;n,1;;n]];ΩB=Ω[[1;;n,n+1;;]];ΩC=Ω[[n+1;;,n+1;;]];solveA=Check[LinearSolve[ΩA],Return[0]];γ=ΩC-1/2Transpose[ΩB].solveA[ΩB];β=1/2Transpose[ΩB].solveA[ΩB];Ξraw=Check[LinearSolve[γ,β],Return[0]];id=IdentityMatrix[Length[Ξraw]];Ξ=Check[Ξraw.MatrixPower[id+MatrixPower[id-Ξraw.Ξraw,1/2],-1],Return[0]];ξvals=Eigenvalues[N[Ξ]];If[VectorQ[ξvals,NumericQ],Total[-Log[1-ξvals]-(ξvals/(1-ξvals))Log[ξvals]],0]]];Stotal=Table[Module[{ℓ=0,Sold=0,Snew=0,converged=False},While[ℓ<=ℓmaxLimit&&!converged,Module[{K,Ω,Sℓ},K=constructK3DDirichlet[Ntot,Λ,ℓ];Ω=MatrixPower[N[K],1/2];Sℓ=computeEntropy[r,Ω];Snew=Sold+(2ℓ+1)*Sℓ;If[ℓ>0&&Abs[Snew-Sold]/Snew<tolerance,converged=True,Sold=Snew;ℓ++];]];Snew],{r,rVals}];entropyData=Prepend[Transpose[{rVals^2,Stotal}],{0,0}];ListLinePlot[entropyData,Frame->True,FrameLabel->{"r² / a²","Entanglement Entropy "},PlotRange->{{0,3500},{0,1200}},PlotLabel->"\!\(\*SuperscriptBox[\(S\), \(total\)]\) vs. r²/a² (Dirichlet BC, convergence to 10⁻⁵)",PlotStyle->Thick,ImageSize->Large]