How to simulate a beam splitter in QuTiP?

Question

I'm trying to mix two states on a beam splitter whose unitary operator is $$\hat{U}=e^{i\theta \left( \hat{a}^†\hat{b} + \hat{a}\hat{b}^† \right)},$$ where $\hat{a}$, $\hat{b}$ are annihilation operators in mode A and B and $\cos^2{\theta}=T$ and T is transmissivity of the beam splitter. Given an input density matrix $\hat{\rho}_\text{in}=\hat{\rho}_\text{A}\otimes\hat{\rho}_\text{B}$, the output of this operation is then $\hat{\rho}_\text{out}=\hat{U}\hat{\rho}_\text{in}\hat{U}^†$. I know what properties the output state should have given that the covariance matrix $\gamma_\text{in}$ and mean vector $\mu_\text{in}$ are transformed using $$B=\begin{pmatrix} t&0&r&0\\0&t&0&r\\-r&0&t&0\\0&-r&0&t\end{pmatrix}$$, where $t=\sqrt{T}$, $r=\sqrt{1-T}$, as $\gamma_\text{out}=B \gamma_\text{in} B^\top$, $\mu_\text{out}=B\mu_\text{in}$.

Since I didn't find an operator for beam splitter in QuTiP, I tried to calculate it using this code:

def beamsplitter(dim, T):
    bs = (
            complex(0, np.arcsin(np.sqrt(1-T))) *
            (qt.tensor(qt.create(dim), qt.qeye(dim)) * qt.tensor(qt.qeye(dim), qt.destroy(dim)) +
             qt.tensor(qt.destroy(dim), qt.qeye(dim)) * qt.tensor(qt.qeye(dim), qt.create(dim)))
    ).expm()
    return bs

Then I tried to test this, using a vacuum in mode A and a coherent state with amplitude $\alpha=a+ib$. Using the transformation above on the input mean vector, which is $$\mu_\text{in}=\begin{pmatrix} 0\\0\\a\sqrt{2}\\b\sqrt{2} \end{pmatrix},$$ I expect to get $$\mu_\text{out}=\begin{pmatrix}ra\sqrt{2}\\rb\sqrt{2}\\ta\sqrt{2}\\tb\sqrt{2}\end{pmatrix}$$ as a result. When I tried to calculate this

def mean_vector(dim, mix):
    A = mix.ptrace(0)
    B = mix.ptrace(1)
    m1 = qt.expect(qt.position(dim), A)
    m2 = qt.expect(qt.momentum(dim), A)
    m3 = qt.expect(qt.position(dim), B)
    m4 = qt.expect(qt.momentum(dim), B)
    return np.array([[m1], [m2], [m3], [m4]])
dim = 50
alfa = complex(-2, -2)
bs = beamsplitter(dim, 0.5)
mix = bs * qt.tensor(qt.fock_dm(dim, 0), qt.coherent_dm(dim, alfa)) * bs.dag()
print(mean_vector(dim, mix))

I found out that there is an issue somewhere - when both real and imaginary parts of $\alpha$ are given, then one of the 4 mean values will havea flipped sign suggesting some rotation in one of the two modes. When only real or imaginary part is given, then one mode is usually displaced correctly, but in the other one the quadratures will be flipped, so displacing $\hat{x}$ when I expect to displace $\hat{p}$ and vice versa. I cannot think of any explanation of this behaviour and would be very thankful if anyone could tell me what I did wrong or how I could fix it.

How to simulate a beam splitter in QuTiP?

0 Answers0