Creating an operator for a polarizing beam splitter I am trying to make a polarizing beam splitter operator but all my approaches have failed so far. I will try to explain what I have done, and hopefully there are people out there that can point out my mistakes and/or suggest better methods.
I do all calculations in QuTiP, but any help regarding the theory is highly appreciated. 
Let's take the case with two polarized photons input to a polarizing beam splitter. They can be in the states $\left|V\right>$  (vertically polarized), $\left|H\right>$ (horizontally polarized) or a superposition of the two. $\left|H\right>$  is defined as a (2x1) vector, a so called "Quantum Object" in QuTiP: Qobj([[0], [1]])  and $\left|V\right>$ is  Qobj([[1], [0]]). The input state is the tensor product of the two states, with the subscripts denoting the ports, for example: $\left| H \right>_a \left| H \right>_b$.
The inputs are at port a and b, and outputs at port c and d.
Using this, I have constructed this unitary PBS operator:
$$
\hat{O}_{PBS} =
\begin{bmatrix}
    0 & 1 & 0 & 0 \\
    1 & 0 & 0 & 0 \\
    0 & 0 & i & 0 \\
    0 & 0 & 0 & i
\end{bmatrix}
$$
The i's are the phases comming from the reflections, and the rows are ordered in the following way: $\left( \left| H \right>_a, \left| H \right>_b, \left| V \right>_a, \left| V \right>_b \right)$ to $\left( \left| H \right>_c, \left| H \right>_d, \left| V \right>_c, \left| V \right>_d \right)$.
At port a, just before the PBS, I put a half-wave plate (HWP) to rotate the polarization of one of the photons. The operator looks like this:
$$
\hat{O}_{HWP} =
\begin{bmatrix}
    \cos(2\theta) & \sin(2\theta)  \\
    \sin(2\theta) & \cos(2\theta)  \\
\end{bmatrix}
$$
where $\theta$ is the angle with which the HWP is rotated.
So if I input two horizontally polarized photons, the state after the PBS is:
$$
\left| \psi \right> = 
\left(\hat{O}_{HWP} \otimes \mathbb{1} \right) \hat{O}_{PBS}
\left| H \right>_a \left| H \right>_b
$$
($\hat{O}_{HWP}$ only acts on the state at port a)
At port c and d, I measure the output the following way:


*

*I make a density matrix from the output state:
$$
\rho = \left| \psi \right> \left< \psi \right|
$$

*and a measurement operator that measures only a certain polarization at both output ports, c and d, for example $\left| H \right>$:


$$
\hat{M} = \left| H \right> \left< H \right|
$$


*The expectation value, I find by tracing:


$$
m = tr \left( \hat{M} \rho \right)
$$
In my calculations, I change the polarization of the photon at port a by iterating $\theta$ from 0 to $\pi$ and then I evaluate the expectation value at each step.
For $\left| H \right>$ input at port a and b, respectively, and measuring  $\left| H \right>$ at both output ports, I get:

But the weird thing happens when I input $\left| V \right>$ at port a and b, respectively, and measure $\left| V \right>$ at both output ports:

Here the expectation value is constant. 
So my questions are:


*

*What am I doing wrong in my calculations?

*Do I construct the PBS operator the correct way?

*If I were to only do a measurement at port c, I would have to project the state at port d onto $\mathbb{1}$. The problem is the fact that the measurement operator is the outer product of two (2x1) vectors, and I cannot make $\mathbb{1}$ as a (2x1) vector. How can I make a measurement operator, that projects the state after the PBS onto for example $\left| H \right>$ at port c and $\mathbb{1}$ at port d?

*When I have a well functioning PBS operator, how do I include losses in it?


Please let me know if I need to elaborate on something. Thank you very much in advance!
 A: 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.
