Number Operator Eigenvectors of Amplifier 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.
 A: I might as well collect my notational simplification in a non-answer, instead of in garlands of comments...
Define $\sqrt{G}\equiv \cosh g \leadsto \sqrt{G-1}= \sinh g$. Moreover, focus on the states $|n\rangle\equiv |n,n\rangle$. Your "hamiltonian" keeps you in this class, and notation convention,
$$
(b^\dagger b - \sinh^2 g)|n\rangle \\ = (\cosh 2g)~ n |n\rangle + \tfrac{1}{2} \sinh 2g \bigl ((n+1)|n+1\rangle +n|n-1\rangle\bigr ) ,
$$
which you may write as the standard Fock infinite-dimensional matrix like this and contemplate the eigenvectors of,
$$
M= \begin{pmatrix}
0& \tfrac{1}{2} \sinh 2g & 0 & ...\\ \tfrac{1}{2} \sinh 2g&
\cosh 2g &   \sinh 2g &0& ... \\
0&   \sinh 2g &  2 \cosh 2g & \tfrac{3}{2} \sinh 2g & 0& ... \\ 0&0&...
\end{pmatrix}.
$$

Edit in response to question's 3rd edit:
To find the eigenvectors of the matrix M, you write the eigenvalue conditions, which  specify a recursive determination of their coefficients, given the escalated structure of the Matrix, and so the (infinity) of eigenvectors will all be infinite-dimensional. (But, unlike the simple eigenvectors of $\hat X$ in Fock space, $ \frac{1}{\pi^{1/4}} e^{-x^2/2} e^{\sqrt{2} xa^\dagger} e^{-a^{\dagger ~2}/2} =\frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle  
$, they may be very complicated...)
The recursive determination of the eigenvectors goes as follows,
$$
M|\psi\rangle= \lambda |\psi\rangle, ~~~ |\psi\rangle = \sum_n c_n|n\rangle ~~\leadsto \\
\lambda c_0=  \tfrac{1}{2} \sinh 2g ~~c_1, \\
\lambda c_1 =  \tfrac{1}{2} \sinh 2g ~~c_0 +  
\cosh 2g ~~ c_1 +    \sinh 2g ~~ c_2, \\ ... ,
$$
so you may recursively solve for $c_1$ in terms of $\lambda, ~g$ and $c_0$, then for $c_2$ in terms of the above, etc... Conceptually, the infinity of $c_n$ are uniquely determinable in terms of  $\lambda, ~g$ and $c_0$, given the stepwise nature of the matrix M ! If the answer were a simple closed form, you could advance to issues of convergence, etc..., as in the linked example., but only then: don't let such qualms stop you.
For a zero eigenvalue, taking $c_0=1$ to work the normalization later, you find $(1,0,-1/2, 2/3 ~ \coth 2g  , 3/8-\coth^2 2g,...)^T$, which is to say
$$
|\psi\rangle= 
|0\rangle  -\tfrac{1}{2}|2\rangle +  \tfrac{2}{3} ~ \coth 2g |3\rangle +(\tfrac{3}{8}-\coth^2 2g)|4\rangle +...
$$
A: We can apply this operator to a general state $$|\psi\rangle=\sum_{mn}\psi_{mn}|m,n\rangle$$ and determine its eigensystem. Expanding, we find
\begin{aligned}
b^\dagger b|\psi\rangle&=\sum_{mn}\psi_{mn}\left\{\left[Gm+(G-1)n\right]|m,n\rangle+\sqrt{G(G-1)}\sqrt{(m+1)(n+1)}|m+1,n+1\rangle+\sqrt{G(G-1)}\sqrt{mn}|m-1,n-1\rangle\right\}\\
&=\sum_{mn}|m,n\rangle\left\{\left[Gm+(G-1)n\right]\psi_{mn}+\sqrt{G(G-1)}(\sqrt{(m+1)(n+1)}\psi_{m+1,n+1}+\sqrt{mn}\psi_{m-1,n-1})\right\}.
\end{aligned} The eigenvalue equation now requires, for all $m$ and $n$,
$$\lambda\psi_{mn}=\left\{\left[Gm+(G-1)n\right]\psi_{mn}+\sqrt{G(G-1)}(\sqrt{(m+1)(n+1)}\psi_{m+1,n+1}+\sqrt{mn}\psi_{m-1,n-1})\right\}.$$ This looks like a recurrence relation:
$$\psi_{m+1,n+1}=\frac{\lambda-Gm-(G-1)n}{\sqrt{G(G-1)}\sqrt{mn}}\psi_{mn}-\frac{\sqrt{mn}}{\sqrt{(m+1)(n+1)}}\psi_{m-1,n-1}.$$ The recurrence relation immediately generates a set of states starting from $m=0$ or $n=0$, of the form $a|0,n_0\rangle+b|1,n_0+1\rangle+c|2,n_0+2\rangle\cdots$, etc. The series does not obviously terminate, so the next step would be to determine whether it is normalizable; i.e., the coefficients must decrease in magnitude sufficiently quickly. If they do, then the infinite-dimensional states are legitimate eigenstates for each eigenvalue $\lambda$; if not, one could say that there are no normalizable eigenstates for this operator.
A: First, I'll consider the Bogoliubov transformation
$$ \hat{a} \rightarrow \hat{b} = \sqrt{G} \hat{a} + \sqrt{G-1} \hat{a}^{\dagger} , $$
which produces a set of Bogoliubov operators that obeys same commutation relation as the original set. Therefore, the Bogoliubov transformation can be seen as an isomorphism, so that the eigenstates of operators defined in terms of $\hat{a}$ are analogous to those of operators defined in terms of $\hat{b}$. Let's then review eigenstates of operators defined in terms of $\hat{a}$ and apply the result to those defined in terms of $\hat{b}$.
The eigenstates of $\hat{a}$ are the coherent states
$$ \hat{a}|\alpha\rangle = \hat{a}\hat{D}(\alpha)|\text{vac}\rangle = |\alpha\rangle \alpha . $$
It follows because
$$ \hat{D}^{\dagger}(\alpha)\hat{a}\hat{D}(\alpha) = \hat{a} + \alpha . $$
Provided that we have at least one solution, the displacement operators then serve to produce all the other solutions. So we can ignore the displacement operator at first and only look for one solution, which in this case is the vacuum state $\hat{a}|\text{vac}\rangle = 0$.
Next, we review the eigenstates of the number operator:
$$ \hat{n} |n\rangle = \hat{a}^{\dagger}\hat{a} |n\rangle = |n\rangle n . $$
This time it follows that
$$ [\hat{n},\hat{a}^{\dagger}]=\hat{a}^{\dagger} . $$
Therefore, if we have one solution, we can generate the other with the aid of the ladder operators. Again, we only need one solution, which is again the vacuum state, because $\hat{n}|\text{vac}\rangle = 0$. So, the vacuum state is common between the two sets of eigenstates for these two sets.
Now we can go and do the same for the Bogoliubov operator. The eigenstates that are analogous to those for the annihilation operators can again be generated with the displacement operator provided that we have at least one solution, because
$$ \hat{D}^{\dagger}(\alpha)\hat{b}\hat{D}(\alpha) = \hat{b} + \alpha \sqrt{G} \alpha + \sqrt{G-1} \alpha^* . $$
The one solution would be the equivalent of the vacuum state that gives $\hat{b}|\text{vac}\rangle_b = 0$. Here $|\text{vac}\rangle_b$ is a squeezed vacuum state associated with the Bogoliubov transformation.
The eigenstates that are analogous to the number states are by analogy generated the Bogoliubov operators, provided we have at least one solution. By analogy it should be the same equivalent of the vacuum state, which is the squeezed vacuum state. In other words, the eigenstate solutions of $\hat{n}_b = \hat{b}^{\dagger}\hat{b}$ are given by
$$ |n\rangle_b = \frac{1}{\sqrt{n!}} \left(\hat{b}^{\dagger}\right)^n \hat{S}_b |\text{vac}\rangle , $$
where $\hat{S}_b$ is the squeezing operator associated with the Bogoliubov transformation that produces $\hat{b}$.
The relationship between the Bogoliubov operator and the squeezing operator is given by
$$ \hat{S}_b^{\dagger} \hat{a} \hat{S}_b = \hat{b} , $$
where
$$ \hat{S}_b 
= \exp(\tfrac{1}{2}\zeta\hat{a}^{\dagger}\hat{a}^{\dagger} 
- \tfrac{1}{2}\zeta^*\hat{a}\hat{a}) . $$
It implies that
$$ \sqrt{G} = \cosh(|\zeta|) , ~~~ \sqrt{1-G} = \sinh(|\zeta|) . $$
Now to generalize this result to the case of two-mode Bogoliubov transformations, we only need to replace the squeezing operator by the equivalent two-mode squeezing operator.
Note that if $\hat{b} = \sqrt{G} \hat{a}_1 + \sqrt{G-1} \hat{a}_2^{\dagger}$, then we still have
$$ [\hat{b},\hat{b}^{\dagger}]
=G[\hat{a}_1,\hat{a}_1^{\dagger}] 
+ (G-1)[\hat{a}_2^{\dagger},\hat{a}_2] = 1 . $$
Moreover,
$$ [\hat{a}_1,\tfrac{1}{2}\zeta\hat{a}_1^{\dagger}\hat{a}_2^{\dagger} 
- \tfrac{1}{2}\zeta^*\hat{a}_1\hat{a}_2]
=[\hat{a}_1,\hat{a}_1^{\dagger}]\tfrac{1}{2}\zeta\hat{a}_2^{\dagger} 
= \tfrac{1}{2}\zeta\hat{a}_2^{\dagger}. $$
and
$$ [\hat{a}_2^{\dagger},\tfrac{1}{2}\zeta\hat{a}_1^{\dagger}\hat{a}_2^{\dagger} 
- \tfrac{1}{2}\zeta^*\hat{a}_1\hat{a}_2]
= - \tfrac{1}{2}\zeta^*\hat{a}_1[\hat{a}_2^{\dagger},\hat{a}_2] 
= \tfrac{1}{2}\zeta^*\hat{a}_1. $$
Using these two relationships repeatedly, one can obtained the equivalent expression for
$$ \hat{S}_b^{\dagger} \hat{a}_1 \hat{S}_b = \hat{b} , $$
where
$$ \hat{S}_b 
= \exp(\tfrac{1}{2}\zeta\hat{a}_1^{\dagger}\hat{a}_2^{\dagger} 
- \tfrac{1}{2}\zeta^*\hat{a}_1\hat{a}_2) . $$
