For the number operator $\hat{N}$, it's eigenvectors are the Fock basis vectors $|n\rangle$, as $\hat{N}|n\rangle = n|n\rangle$. Let us suppose we have a bipartite set of basis vectors $\{|n,n\rangle\}$. I wish to pass a state into a parametric amplifier and measure the number operator of one of it's output modes, $$\hat{b} = \sqrt{G} \hat{a}_1 + \sqrt{G-1} \hat{a}_2^{\dagger}~,$$ where $G$ is the gain of the amplifier. We can measure the number of photons at the output by computing: $$\hat{b}^{\dagger}\hat{b} = G \hat{a}_1^{\dagger}\hat{a}_1 + (G-1)\hat{a}_2^{\dagger}\hat{a}_2 + G-1\\ +\sqrt{G(G-1)}\left[ \hat{a}_1^{\dagger}\hat{a}_2^{\dagger} + \hat{a}_2\hat{a}_1 \right].$$
My question is: is there a procedure to determine the eigenvectors of this new number operator? As this operator no longer appears to have the Fock basis as it's eigenvectors, since $\hat{b}^{\dagger}\hat{b}|n,n\rangle \ne \alpha |n,n \rangle$ with $\alpha \in \mathbb{C}$. It will be a linear combination of the Fock state vectors $\{ |n,n\rangle , |n-1,n-1\rangle , |n+1,n+1\rangle \}$.
EDIT: I forgot to mention my confusion comes from the fact that the matrix representation of the number operator is infinite dimensional, so conventional methods to find the eigenvectors will not work. From my research, it appears as it because of this infinite nature, its completely possible it has no eigenvectors / eigenvalues. Is this correct?
EDIT 2: @Quantum Mechanic
I started to proceed in the manner with which you answered. I didn't quite get the same recursion relation as you. Namely, I got: $${\psi _{m + 1,n + 1}} = \frac{{\left( {\lambda + (1 - G)(n + 1) - Gm} \right)}}{{\sqrt {G(m + 1)(n + 1)(G - 1)} }}{\psi _{m,n}} - \frac{{\sqrt {mn} }}{{\sqrt {(m + 1)(n + 1)} }}{\psi _{m - 1,n - 1}}$$
It's a little bit difficult to figure out how to proceed. I have been trying to simply evaluate the coefficients, find a pattern so I can put in closed form, and then evaluate the infinite sums. Mathematica has been a god-send but I am wondering if there is a better way. If I limit myself to the case there $m=n$ I can extract a (very complicated) pattern of the terms, but any more general case gets confusing. I will keep trying.
EDIT 3: @Quantum Mechanic @Cosmas Zachos
I think I am almost there in regard to understanding what to do, but I am still a bit confused about the matrix and finding the coefficients. I don't think I get what the matrix represents. I know that it should represent quantities such as $\langle \psi_m | \hat{O} |\psi_n \rangle$, but I just can't quite make that make sense because I am confused on what the eigenvectors should be. Let me explain more.
So, we want:
$$\hat{H} |\psi\rangle = \lambda |\psi \rangle = 0$$
The most general state to begin with is $|\psi_0\rangle = \sum_n \psi_n |n\rangle$, (where the subscript $0$ is for the $\lambda = 0$ case) is this the form we are assuming? Or is it simply $|\psi_n\rangle = \psi_n|n,n\rangle = \psi_n|n\rangle$? (where the subscript now represents the $nth$ eigenvector, being more akin the the wiki link given on ladder operators) There was a lot of talk earlier about the requisite states being infinite in nature so I think its the former, but the matrix given doesn't make sense if this is the case, and would not represent $\langle \psi_m | \hat{O} | \psi_n \rangle$. I can show this by the following:
Applying the general state to the right-hand side of the "Hamiltonian" equation gives:
$$\sum_n \psi_n\left[ (\cosh 2g) n |n\rangle + \frac{1}{2} \sinh 2g \left( (n+1)|n+1\rangle + n|n-1\rangle \right) \right] = 0$$
Since we cannot be adding together different kets together to get zero, the only way I can make sense of Quantum Mechanics comment is to calculate an inner product with the same $\lambda = 0$ eigen-vector. This gives:
$$\sum_m\sum_n \psi_m\psi_n\left[ (\cosh 2g) n \delta_{m,n} + \frac{1}{2} \sinh 2g \left( (n+1)\delta_{m,n+1} + n\delta_{m,n-1} \right) \right] = 0$$
Thus, the elements in the matrix should have products of coefficients. So for example, if $m=0$ and $n=1$, this gives: $\frac{1}{2} (\sinh 2g) \psi_0 \psi_1 = 0$. This makes solving for the coefficients more difficult.
Am I off base?? Thanks for all of your help so far!
EDIT 4: @flippiefanus
Here is my calculation showing that $[\hat{b}^{\dagger}\hat{b},\hat{b}]$ isn't equal to a scalar multiple of $\hat{b}$, where I am writing $\hat{b}=\cosh(g)\hat{a}_1+\sinh(g)\hat{a}_2^{\dagger}$:
EDIT 5: Nevermind! I was able to make the calculation work and have deleted my erroneous work.