I'm trying to reproduce the computations of this paper and I'm running into some troubles because I'm rather new to Mathematica.
In the paper's page 4, in figure 2 the authors show a plot of a function. The function is as follows.
First define
$$\mathcal{D}(t,\theta)=1+\frac{1-t^2}{2}\sum_{n=0}^\infty t^{2n}(1+(n+1)(1-t^2))\log\left[\frac{1-t^2}{2}t^{2n}(1+(n+1)(1-t^2))\right]-\frac{1}{2}\sum_{i=\pm}\operatorname{Tr}(\rho_{R|i}\log \rho_{R|i})$$
The plot I'm trying to construct is $\mathcal{D}(\tanh r,\pi/2)$. In the above we have:
$$\rho_{R|\pm}=\dfrac{1-t^2}{2}\left(1\pm\cos\theta)M_{00}+(1\mp\cos\theta)M_{11}\pm\sin\theta M_{10}\pm\sin\theta M_{01}\right)$$
The matrices $M_{ij}$ above as well as the matrices $\rho_{R|i}$ were discussed in this question and the Mathematica code that generates them is given by @Roman's answer:
M00[nmax_Integer,t_] := SparseArray[Band[{1,1}] -> Table[t^(2n), {n,0,nmax}]]
M11[nmax_Integer,t_] := (1-t^2)*SparseArray[Band[{1,1}] -> Table[n*t^(2(n-1)), {n,0,nmax}]]
M01[nmax_Integer,t_] := Sqrt[1-t^2]*SparseArray[Band[{1,2}] -> Table[Sqrt[n+1]*t^(2n), {n,0,nmax-1}], {nmax+1,nmax+1}]
M10[nmax_Integer,t_] := Transpose[M01[nmax,t]]
rhoplus[nmax_Integer,t_,th_] := (1-t^2)/2*((1+Cos[th])*M00[nmax,t]+(1-Cos[th])*M11[nmax,t]+Sin[th]*(M10[nmax,t]+M01[nmax,t]))
rhominus[nmax_Integer,t_,th_]:= (1-t^2)/2*((1-Cos[th])*M00[nmax,t]+(1+Cos[th])*M11[nmax,t]-Sin[th]*(M10[nmax,t]+M01[nmax,t]))
With this we can compute the traces from the eigenvalues of these matrices. As also discussed in the question, we do:
s[x_] = Piecewise[{{x*Log[2,x], 0<x<1}}]
EntropyPlus[nmax_Integer, t_, th_] := Total[s /@ Eigenvalues[rhoplus[nmax,t,th]]]
EntropyMinus[nmax_Integer, t_, th_] := Total[s /@ Eigenvalues[rhominus[nmax,t,th]]]
I then computed the other sum with a cutoff of $N = 100$:
SumAuxElement[t_][n_Integer] := ((1 - t^2)/2)*t^(2 n)*(1 + (n + 1) (1 - t^2))
SumAuxList[nmax_Integer, t_] := Array[SumAuxElement[t], nmax + 1, {0, nmax}]
SumAux[t_,nmax_Integer]:=Total[s /@ SumAuxList[nmax, t]]
First I defined one auxilairy $n$-th element of the sum, assembled those in a list and called total on it using the $s$ function as with the eigenvalues. The minus sign is because on the function $s$ we have one minus that doesn't appear on the second term of $\mathcal{D}$.
I finally defined the function to be plotted:
Discord[nmax_Integer, t_, th_] := 1 + SumAux[t,nmax] - 1/2 (EntropyPlus[nmax, t, th] + EntropyMinus[nmax, t, th])
Plot[Discord[100, Tanh[r], Pi/2], {r, 0, 2.5}]
And the graph I get is:
Which is not the correct one (see below the green line):
So what am I doing wrong? I'm confident on the answer to the other question so that the last term in $\mathcal{D}$ is probably right. I believe my issue is with that infinite sum.
Why am I not getting the right plot here? What should I modify in this code?
SumAux
instead of using thenmax
parameter? $\endgroup$