# 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:

1. 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.
2. 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).

• Consider starting here: arxiv.org/abs/1711.00080 and change the indices $i$ and $j$ by $i,\alpha$ and $j,\beta$ where $\alpha,\beta$ refer to polarization. Alternatively, work with $\hat a^\dagger_\alpha$ and $\hat b^\dagger_\beta$ for the same reasons. Aug 6, 2018 at 13:28
• Second quantization is an approach that was developed for quantum field theory. However, you don't need QFT for this problem. The quantum information community (in particular those working with so-called continuous variables) use this term to refer to an approach based on the use of creation and annihilation operators instead of bra and ket vector. This is especially convenient for coherent states (or more generally, Gaussian states). So the idea is to represent the processes as operators expressed in term of creation and annihilation operators. Aug 7, 2018 at 4:22

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).