A formalism like the one you use in your question, once a few mistakes are amended, will allow you to correctly describe the evolution of single-photon states.
In particular, the state of a single photon with two spatial degrees of freedom and two possible polarization states is described as a vector in a four-dimensional space.
A possible convention is to use
$$
|V\rangle_a=(1,0,0,0)^T, \\
|H\rangle_a= (0,1,0,0)^T, \\
|V\rangle_b=(0,0,1,0)^T, \\
|H\rangle_b= (0,0,0,1)^T.
$$
You can then use the evolution matrix you gave to describe the evolution of any input photon through the PBS.
However, this will not work as soon as you have multiple indistinguishable photons as inputs.
The reason is that the space of possible modes of many indistinguishable photons (or, more generally, bosons) is smaller than the tensor product of the spaces of the single photons.
Roughly speaking, this is because if the photons are indistinguishable, states like $|H_a\rangle|V_b\rangle$ and $|V_b\rangle|H_a\rangle$ are actually the same state.
To properly describe the evolution of many-boson states you need to take into account their indistinguishability.
This can be done in several equivalent ways:
1) using second quantization formalism, which automatically takes into account the symmetry properties of the states, 2) using the unitary evolution that properly describes how many-boson basis states evolve, or 3) keep using the standard formalism, but only computing the amplitudes between symmetrized input and output states.
The input coherent state is : $$|\alpha\rangle\otimes|\beta\rangle = \exp\left(-\frac{|\alpha|^2+|\beta|}{2}\right)\exp\left(\alpha \hat a_i^\dagger + \beta \hat b_i^\dagger\right)|0\rangle$$
Then, since $\hat T$ is unitary and $\hat T^\dagger|0\rangle=|0 \rangle $ we have :
\begin{align}
\hat T|\alpha\rangle\otimes|\beta\rangle &= \exp\left(-\frac{|\alpha|^2+|\beta|}{2}\right) \exp\left(\hat T\left(\alpha \hat a_i^\dagger + \beta \hat b_i^\dagger\right)\hat T^\dagger\right)|0\rangle \\
&=\exp\left(-\frac{|\alpha|^2+|\beta|}{2}\right) \exp\left(\alpha \hat a_o^\dagger + \beta \hat b_o^\dagger\right)|0\rangle \\
&=\exp\left(-\frac{|\alpha|^2+|\beta|}{2}\right) \exp\left( (\alpha\cos(\theta)-i\beta\sin(\theta))\hat a_i^\dagger + (\beta\cos(\theta) - i\alpha\sin(\theta)) \hat b_i^\dagger\right)|0\rangle \\
&=|\alpha\cos(\theta)-i\beta\sin(\theta)\rangle\otimes|\beta\cos(\theta) - i\alpha\sin(\theta) \rangle
\end{align}
Best Answer
The operator for the beam splitter is simply a unitary transformation that gives the output of the beam splitter directly from the input. It does not incorporate any time evolution. Therefore it is not a Hamiltonian for the system. The Hamiltonian for such a linear system is trivial, because its time evolution is trivial.