What am I doing wrong when trying to plot this function?

Question

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?

Answer 1

It seems there is a problem with the signs in your expressions. Try using this (although the Piecewise is not needed)

s[x_] = Piecewise[{{x*Log[2,x], 0<x<1}}];

and then

SumAux[t_, nmax_Integer] := Total[s /@ SumAuxList[100, t]]

These should give you the paper's graph. Note you should type "SumAux[t, nmax]" instead of just "SumAux[t]", but since you got the first graph I suppose you already have! :)

Discord[nmax_Integer, t_, th_] := 1 + SumAux[t,nmax] - 1/2 (EntropyPlus[nmax, t, th] + EntropyMinus[nmax, t, th])

Trying for $nmax = 300$, you get the following output: