I am working with an infinite-dimensional density matrix in Mathematica, where the basis states are:

∣000n⟩, ∣100n⟩, ∣010n⟩, ∣001n⟩, ∣110n⟩, ∣101n⟩, ∣000(n+1)⟩, ∣100(n+1)⟩

The density matrix is given by:

X = Sum[ (tanh(r)^(2n) / cosh(r)^2) [
    |100n⟩ ⟨100n| + 
    |010n⟩ ⟨010n| + 
    |001n⟩ ⟨001n| + 
    |001n⟩ ⟨010n| + 
    |010n⟩ ⟨001n| + 
    |110n⟩ ⟨000n| + 
    |101n⟩ ⟨000n| + 
    |000n⟩ ⟨110n| + 
    |000n⟩ ⟨101n| + 
    (sqrt(n+1) / cosh(r)) [
        |000n⟩ ⟨100(n+1)| + 
        |010n⟩ ⟨000(n+1)| + 
        |001n⟩ ⟨000(n+1)| + 
        |100(n+1)⟩ ⟨000n| + 
        |000(n+1)⟩ ⟨010n| + 
        |000(n+1)⟩ ⟨001n|
    ] + 
    (n+1 / cosh(r)^2) |000(n+1)⟩ ⟨000(n+1)|
]
,{n,0,Infinity}]

where n runs from 0 to infinity, meaning the full matrix is infinite-dimensional.

My Questions:

1- How can I define this density matrix in Mathematica?

2- Is there a more efficient way to represent it symbolically?

3-What is the best approach to compute its eigenvalues numerically?