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

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?

• Why do you hardcode the cutoff at 100 in SumAux instead of using the nmax parameter? – Roman Jan 26 '19 at 7:58
• @Roman, I did it when testing, later I've removed it and used the variable instead. I've updated the code in the post. – user1620696 Jan 26 '19 at 13:06

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:

• The Piecewise is needed because sometimes numerical inaccuracies give slightly negative eigenvalues. Take, for example, a vector v = RandomReal[{-1, 1}, 10] and its density matrix R = KroneckerProduct[v, v]. The eigenvalues of R should be all zero except one of them. Numerically, though, they fluctuate around zero: Eigenvalues[R]. This creates trouble for the s function if defined without Piecewise. – Roman Jan 26 '19 at 7:55
• I think the paper uses the natural logarithm, not the two's logarithm, in the definition of the entropy: use s[x_] = Piecewise[{{x*Log[x], 0 < x < 1}}] instead of what you've used. – Roman Jan 26 '19 at 8:19
• There indeed was one sign issue. The last term in the equation has a minus in front. Only when the minus is included in the sum we get the sum of two entropies. I left the minus sign together with the entropies in the code. I've corrected these issues. The function s now doesn't include the sign and it is placed where it is required. The graph got better, but it is still wrong. I've tried setting a higher cutoff (300) but it took forever to build the graph and didn't change almost. – user1620696 Jan 26 '19 at 13:29
• I see @Roman! Indeed these expressions numerically are rather problematic. Mathematica complains about a precision loss even for nmax as low as 40. @user1620696, the notch you notice is due to convergence issues. You will notice that as you move the nmax higher, the notch moves towards the left. Trying nmax=300 gives the output I have uploaded, which I think is pretty similar to the one in the paper. – Sotiris Jan 26 '19 at 15:23
• @Sotiris thanks, I've found what was wrong all along (apart from the sign mistake): the plot I got had the $y$ axis starting at $0.7$ and I thought it was starting at $0.0$ giving me the impression the graph was totally wrong. I've made the axis start at $0.0$ with PlotRange and got the same result as you, which agrees with the paper. Thanks a lot ! – user1620696 Jan 28 '19 at 0:02