# Minimizing the amount of correlations¶

In this paper by Jiang, Lang, and Caves , the authors show that if one has two qumodes states $$\left|\psi \right\rangle$$ and $$\left|\phi \right\rangle$$ and a beamsplitter $$\text{BS}(\theta)$$, then the only way no entanglement is generated when the beamsplitter acts on the product of the two states

$\left|\Psi \right\rangle = \text{BS}(\theta) \ \left|\psi \right\rangle \otimes \left|\phi \right\rangle,$

is if the states $$\left|\psi \right\rangle$$ and $$\left|\phi \right\rangle$$ are squeezed states along the same quadrature and by the same amount.

Now imagine the following task: > Given an input state $$\left|\psi \right\rangle$$, which is not necessarily a squeezed state, what is the optimal state $$\left|\phi \right\rangle$$ incident on a beamsplitter $$\text{BS}(\theta)$$ together with $$\left|\psi \right\rangle$$ such that the resulting entanglement is minimized?

In our paper we showed that if $$\theta \ll 1$$ the optimal state $$\left|\phi \right\rangle$$, for any input state $$\left|\psi \right\rangle$$, is always a squeezed state. We furthermore conjectured that this holds for any value of $$\theta$$.

Here, we numerically explore this question by performing numerical minimization over $$\left|\phi \right\rangle$$ to find the state that minimizes the entanglement between the two modes.

First, we import the libraries required for this analysis; NumPy, SciPy, TensorFlow, and StrawberryFields.

:

import numpy as np
from scipy.linalg import expm
import tensorflow as tf

:

import strawberryfields as sf
from strawberryfields.ops import *
from strawberryfields.backends.tfbackend.ops import partial_trace


Now, we set the Fock basis truncation; in this case, we choose $$cutoff=30$$:

:

cutoff = 30


## Creating the initial states¶

