Simulating the evolution of many-boson states My task is to simulate the scheme presented in this paper: https://journals.aps.org/pra/abstract/10.1103/PhysRevA.77.062316.
In this question: Creating an operator for a polarizing beam splitter, I asked how to model two photons interfering at a polarizing beam splitter (PBS), but with four photons, things get even more complicated (which was also mentioned in the reply to the above question).
The scheme: 
The first figure is an illustration of the scheme with four single photons input at A1, A2, B1 and B2, respectively. All four photons are diagonally polarized at input: 
$\left|D\right> = \frac{1}{\sqrt{2}} \left( \left|H\right> + \left|V\right> \right)$, where $\left|V\right>$ is vertically polarized light and $\left|H\right>$ is horizontally polarized light.
A1 and A2 interfere at one PBS, say PBS1, while B1 and B2 interfere at another, PBS2. 
Now, interference also occurs at the rotated PBS where mode A2' and B2' meet. The rotated PBS is just a PBS (PBS3) positioned as the other PBS' but with one half wave plate (HWP) at each input and output ports, respectively, with all four HWPs oriented at $45^{\circ}$.
Projective measurements are done using a HWP and a PBS that are placed just before the detectors.
Now, my questions are: 


*

*It was suggested in my previous question, that I used second quantization in order to take the symmetry properties of the states into consideration. Since I am unfamiliar with this formalism, I would like to hear if anyone have ideas about how to get started, what literature I should have a look at etc. 

*Depending on the method you suggest, which tools should I use to simulate the model? I would prefer to do the simulation in QuTiP (http://qutip.org/) but Matlab or C++ is also an option, as long as your suggestions can help me overcome the problems, I am stuck with.


Any hints and guidance as to how I should get started with this is very much appreciated! (I would love to dig deep into and understand all aspects of the theory behind but as you might have guessed, I am not a theoretical physicist, thus I am lacking some of the math skills that might be needed).
 A: I believe that the first two sections of Tichy et al. 2012 do a nice job at explaining the basics of how the evolution of many indistinguishable bosons works.
Regarding tools to simulate the evolution of many-boson states, it mostly depends on what kind of simulations you want to perform.
For simple things I've used the python code you can find here. It's a very simple package that handles the basic operations.
It uses qutip for some things but is mostly built with just numpy.
This will probably not be enough for more complex operations, and you might be better off building your own code in that case.
One thing to notice is that many-boson spaces get very big very soon with the number of bosons (the Hilbert space of an $n$-boson space grows exponentially with $n$, if the number of modes $m$ is much bigger than $n$), so that simulating the evolution of many bosons over many modes is in general unfeasible (though there might be tricks to do it in special cases).
Another thing to notice is that if you simulate bosons having an internal degree of freedom (like photons with polarization), and you want to use a code such as the one linked above, then you will need to double the space to take into account the polarization, and the operations to perform will need to be adjusted accordingly. A simple example of this is provided in my answer to the other question you mentioned.
