# Numerically computing the eigenvalues of an infinite-dimensional tridiagonal matrix

I have one infinite dimensional tridiagonal matrix whose eigenvalues I have to compute. How can that be done numerically using Mathematica?

Let me expose the concrete case I want to do it. I shall use Dirac's notation to write down the matrix. It is:

$$\rho= \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)$$

Where $$M_{00},M_{01},M_{10},M_{11}$$ are the infinite dimensional matrices:

$$M_{00}=\sum_{n=0}^\infty t^{2n}|n\rangle\langle n|,\quad M_{11}=(1-t^2)\sum_{n=0}^\infty (n+1)t^{2n}|n+1\rangle\langle n+1|$$

$$M_{01}=\sqrt{1-t^2}\sum_{n=0}^\infty\sqrt{n+1} t^{2n}|n\rangle\langle n+1|,\quad M_{10}=M_{01}^\dagger$$

The notation is such that $$M=\sum_{nm}M_{nm}|n\rangle\langle m|$$

where $$n$$ labels the row and $$m$$ the colum. So $$M_{00}$$ and $$M_{11}$$ are diagonal and the other two parts gives the off-diagonal terms just above and below the diagonal.

According to this paper:

The density matrices $$\rho$$ are tridiagonal, whose eigenvalues can be obtained easily numerically.

So they claim it is easy to find numerically the eigenvalues of such matrix.

I just have no idea how to do it because first of all these are infinite dimensional matrices which I don't know how to define on Mathematica.

So can I use Mathematica to compute the eigenvalues of such a matrix? If so, how can I do it?

I want in the end to have the eigenvalues defined as functions of $$\theta$$.

• What are these matrix components in Mathematica code? (I'm not asking for infinite dimensional vectors, just a clear indication of how one might form a finite upper left submatrix). Jan 24, 2019 at 16:22
• Have you tried summing only to a finite $n_{\text{max}}$ and seeing how the numerical eigenvalues converge as this upper limit increases? Jan 24, 2019 at 16:44
• Not familiar with the notation...are your matrices Toeplitz? Or does Roman's code generate the correct matrices? Jan 24, 2019 at 18:34
• The matrices are the ones generated by @Roman's code. Jan 25, 2019 at 14:18
• Maybe I'll comment here that such matrices are specific to quantum mechanics and have properties that guarantee that in a low-energy subspace (i.e., a finite cutoff $n\le n_{\text{max}}$) the results are fairly accurate; there is no need for really taking $n$ to infinity. This can be motivated with arguments from physics. In a purely mathematical context, however, my recommendations and the code below would be insufficient in general. Jan 25, 2019 at 15:12

For a finite upper cutoff $$n\le n_{\text{max}}$$ you can define the matrices with

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]]


(assuming here that $$t\in\mathbb{R}$$ so that the Hermitian transpose is just the transpose) and the density matrices with

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]))


Then you get the eigenvalues numerically. For example, $$n_{\text{max}}=10$$, $$t=0.2$$, $$\theta=0.5$$:

Eigenvalues[rhoplus[10, 0.2, 0.5]]
(* list of eigenvalues *)


Then you can study the convergence of these eigenvalues as $$n_{\text{max}}$$ becomes large.

To find out whether your $$n_{\text{max}}$$ is large enough, you can look at the trace of the density matrix: it approaches 1 as $$n_{\text{max}}\to\infty$$. If it is much smaller than 1, then you should probably increase $$n_{\text{max}}$$.

Tr[rhoplus[3, 0.2, 0.5]]
(* 0.9999823973452051 *)
Tr[rhoplus[10, 0.2, 0.5]]
(* 0.9999999999999929 *)
Tr[rhoplus[10, 0.9, 0.5]]
(* 0.8859700675762211 *)
Tr[rhoplus[100, 0.9, 0.5]]
(* 0.9999999985998526 *)

• Thanks Roman ! So for instance, if I wanted to find the entropy of the density matrix, I would have to use the function Entropy on the list returned by Eigenvalues, which would depend on $n$, and then take the limit as $n\to \infty$ using DiscreteLimit? Jan 25, 2019 at 14:21
• I think the function Entropy does something different. For the von Neumann entropy I'd define a helper function s[x_] = Piecewise[{{-x*Log[2,x], 0<x<1}}] (or any other logarithm) and then calculate Total[s /@ Eigenvalues[rhoplus[10, 0.2, 0.5]]]. The role of the function s is to make sure that spurious negative eigenvalues (due to numerical inaccuracies) do not create trouble. Jan 25, 2019 at 14:46
• Convergence looks exponential in $n_{\text{max}}$ so just keep increasing $n_{\text{max}}$ until the result no longer changes. I don't think DiscreteLimit or NumericalMathNSequenceLimit` are necessary here. Jan 25, 2019 at 14:53