Simulating a TMSV and environmental effects
Hello! I am trying to simulate just a TMSV with one of the modes with Thermal Loss
and then a rotation on this mode. To do this I use strawberry fields reduced_gaussian([0,1]) function to obtain a covariance matrix and a vector of expected values. Using these I input them in a multivariate distribution sampling function; multivariate_normal(mu,cov,sample_size), to get samples. I do this because it is faster if I want to sample like 1e6 samples and I am assuming gaussian systems. Also I sample from the Husimi covariance matrix just to emulate heterodyne measurments.
import matplotlib.pyplot as plt
import numpy as np
from sympy import *
from sympy.physics.quantum import TensorProduct
import scipy.constants as constants
import strawberryfields as sf
from strawberryfields.ops import *
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
import scipy.stats as stats
import scipy.constants as constants
builds the circuit:
def Two_mode(eta = 1, r2 =0 , gamma = 0, phi1 = 0, phi2 = 0, nbar = 0, hbar = 2):
"This function builds the circuit I want to simulate"
eng = sf.Engine("gaussian")
prog=sf.Program(2)
sf.hbar=hbar
with prog.context as q:
S2gate(r2,gamma) | (q[0],q[1])
ThermalLossChannel(eta,nbar) | (q[0])
ThermalLossChannel(1,0) | (q[1])
Rgate(phi1) | (q[0])
Rgate(phi2) | (q[1])
result=eng.run(prog)
state=result.state
return state.reduced_gaussian([0,1])
Samples from multivariate normal distribution:
def Statistics2( et = 1, r2squeeze = 0, gam = 0, ph1 = 0, ph2 = 0, nth = 0, h_bar = 2, sample_size = int(1e6)):
"Samples from multivariate normal distribution (Husimi)"
#Gives wigner matrix and mean vector:
mu, cov = Two_mode(eta=et, r2=r2squeeze, gamma=gam, phi1=ph1 ,phi2=ph2 ,nbar=nth ,hbar=h_bar)
#Theoretical Husimi:
cov += 0.5*h_bar*np.eye(4)
#Sampled Husimi
rng=np.random.default_rng()
samples=rng.multivariate_normal(mu,cov,sample_size)
X1=np.array([samples[i][0] for i in range(sample_size)])
X2=np.array([samples[i][1] for i in range(sample_size)])
P1=np.array([samples[i][2] for i in range(sample_size)])
P2=np.array([samples[i][3] for i in range(sample_size)])
covH = np.matrix(np.cov((X1,P1,X2,P2), bias=True))
meanX1 = np.mean(X1)
meanP1 = np.mean(P1)
meanX2 = np.mean(X2)
meanP2 = np.mean(P2)
return covH, meanX1, meanP1, meanX2, meanP2
Calculates the average number of thermal photons
def nbar_thermal(T,f):
lamb = (constants.hbar*2*np.pi*f)/(constants.Boltzmann*T)
return(1/(np.exp(lamb) - 1))
Example (12 seconds):
T = 300
f = 10e9
transmittance = 0.3
rsqueeze = np.linspace(0,3.001,100)
nbar = nbar_thermal(T,f)
sample_size=int(1e4)
phi1 = np.pi*np.array([1/4,1/3,1/2,1])
Sampled_covH = np.empty([len(phi1),len(rsqueeze),4,4])
for i, p in enumerate(phi1):
for j, r in enumerate(rsqueeze):
Sampled_covH[i,j], meanX1, meanP1, meanX2, meanP2 = Statistics2( et = transmittance, r2squeeze = r, ph1 = p, nth = nbar, sample_size = sample_size)
Now I am doing this because I want to compare these results with theoretical covariance matrices I derived using the symplectic transformation formalism. Now I will show the covariance matrix I obtained myself by expressing all the symplectic matrices for differents operations respecting the conventions in strawberry fields and also by using thewalrus.symplectic. Both of these methods gave me the same thing:
V_{\text{Husimi}}=\frac{\hbar}{2}
\begin{pmatrix}
\eta \text{cosh}(2r) + (1-\eta)(2 n_{th} + 1) + 1 & 0 & \sqrt{\eta} \text{cos}(\phi)\text{sinh}(2r) & \sqrt{\eta} \text{sin}(\phi)\text{sinh}(2r)\\
0 & \eta \text{cosh}(2r) + (1-\eta)(2 n_{th} + 1) + 1 & \sqrt{\eta} \text{sin}(\phi)\text{sinh}(2r) & -\sqrt{\eta} \text{cos}(\phi)\text{sinh}(2r)\\
\sqrt{\eta}\text{cos}(\phi)\text{sinh}(2r) & \sqrt{\eta} \text{sin}(\phi)\text{sinh}(2r) & \text{cosh}(2r) + 1 & 0\\
\sqrt{\eta} \text{sin}(\phi)\text{sinh}(2r) & - \sqrt{\eta} \text{cos}(\phi)\text{sinh}(2r) & 0 & \text{cosh}(2r) + 1
\end{pmatrix}
Here \eta is the transmittance or transmissivity of the beamsplitter in the ThermalLossChannel. I set \hbar = 2 in the functions above. Now I compared was the variances on the diagonal of V_{\text{Husimi}} with the variances I would get from the MonteCarlo simulations I did above.
Comparison with theory
"Theoretical variances of mode 1 and mode 2 in the covariance matrix"
def var_mode1(eta,r,nth,hbar=2):
return hbar*0.5*(eta*np.cosh(2*r) + (1-eta)*(2*nth + 1) + 1)
def var_mode2(eta,r,nth,hbar=2):
return hbar*0.5*(np.cosh(2*r) + 1)
"We calculate the variances"
#theoretical
varmode1=[[var_mode1(eta=transmittance, r=r, nth=nbar) for r in rsqueeze] for p in phi1]
varmode2=[[var_mode2(eta=transmittance, r=r, nth=nbar) for r in rsqueeze] for p in phi1]
#Sampled
sampled_VarI1=[[Sampled_covH[i,j][0,0] for j in range(len(rsqueeze))] for i in range(len(phi1))]
sampled_VarQ1=[[Sampled_covH[i,j][1,1] for j in range(len(rsqueeze))] for i in range(len(phi1))]
sampled_VarI2=[[Sampled_covH[i,j][2,2] for j in range(len(rsqueeze))] for i in range(len(phi1))]
sampled_VarQ2=[[Sampled_covH[i,j][3,3] for j in range(len(rsqueeze))] for i in range(len(phi1))]
colors_plot = ['slateblue', 'tab:red', 'limegreen','black']
colors_scatter = ['cornflowerblue', 'lightcoral', 'springgreen','grey']
legend = [r'$\frac{\pi}{4}$', r'$\frac{\pi}{3}$', r'$\frac{\pi}{2}$', r'$\pi$']
fig, axs = plt.subplots(2,2,figsize=(15, 10))
for i in range(len(phi1)):
axs[0,0].scatter(rsqueeze, sampled_VarI1[i], color=colors_scatter[i], s=20)
axs[0,0].plot(rsqueeze, varmode1[i], label = r'$\phi$' +' = ' + legend[i], color=colors_plot[i], linewidth=3, alpha=0.9)
axs[0,0].set_xlabel('squeezing parameter ' + r'$r_{2}$', fontsize=15)
axs[0,0].set_ylabel( r'$VAR[I_{1}]$', fontsize=15)
axs[0,0].set_title('Variance of ' + r'$I_{1}$')
axs[0,0].legend(fontsize=12)
axs[0,1].scatter(rsqueeze, sampled_VarQ1[i], color=colors_scatter[i], s=20)
axs[0,1].plot(rsqueeze, varmode1[i], label = r'$\phi$' +' = ' + legend[i], color=colors_plot[i], linewidth=3, alpha=0.9)
axs[0,1].set_xlabel('squeezing parameter ' + r'$r_{2}$', fontsize=15)
axs[0,1].set_ylabel( r'$VAR[Q_{1}]$', fontsize=15)
axs[0,1].set_title('Variance of ' + r'$Q_{1}$')
axs[0,1].legend(fontsize=12)
axs[1,0].scatter(rsqueeze, sampled_VarI2[i], color=colors_scatter[i], s=20)
axs[1,0].plot(rsqueeze, varmode2[i], label = r'$\phi$' +' = ' + legend[i], color=colors_plot[i], linewidth=3, alpha=0.9)
axs[1,0].set_xlabel('squeezing parameter ' + r'$r_{2}$', fontsize=15)
axs[1,0].set_ylabel( r'$VAR[I_{2}]$', fontsize=15)
axs[1,0].set_title('Variance of ' + r'$I_{2}$')
axs[1,0].legend(fontsize=12)
axs[1,1].scatter(rsqueeze, sampled_VarQ2[i], color=colors_scatter[i], s=20)
axs[1,1].plot(rsqueeze, varmode2[i], label = r'$\phi$' +' = ' + legend[i], color=colors_plot[i], linewidth=3, alpha=0.9)
axs[1,1].set_xlabel('squeezing parameter ' + r'$r_{2}$', fontsize=15)
axs[1,1].set_ylabel( r'$VAR[Q_{2}]$', fontsize=15)
axs[1,1].set_title('Variance of ' + r'$Q_{2}$')
axs[1,1].legend(fontsize=12)
plt.tight_layout()
As you will see the variances for the quadratures of the second mode (idler), which was not applied the same ThermalLossChannel as the first mode, is off from the theoretical value from the covariance matrix. Now I found the term to correct this but I don’t understand where it comes from. Here it is:
"Theoretical variances of mode 1 and mode 2 corrected"
def var_mode1(eta,r,nth,hbar=2):
return hbar*0.5*(eta*np.cosh(2*r) + (1-eta)*(2*nth + 1) + 1)
def var_mode2(eta,r,nth,hbar=2):
return hbar*0.5*(np.cosh(2*r) + (1-eta)*(2*nth + 1) + 1)
"We calculate the variances"
#theoretical
varmode1=[[var_mode1(eta=transmittance, r=r, nth=nbar) for r in rsqueeze] for p in phi1]
varmode2=[[var_mode2(eta=transmittance, r=r, nth=nbar) for r in rsqueeze] for p in phi1]
#Sampled
sampled_VarI1=[[Sampled_covH[i,j][0,0] for j in range(len(rsqueeze))] for i in range(len(phi1))]
sampled_VarQ1=[[Sampled_covH[i,j][1,1] for j in range(len(rsqueeze))] for i in range(len(phi1))]
sampled_VarI2=[[Sampled_covH[i,j][2,2] for j in range(len(rsqueeze))] for i in range(len(phi1))]
sampled_VarQ2=[[Sampled_covH[i,j][3,3] for j in range(len(rsqueeze))] for i in range(len(phi1))]
colors_plot = ['slateblue', 'tab:red', 'limegreen','black']
colors_scatter = ['cornflowerblue', 'lightcoral', 'springgreen','grey']
legend = [r'$\frac{\pi}{4}$', r'$\frac{\pi}{3}$', r'$\frac{\pi}{2}$', r'$\pi$']
fig, axs = plt.subplots(2,2,figsize=(15, 10))
for i in range(len(phi1)):
axs[0,0].scatter(rsqueeze, sampled_VarI1[i], color=colors_scatter[i], s=20)
axs[0,0].plot(rsqueeze, varmode1[i], label = r'$\phi$' +' = ' + legend[i], color=colors_plot[i], linewidth=3, alpha=0.9)
axs[0,0].set_xlabel('squeezing parameter ' + r'$r_{2}$', fontsize=15)
axs[0,0].set_ylabel( r'$VAR[I_{1}]$', fontsize=15)
axs[0,0].set_title('Variance of ' + r'$I_{1}$')
axs[0,0].legend(fontsize=12)
axs[0,1].scatter(rsqueeze, sampled_VarQ1[i], color=colors_scatter[i], s=20)
axs[0,1].plot(rsqueeze, varmode1[i], label = r'$\phi$' +' = ' + legend[i], color=colors_plot[i], linewidth=3, alpha=0.9)
axs[0,1].set_xlabel('squeezing parameter ' + r'$r_{2}$', fontsize=15)
axs[0,1].set_ylabel( r'$VAR[Q_{1}]$', fontsize=15)
axs[0,1].set_title('Variance of ' + r'$Q_{1}$')
axs[0,1].legend(fontsize=12)
axs[1,0].scatter(rsqueeze, sampled_VarI2[i], color=colors_scatter[i], s=20)
axs[1,0].plot(rsqueeze, varmode2[i], label = r'$\phi$' +' = ' + legend[i], color=colors_plot[i], linewidth=3, alpha=0.9)
axs[1,0].set_xlabel('squeezing parameter ' + r'$r_{2}$', fontsize=15)
axs[1,0].set_ylabel( r'$VAR[I_{2}]$', fontsize=15)
axs[1,0].set_title('Variance of ' + r'$I_{2}$')
axs[1,0].legend(fontsize=12)
axs[1,1].scatter(rsqueeze, sampled_VarQ2[i], color=colors_scatter[i], s=20)
axs[1,1].plot(rsqueeze, varmode2[i], label = r'$\phi$' +' = ' + legend[i], color=colors_plot[i], linewidth=3, alpha=0.9)
axs[1,1].set_xlabel('squeezing parameter ' + r'$r_{2}$', fontsize=15)
axs[1,1].set_ylabel( r'$VAR[Q_{2}]$', fontsize=15)
axs[1,1].set_title('Variance of ' + r'$Q_{2}$')
axs[1,1].legend(fontsize=12)
plt.tight_layout()
From this you should see that everything works. So I don’t understand where this additive thermal noise comes from for the second mode. There is a dependance with 1-\eta but I’m not sure if that makes because I don’t see how thermal photons from the first mode can propagate to the second mode. Is itreally just because of the entanglement? And if so, why can I not see it when I calculate the theoretical covariance matrix?
For the sake of having everything in one place I will show how I found the covariance matrix:
Symplectic transformations
"Symbolic notations are used in this notebook"
eta,alpha,theta,r,gamma,phi,nth,hbar=symbols('\u03B7 \u03B1 \u0398 r \u03B3 \u03D5 n\u209C\u2095 \u045B ')
"Let's define some of the symplectic matrices"
#Beam splitter
def BS(eta):
S_BS=Matrix([[ sqrt(eta) , 0 , -sqrt(1-eta) , 0 ],
[ 0 , sqrt(eta) , 0 , -sqrt(1-eta) ],
[ sqrt(1-eta) , 0 , sqrt(eta) , 0 ],
[ 0 , sqrt(1-eta) , 0 , sqrt(eta) ]])
return S_BS
#single_mode squeezing with gamma=0
def SQ1(r):
SQ1=Matrix([[ exp(-r) , 0],
[ 0 , exp(r) ]])
return SQ1
#rotation of a single mode (anticlockwise)
def Srot(phi):
S_rot=Matrix([[ cos(phi) , -sin(phi)],
[ sin(phi) , cos(phi) ]])
return S_rot
#Two mode squeezing with gamma=0
def SQ2(r):
SQ2=Matrix([[ cosh(r) , 0 , sinh(r) , 0 ],
[ 0 , cosh(r) , 0 , -sinh(r) ],
[ sinh(r) , 0 , cosh(r) , 0 ],
[ 0 , -sinh(r) , 0 , cosh(r) ]])
return SQ2
#Covariance for vacuum which is the same for a coherent state
def CM_vac():
V=Matrix([[hbar/2,0],[0,hbar/2]])
return V
#Covariance for a Thermal state
def CM_Thermal(nth):
V=(2*nth + 1)*Matrix([[hbar/2,0],[0,hbar/2]])
return V
#Covariance for a single-mode squeezed state gamma=0
def CM_SQ1(r):
V=Matrix([[ (hbar/2)*exp(-2*r) , 0],
[ 0 , (hbar/2)*exp(2*r) ]])
return V
#Covariance matrix for a two-mode vacuum system
def CM_2vac():
V_vac=CM_vac()
V2=Matrix.diag(V_vac,V_vac)
return V2
V_system = diag(V2vac,CM_Thermal(nth),CM_vac())
display(Matrix(V_system))
g4 = Matrix([[ 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ],
[ 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ],
[ 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 ],
[ 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 ],
[ 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ],
[ 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 ],
[ 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 ],
[ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 ]])
Squeezed2_sys = diag(SQ2(r), np.eye(4))
S_BS_sys = diag(BS(eta),np.eye(4))
S_rotmode1 = diag(Srot(phi),np.eye(6))
V_withnoise = S_rotmode1 @ (g4.T @ (S_BS_sys @ (g4 @ (Squeezed2_sys @ V_system @ Squeezed2_sys.T) @ g4.T) @ S_BS_sys.T) @ g4) @ S_rotmode1.T
display(Matrix(V_withnoise))
V_systemA = V_withnoise[:4,:4]
display(Matrix(V_systemA))
from which we use identities such as \text{cosh}(2x) = \text{sinh}^{2}(x) + \text{cosh}^{2}(x) and also \text{sinh}(2x) = 2\text{sinh}(x)\text{cosh}(x) to simplify. I also verified that S \Omega S^{T}=\Omega.
Sorry if this is a very long question I am just really curious as to what I don’t understand about this problem.