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

  • 3
    $\begingroup$ 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). $\endgroup$ Jan 24, 2019 at 16:22
  • $\begingroup$ Have you tried summing only to a finite $n_{\text{max}}$ and seeing how the numerical eigenvalues converge as this upper limit increases? $\endgroup$
    – Roman
    Jan 24, 2019 at 16:44
  • $\begingroup$ Not familiar with the notation...are your matrices Toeplitz? Or does Roman's code generate the correct matrices? $\endgroup$
    – MikeY
    Jan 24, 2019 at 18:34
  • $\begingroup$ The matrices are the ones generated by @Roman's code. $\endgroup$ Jan 25, 2019 at 14:18
  • 1
    $\begingroup$ 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. $\endgroup$
    – Roman
    Jan 25, 2019 at 15:12

1 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` *)
  • $\begingroup$ 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? $\endgroup$ Jan 25, 2019 at 14:21
  • $\begingroup$ 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. $\endgroup$
    – Roman
    Jan 25, 2019 at 14:46
  • $\begingroup$ 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 NumericalMath`NSequenceLimit are necessary here. $\endgroup$
    – Roman
    Jan 25, 2019 at 14:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.