We define our input state $$\newcommand{ket}{\left|#1\right\rangle}\ket{\psi}$$, an equal superposition of $$\ket{0}$$ and $$\ket{1}$$:

$\ket{\psi}=\frac{1}{\sqrt{2}}\left(\ket{0}+\ket{1}\right)$
:

psi = np.zeros([cutoff], dtype=np.complex128)
psi = 1.0
psi = 1.0
psi /= np.linalg.norm(psi)


We can now define our initial random guess for the second state $$\ket{\phi}$$:

:

phi = np.random.random(size=[cutoff]) + 1j*np.random.random(size=[cutoff])
phi[10:] = 0.
phi /= np.linalg.norm(phi)


Next, we define the creation operator $$\hat{a}$$,

:

a = np.diag(np.sqrt(np.arange(1, cutoff)), k=1)


the number operator $$\hat{n}=\hat{a}^\dagger \hat{a}$$,

:

n_opt = a.T @ a


and the quadrature operators $$\hat{x}=a+a^\dagger$$, $$\hat{p}=-i(a-a^\dagger)$$.

:

# Define quadrature operators
x = a + a.T
p = -1j*(a-a.T)


We can now calculate the displacement of the states in the phase space, $$\alpha=\langle \psi \mid\hat{a}\mid\psi\rangle$$. The following function calculates this displacement, and then displaces the state by $$-\alpha$$ to ensure it has zero displacement.

:

def recenter(state):
alpha = state.conj() @ a @ state
disp_alpha = expm(alpha.conj()*a - alpha*a.T)
out_state = disp_alpha @ state
return out_state


First, let’s have a look at the displacement of state $$\ket{\psi}$$ and state $$\ket{\phi}$$:

:

psi.conj().T @ a @ psi

:

(0.4999999999999999+0j)

:

phi.conj().T @ a @ phi

:

(1.5466747558600655+0.10715327000748875j)


Now, let’s center them in the phase space:

:

psi = recenter(psi)
phi = recenter(phi)


Checking they now have zero displacement:

:

np.round(psi.conj().T @ a @ psi, 9)

:

0j

:

np.round(phi.conj().T @ a @ phi, 9)

:

(6.3e-08+4e-09j)


## Performing the optimization¶

We can construct the variational quantum circuit, using Strawberry Fields:

:

prog = sf.Program(2)
eng = sf.Engine("tf", backend_options={"cutoff_dim": cutoff})

psi = tf.cast(psi, tf.complex64)
phi_var = tf.cast(tf.Variable(phi.real),tf.complex64) \
+ 1j*tf.cast(tf.Variable(phi.imag), tf.complex64)

in_state = tf.einsum('i,j->ij', psi, phi_var)

with prog.context as q:
Ket(in_state) | q
BSgate(np.pi/4, 0) | q

result = eng.run(prog, run_options={"eval": False, "modes": })
state = result.state


Here, we are initializing a TensorFlow variable phi_var representing the initial state of mode q, which we will optimize over. Note that we take the outer product $$\ket{in}=\ket{\psi}\otimes\ket{\phi}$$, and use the Ket operator to initialise the circuit in this initial multimode pure state.

Finally, when we run the engine, we use the argument modes= to return the state of mode q.

We can now extract the density matrix of this state:

:

rhoB = state.dm()


The cost function contains the purity of the reduced density matrix

$\text{Tr}(\rho_B^2),$

and an extra penalty that forces the optimized state to have zero displacement; that is, we want to minimise the value

$\langle \hat{x}\rangle=\text{Tr}(\rho_B\hat{x}).$

Finally, we divide by the $$\text{Tr}(\rho_B)^2$$ so that the state is always normalized.

:

penalty_strength = 10

:

cost = tf.cast(tf.real((tf.trace(rhoB @ rhoB)
-penalty_strength*tf.trace(rhoB @ x)**2
-penalty_strength*tf.trace(rhoB @ p)**2)
/(tf.trace(rhoB))**2), tf.float64)


We can now set up the optimization, to minimise the cost function:

:

optimizer = tf.train.AdamOptimizer()
minimize_op = optimizer.minimize(-cost)

:

sess = tf.Session()
sess.run(tf.global_variables_initializer())


Running the optimization process for 1201 reps:

:

reps = 1201

cost_progress = []

for i in range(reps):

[_, cost_val, ket_val] = sess.run([minimize_op, cost, phi_var])
# Stores cost at each step
cost_progress.append(cost_val)

# Prints progress
if i % 50 == 0:
print("Rep: {} Cost: {}".format(i, cost_val))

Rep: 0 Cost: 0.529472827911377
Rep: 50 Cost: 0.6242330074310303
Rep: 100 Cost: 0.7068203687667847
Rep: 150 Cost: 0.7573347687721252
Rep: 200 Cost: 0.797258198261261
Rep: 250 Cost: 0.8431282043457031
Rep: 300 Cost: 0.8931565284729004
Rep: 350 Cost: 0.9057941436767578
Rep: 400 Cost: 0.9095619916915894
Rep: 450 Cost: 0.9111961722373962
Rep: 500 Cost: 0.9118565917015076
Rep: 550 Cost: 0.9121061563491821
Rep: 600 Cost: 0.9121946096420288
Rep: 650 Cost: 0.9122240543365479
Rep: 700 Cost: 0.9122328162193298
Rep: 750 Cost: 0.9122357964515686
Rep: 800 Cost: 0.9122360944747925
Rep: 850 Cost: 0.9122363328933716
Rep: 900 Cost: 0.9122363328933716
Rep: 950 Cost: 0.9122363328933716
Rep: 1000 Cost: 0.9122363924980164
Rep: 1050 Cost: 0.9122364521026611
Rep: 1100 Cost: 0.9122364521026611
Rep: 1150 Cost: 0.9122365117073059
Rep: 1200 Cost: 0.9122363924980164


We can see that the optimization converges to the optimum purity value of 0.9122365.

## Visualising the optimum state¶

We can now calculate the density matrix of the input state $$\ket{\phi}$$ which minimises entanglement:

$\rho_{\phi} = \ket{\phi}\left\langle \phi\right|$
:

out_rhoB = np.outer(ket_val, ket_val.conj())


Next, we can use the following function to plot the Wigner function of this density matrix:

:

import copy
def wigner(rho, xvec, pvec):
# Modified from qutip.org
Q, P = np.meshgrid(xvec, pvec)
A = (Q + P * 1.0j) / (2 * np.sqrt(2 / 2))

Wlist = np.array([np.zeros(np.shape(A), dtype=complex) for k in range(cutoff)])

# Wigner function for |0><0|
Wlist = np.exp(-2.0 * np.abs(A) ** 2) / np.pi

# W = rho(0,0)W(|0><0|)
W = np.real(rho[0, 0]) * np.real(Wlist)

for n in range(1, cutoff):
Wlist[n] = (2.0 * A * Wlist[n - 1]) / np.sqrt(n)
W += 2 * np.real(rho[0, n] * Wlist[n])

for m in range(1, cutoff):
temp = copy.copy(Wlist[m])
# Wlist[m] = Wigner function for |m><m|
Wlist[m] = (2 * np.conj(A) * temp - np.sqrt(m)
* Wlist[m - 1]) / np.sqrt(m)

# W += rho(m,m)W(|m><m|)
W += np.real(rho[m, m] * Wlist[m])

for n in range(m + 1, cutoff):
temp2 = (2 * A * Wlist[n - 1] - np.sqrt(m) * temp) / np.sqrt(n)
temp = copy.copy(Wlist[n])
# Wlist[n] = Wigner function for |m><n|
Wlist[n] = temp2

# W += rho(m,n)W(|m><n|) + rho(n,m)W(|n><m|)
W += 2 * np.real(rho[m, n] * Wlist[n])

return Q, P, W / 2

:

# Import plotting
from matplotlib import pyplot as plt
%matplotlib inline

:

x = np.arange(-5, 5, 0.1)
p = np.arange(-5, 5, 0.1)
X, P, W = wigner(out_rhoB, x, p)
plt.contourf(X, P, np.round(W,3), cmap="PiYG")

:

<matplotlib.contour.QuadContourSet at 0x7fbc7368cdd8> We see that the optimal state is indeed a (mildly) squeezed state. This can be confirmed by visuallising the Fock state probabilities of state $$\ket{\phi}$$:

:

plt.bar(np.arange(cutoff), height=np.abs(ket_val)**2)

:

<BarContainer object of 30 artists> Finally, printing out the mean number of photons $$\bar{n} = \left\langle \phi \mid \hat{n} \mid \phi\right\rangle$$, as well as the squeezing magnitude $$r=\sinh^{-1}\left(\sqrt{\bar{n}}\right)$$ of this state:

:

nbar = ((ket_val.conj()).T @ n_opt @ ket_val).real
print("mean number of photons =",nbar)
print("squeezing parameter =",np.arcsinh(np.sqrt(nbar)))

mean number of photons = 0.08448536543767962
squeezing parameter = 0.2867190599872425