Numerically computing the eigenvalues of an infinite-dimensional tridiagonal matrix

Question

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$.

Answer 1

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` *)