Mathematics and linear algebra

(*Symplectic Structure*)
J[m_]:=Join[Join[Table[0,{i,m},{j,m}],Table[-Boole[i==j],{i,m},{j,m}],2],Join[Table[Boole[i==j],{i,m},{j,m}],Table[0,{i,m},{j,m}],2]];

(*Useful matrices*)
(**Identity matrix**)
Id[m_]:=IdentityMatrix[2m];
(**Projector on vector g**)
P[g_]:=TensorProduct[g,g];

(*Partitioner*)
(**Useful tool for calculating pair partitions that are used in truncated correlation functions**)
idx[{a_,b_}]:={{a,b}}
idx[list_List]:=Flatten[Function[{row},Join[First@row,#]&/@Last@row]/@(({#,idx[Complement[list,#]]}&)/@(list[[{1,#}]]&/@Range[2,Length@list])),1]
idx[n_?EvenQ]:=idx[Range@n]
part[list_List/;EvenQ[Length@list]]:=Fold[Partition,list[[Flatten[idx[Length@list]]]],{2,Length@list/2}];

For the wigner function

Gaussian states

(*Create some useful covariance matrices*)
(**Squeezed vacuum**)
SqueezedVacuum[Squeezing_List]:=Module[{s,v,m},
m=Length[Squeezing];
s=Join[Table[10^(Squeezing[[j]]/10),{j,m}],Table[10^(-Squeezing[[j]]/10),{j,m}]];
v=DiagonalMatrix[s];
v];

(*The Wigner function*)
W0[β_,V_]:=1/(2^(Length[V]/2)*Pi^(Length[V]/2)*Sqrt[Det[V]])*Exp[-1/2*β.Inverse[V].β];

(*The Wigner function for a subset of modes listed in flist*)
W0Red[β_,V_,flist_]:=Module[{Vsub,Wig,m,p,U,int,bla},m=Length[V]/2;U=Join[flist,Table[J[m].flist[[k]],{k,Length[flist]}]];
Vsub=Table[U[[x]].(V.U[[y]]),{x,Length[U]},{y,Length[U]}];
Wig=W0[β,Vsub,{{1,{1,0}}}];
Return[Wig];];

Photon subtraction

(*We typically follow the notation of the papers*)

(**Compute the matrix A**)
(***Note that we have a mixture of g's. The format should be gmix={{p1,g1},{p2,g2},...} here the pi's sum to 1 and gi should be a 2m dimensional real vector.***)
As[V_,gmix_]:=2*(V-IdentityMatrix[Length[V]]).(Sum[gmix[[i,1]]*(P[gmix[[i,2]]]+P[J[Length[V]/2].gmix[[i,2]]]),{i,Length[gmix]}]) .(V-IdentityMatrix[Length[V]])/(Sum[gmix[[i,1]]*Tr[(V-IdentityMatrix[Length[V]]).(P[gmix[[i,2]]]+P[J[Length[V]/2].gmix[[i,2]]])],{i,Length[gmix]}]);

(**Wigner function**)
Ws[β_,V_,gmix_]:=1/(2^(Length[V]/2+1)*Pi^(Length[V]/2)*Sqrt[Det[V]])*(β.Inverse[V].As[V,gmix].Inverse[V].β - Tr[Inverse[V].As[V,gmix]]+2)*Exp[-1/2*β.Inverse[V].β];
(**Wigner function for a subset of modes, grouped in flist**)
WsRed[β_,V_,gmix_,flist_]:=Module[{Vsub,Wig,Asub,m,p,U,int,bla},m=Length[V]/2;U=Join[flist,Table[J[m].flist[[k]],{k,Length[flist]}]];
Vsub=Table[U[[x]].(V.U[[y]]),{x,Length[U]},{y,Length[U]}];
Asub=Table[U[[x]].(As[V,gmix].U[[y]]),{x,Length[U]},{y,Length[U]}];
Wig=1/(2)*(β.Inverse[Vsub].Asub.Inverse[Vsub].β - Tr[Inverse[Vsub].Asub]+2)*W0[β,Vsub];
Return[Wig];];

(**The probability density function to obtain the result x upon measuring in mode measuredMode**)
MeasurementProbability[x_,V_,gmix_,measuredMode_]:=Module[{β,Vsub,Wig,Asub,m,p,U,int,bla},m=Length[V]/2;
Vsub=measuredMode.(V.measuredMode);
Asub=measuredMode.(As[V,gmix].measuredMode);
Wig=1/(2)*((x/Vsub)^2*Asub - Asub/Vsub+2)*Exp[-1/2*(x^2/Vsub)]/Sqrt[2*Pi*Vsub];
Return[Wig];];

(**
Wigner Function in the modes "remainingModes" that is obtained when postselecting on a homodyne measurement outcome in the "measuredMode" with particular result "outcome". Note that if we have a total of m modes, V is a 2m x 2m matrix, measuredMode is a 2m-dimensional vector, outcome is a real number,and remaining modes is a list of N ≤ m-1 2m-dimesnional vectors (these should be orthogonal to the measuredMode), and x is and 2N-dimensional vector, organised as (x1,x2,...,xN,p1,...,pN)  
**)
WsPostSelected[x_,V_,gmix_,remainingModes_,measuredMode_,outcome_]:=Module[{β,Vsub,Wig,Asub,m,p,U,int,bla,norm},
β=Join[x,{outcome}];
m=Length[V]/2;
U=Join[Join[remainingModes,Table[J[m].remainingModes[[k]],{k,Length[remainingModes]}]],{measuredMode}];
Vsub=Table[U[[k]].(V.U[[l]]),{k,Length[U]},{l,Length[U]}];
Asub=Table[U[[k]].(As[V,gmix].U[[l]]),{k,Length[U]},{l,Length[U]}];
norm=MeasurementProbability[outcome,V,gmix,measuredMode];
Wig=1/(2)*(β.Inverse[Vsub].Asub.Inverse[Vsub].β - Tr[Inverse[Vsub].Asub]+2)*W0[β,Vsub];
Return[Wig/norm];
];

PrPostSelected[x1_,V_,gmix_,remainingModes_,measuredMode_,outcome_,condquadratureC_,condquadrature_]:=Block[{Wig,prob,norma},
Wig=WsPostSelected[x1,V,gmix,remainingModes,measuredMode,outcome];
prob=Integrate[Wig,{x1[[condquadratureC]],-Infinity,Infinity}];
(*norma=Integrate[prob,{x[[condquadrature]],-Infinity,Infinity}];*)
Return[prob];
]

For the computation of the witness on photon subtracted states

(*Conditional Fisher Information using displacements as generators *)
CondF[{x_,y_},V_,gmix_,remainingModes_,measuredMode_,outcome_,condquadratureC_,condquadrature_,θ_]:=Block[{prob,derp,measurprob,Qf},
prob=PrPostSelected[{x-θ,y},V,gmix,remainingModes,measuredMode,outcome,condquadratureC,condquadrature];
derp=D[prob,θ];
measurprob=MeasurementProbability[outcome,V,gmix,measuredMode];
Qf=NIntegrate[((1/prob) (derp^2) measurprob)/.{θ->0.},{x,-100,100},{outcome,-100,100}];
Return[Qf];
]
(*Conditional variance*)
CondQVar[x1_,V_,gmix_,remainingModes_,measuredMode_,outcome_,condquadratureC_,condquadrature_]:=Block[{prob,derp,measurprob,Qvar,var},
prob=PrPostSelected[x1,V,gmix,remainingModes,measuredMode,outcome,condquadratureC,condquadrature];
measurprob=MeasurementProbability[outcome,V,gmix,measuredMode];
var=Integrate[(x1[[condquadrature]]^2 prob),{x1[[condquadrature]],-Infinity,Infinity}]-(((Integrate[(x1[[condquadrature]] prob),{x1[[condquadrature]],-Infinity,Infinity}]))^2);
Qvar=NIntegrate[var measurprob,{outcome,-30,30}];
Return[Qvar];
]

For the computation using binned data

PrPostSelectedBinned[x1_,V_,gmix_,remainingModes_,measuredMode_,outcome_,li_,ls_,condquadratureC_,condquadrature_]:=Module[{β,Vsub,Wig,Asub,m,p,U,int,bla,norm},
β=Join[x1,{outcome}];
m=Length[V]/2;
U=Join[remainingModes,{measuredMode}];
Vsub=Table[U[[k]].(V.U[[l]]),{k,Length[U]},{l,Length[U]}];
Asub=Table[U[[k]].(As[V,gmix].U[[l]]),{k,Length[U]},{l,Length[U]}];
norm=Integrate[MeasurementProbability[outcome,V,gmix,measuredMode],{outcome,li,ls}];
Wig=Integrate[1/(2)*(β.Inverse[Vsub].Asub.Inverse[Vsub].β - Tr[Inverse[Vsub].Asub]+2)*W0[β,Vsub],{outcome,li,ls}];
Return[Wig/norm];
];
CondFBinned[{x_,y_},V_,gmix_,remainingModes_,measuredMode_,outcome_,intervals_,condquadratureC_,condquadrature_,θ_]:=Block[{prob,derp,measurprob,Qf},
prob=PrPostSelectedBinned[{x},V,gmix,remainingModes,measuredMode,outcome,#[[1]],#[[2]],condquadratureC,condquadrature]&/@intervals;
prob=prob/.{x->(x-θ)};
derp=D[#,θ]&/@prob;
prob=prob/.{θ->0};
derp=derp/.{θ->0};
measurprob=NIntegrate[MeasurementProbability[outcome,V,gmix,measuredMode],{outcome,#[[1]],#[[2]]}]&/@intervals;
Qf=Table[NIntegrate[((1/((prob[[i]])))((( derp[[i]]))^2) measurprob[[i]]),{x,-10,10}],{i,1,Length[prob]}];
Return[{Total[Qf],measurprob,Qf/measurprob}];
]
CondQVarBinned[{x_,y_},V_,gmix_,remainingModes_,measuredMode_,outcome_,intervals_,condquadratureC_,condquadrature_]:=Block[{prob,derp,measurprob,Qvar,var},
prob=PrPostSelectedBinned[{x},V,gmix,remainingModes,measuredMode,outcome,#[[1]],#[[2]],condquadratureC,condquadrature]&/@intervals;
measurprob=Integrate[MeasurementProbability[outcome,V,gmix,measuredMode],{outcome,#[[1]],#[[2]]}]&/@intervals;
var=(NIntegrate[((x^2) #),{x,-20,20}]-(((NIntegrate[(x #),{x,-20,20}]))^2))&/@prob;
Qvar=Total[Table[var[[i]] measurprob[[i]],{i,1,Length[var]}]];
Return[Qvar];
]

For the computation of the witness on Gaussian states

GaussPrPostSelect[x_,V_,remainingModes_,measuredMode_,outcome_]:=Block[{β,Vsub,Wig,Asub,m,p,U,int,bla,norm},
β=Join[x,{outcome}];
m=Length[V]/2;
U=Join[remainingModes,{measuredMode}];
Vsub=Table[U[[k]].(V.U[[l]]),{k,Length[U]},{l,Length[U]}];
norm=GaussMeasProb[outcome,V,measuredMode];
Wig=W0[β,Vsub];
Return[Wig/norm];
]
GaussMeasProb[x_,V_,measuredMode_]:=Block[{Vsub,Wig,m,p,U,int,bla},
m=Length[V]/2;
Vsub=measuredMode.(V.measuredMode);
Wig=Exp[-1/2*(x^2/Vsub)]/Sqrt[2*Pi*Vsub];
Return[Wig];
]

GaussCondF[{x_},V_,remainingModes_,measuredMode_,outcome_,θ_]:=Block[{prob,derp,measurprob,Qf},
prob=GaussPrPostSelect[{x-θ},V,remainingModes,measuredMode,outcome];
derp=D[prob,θ];
measurprob=GaussMeasProb[outcome,V,measuredMode];
Qf=NIntegrate[((1/prob) (derp^2) measurprob)/.{θ->0},{x,-30,30},{outcome,-30,30}];
Return[Qf];
]
GaussCondQVar[x1_,V_,remainingModes_,measuredMode_,outcome_]:=Block[{prob,derp,measurprob,Qvar,var},
prob=GaussPrPostSelect[x1,V,remainingModes,measuredMode,outcome];
measurprob=GaussMeasProb[outcome,V,measuredMode];
var=Integrate[(x^2 prob),{x,-Infinity,Infinity}]-(((Integrate[(x prob),{x,-Infinity,Infinity}]))^2);
Qvar=NIntegrate[var measurprob,{outcome,-100,100}];
Return[Qvar];
]

For computation of the conditional enntropies

CondEntropy[{x_,y_},V_,gmix_,remainingModes_,measuredMode_,outcome_,condquadratureC_,condquadrature_]:=Block[{prob,condh,measurprob},
prob=PrPostSelected[{x,y},V,gmix,remainingModes,measuredMode,outcome,condquadratureC,condquadrature];
measurprob=MeasurementProbability[outcome,V,gmix,measuredMode];
condh=-NIntegrate[( measurprob prob Log[prob]),{x,-20,20},{outcome,-20,20}];
Return[condh];
]
(*linearized witness, comparable to the Reid's one -equal for Gaussian states-*)
Sh[hX_,hP_]:=Max[0,E^(1-2 hX) 2 Pi-E^(2hP-1)/(2 Pi)]

Useful tools

(*Beam splitter matrix*)
M[τ_]:=(({{(({{Sqrt[1-τ], Sqrt[τ]}, {-Sqrt[τ], Sqrt[1-τ]}})), ConstantArray[0,{2,2}]}, {ConstantArray[0,{2,2}], (({{Sqrt[1-τ], Sqrt[τ]}, {-Sqrt[τ], Sqrt[1-τ]}}))}}))//ArrayFlatten
(*Action of the BS in a covariance matrix V1*)
VBS[τ_,V1_]:=Transpose[M[τ]].(V1.M[τ])

AddLosses[Vi_,η_]:=(1-η)Vi+η IdentityMatrix[Length[Vi]]

Gaussian Entangled states (Fig. 3):

(*Here we calculate the Quantum Fisher Information for the two modes squeezed state (i.e two modes initially squeezed in opposed quadratures, then mixed in a Beam splitter)*)
GaussQFsq=GaussCondF[{x},VBS[0.5,DiagonalMatrix[{1/#,#,#,1/#}]],{{1,0,0,0}},{0,1,0,0},q,θ]&/@Table[s,{s,1.01,4,0.2}];

(*Here the quantum conditional variance (calc the variance of quadrature x_B after measuring quadrature x_A, it is equivalent to p_a given p_b in this case)*)
GaussCVarsq=GaussCondQVar[{x},VBS[0.5,DiagonalMatrix[{1/#,#,#,1/#}]],{{0,0,1,0}},{0,0,0,1},q]&/@Table[s,{s,1.01,4,0.2}];

(**Arranging the data**)
squeezing=Table[s,{s,1.01,4,0.2}];
dataGaussQFsq=Transpose[{squeezing,GaussQFsq}];
dataGaussCVarsq=Transpose[{squeezing,GaussCVarsq}];

dataGaussQFsqdB={10 Log10[#[[1]]],#[[2]]}&/@dataGaussQFsq;
dataGaussCVarsqdB={10 Log10[#[[1]]],#[[2]]}&/@dataGaussCVarsq;
dataGaussSteeringsq=Transpose[{10 Log10[squeezing],Max[0,#]&/@((#[[2]]&/@dataGaussQFsq)-(#[[2]]&/@dataGaussCVarsq))}];
Export["./steering_metr_ineq_Gaussian_vs_sq.dat",dataGaussSteeringsq];

dataGaussSteeringsq=Import["./steering_metr_ineq_Gaussian_vs_sq.dat"]

(*This line generates Fig.3*)
ListLinePlot[{dataGaussQFsqdB,dataGaussCVarsqdB,dataGaussSteeringsq},PlotRange->{0,2.5},PlotLegends->Placed[{Style[ToExpression[RowBox[{SuperscriptBox["F","B|A"],"[",RowBox[{"𝒜",",",OverscriptBox["p","^"]}],"]"}]],FontFamily->"Times",FontSize->16],Style[4 ToExpression[RowBox[{SuperscriptBox["Var","B|A"],"[",RowBox[{"𝒜",",",OverscriptBox["p","^"]}],"]"}]],FontFamily->"Times",FontSize->16],Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{"𝒜"}],"]"}]],FontFamily->"Times",FontSize->16]},{{0.3,0.8},Center}],AxesLabel->{"s(dB)"},AxesStyle->{Directive[Black,19],Directive[Black,19]},BaseStyle->{Large,FontFamily->"Times"},PlotStyle->{Continue,Dashed,Dotted}]

Considering Loses (at a squeezing of around 3dB) (Fig. 4)

gaussQFLosses=GaussCondF[{x},AddLosses[VBS[0.5,DiagonalMatrix[{1/2,2,2,1/2}]],#],{{1,0,0,0}},{0,1,0,0},q,θ]&/@Table[s,{s,0,1,0.02}];
gaussCVarLosses=GaussCondQVar[{x},AddLosses[VBS[0.5,DiagonalMatrix[{1/2,2,2,1/2}]],#],{{0,0,1,0}},{0,0,0,1},q]&/@Table[s,{s,0,1,0.02}];

lossesAdd=Table[s,{s,0,1,0.02}];
dataGaussQFLosses=Transpose[{lossesAdd,gaussQFLosses}];
Export["./QFI_Gaussian_vs_Losses3dB.dat",dataGaussQFLosses];
dataGaussCVarLosses=Transpose[{lossesAdd,gaussCVarLosses}];
Export["./CVar_Gaussian_vs_Losses3dB.dat",dataGaussCVarLosses];
dataGaussSteeringLosses=Transpose[{lossesAdd,Max[0,#]&/@(gaussQFLosses-gaussCVarLosses)}];
Export["./steering_metr_ineq_Gaussian_vs_losses_sq_3dB.dat",dataGaussSteeringLosses];

(*This line generates Fig.4 *)
ListLinePlot[{dataGaussSteeringLosses,dataGaussQFLosses,dataGaussCVarLosses},PlotRange->All,Frame->True,FrameTicks->{{{0.0,0.4,0.8,1.2},Automatic},{{0.0,0.2,0.4,0.6,0.8},Automatic}},PlotLegends->Placed[{Style[ToExpression[RowBox[{SuperscriptBox["F","B|A"],"[",RowBox[{SubscriptBox["𝒜","G"],",",OverscriptBox["p/2","^"]}],"]"}]],FontFamily->"Times",FontSize->16],Style[4 ToExpression[RowBox[{SuperscriptBox["Var","B|A"],"[",RowBox[{SubscriptBox["𝒜","G"],",",OverscriptBox["p/2","^"]}],"]"}]],FontFamily->"Times",FontSize->16],Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{SubscriptBox["𝒜","G"]}],"]"}]],FontFamily->"Times",FontSize->16]},{{0.67,0.32}}],FrameLabel->{"η"},FrameStyle->{Directive[Black,19],Directive[Black,19]},BaseStyle->{Large,FontFamily->"Times"},PlotStyle->{Directive[Dashed],Directive[Dotted]}]

Steering Photon subtracted states correlated basis (For Figure 5 and Fig.6)

Steering from the complementary mode

qfSqueezPhtSubtrCompl=Table[CondF[{x,y},VBS[0.5,DiagonalMatrix[{s1,1/s1,1/s1,s1}]],{{1,{1,0,0,0}}},{{1,0,0,0}},Normalize[{0,1,0,0}],q,2,1,θ],{s1,1.01,4,0.2}];
cVarSqueezPhtSubtrCompl=Table[CondQVar[{x,y},VBS[0.5,DiagonalMatrix[{s1,1/s1,1/s1,s1}]],{{1,{1,0,0,0}}},{{0,0,1,0}},{0,0,0,1},q,2,1],{s1,1.01,4,0.2}];

qfdataPhtSubtrCompl=Transpose[{10Log10[squeezing],qfSqueezPhtSubtrCompl}];
Export["./QFI_SubtractedPh_State_Corr_basis_vs_sq.dat",qfdataPhtSubtrCompl];
qVardataPhtSubtrCompl=Transpose[{10Log10[squeezing],cVarSqueezPhtSubtrCompl}];
Export["./CVar_SubtractedPh_State_Corr_basis_vs_sq.dat",qVardataPhtSubtrCompl];
steeringMetrCompl=Transpose[{10Log10[squeezing],Max[#,0]&/@(qfSqueezPhtSubtrCompl-cVarSqueezPhtSubtrCompl)}];
Export["./steering_metr_ineq_Ph_subtracted_corr_basis_vs_sq.dat",steeringPhtSubtrCompl]

steeringMetrCompl=Import["./steering_metr_ineq_Ph_subtracted_corr_basis_vs_sq.dat"];

Steering from the photon subtracted mode:

qfSqueezPhtSubtrPSmode=Table[CondF[{x,y},VBS[0.5,DiagonalMatrix[{s1,1/s1,1/s1,s1}]],{{1,{0,1,0,0}}},{{1,0,0,0}},Normalize[{0,1,0,0}],q,2,1,θ],{s1,1.01,4,0.2}];
cVarSqueezPhtSubtrPSmode=Table[CondQVar[{x,y},VBS[0.5,DiagonalMatrix[{s1,1/s1,1/s1,s1}]],{{1,{0,1,0,0}}},{{0,0,1,0}},{0,0,0,1},q,2,1],{s1,1.01,4,0.2}];

squeezing=Table[s1,{s1,1.01,4,0.2}];
qfdataPhtSubtrPSmode=Transpose[{10 Log10[squeezing],qfSqueezPhtSubtrPSmode}];
Export["./QFI_SubtractedPh_State_Corr_basis_vs_sq_PS_Mode.dat",qfdataPhtSubtrPSmode];
qVardataPhtSubtrPSmode=Transpose[{10 Log10[squeezing],cVarSqueezPhtSubtrPSmode}];
Export["./CVar_SubtractedPh_State_Corr_basis_vs_sq_PS_Mode.dat",qVardataPhtSubtrPSmode];
steeringMetrPSmodedB=Transpose[{10 Log10[squeezing],Max[#,0]&/@(qfSqueezPhtSubtrPSmode-cVarSqueezPhtSubtrPSmode)}];
Export["./steering_metr_ineq_Ph_subtracted_ps_mode_corr_basis_vs_sq.dat",steeringMetrPSmodedB];

steeringMetrPSmodedB=Import["./steering_metr_ineq_Ph_subtracted_ps_mode_corr_basis_vs_sq.dat"];

Computation of the Reid criteria in photon subtracted in the correlated basis

(**Steering from the mode in which the photon is subtracted**)

CorrBasisQVarsq=CondQVar[{x,y},VBS[0.5,DiagonalMatrix[{#,1/#,1/#,#}]],{{1,Normalize[{0,1,0,0}]}},{{1,0,0,0}},{0,1,0,0},q,2,1]&/@Table[s,{s,1.01,4,0.2}];
CorrBasisPVarsq=CondQVar[{x,y},VBS[0.5,DiagonalMatrix[{#,1/#,1/#,#}]],{{1,Normalize[{0,1,0,0}]}},{{0,0,1,0}},{0,0,0,1},q,2,1]&/@Table[s,{s,1.01,4,0.2}];

squeezing=Table[s,{s,1.01,4,0.2}];
dataCorrBasisQVarsq=Transpose[{10 Log10[squeezing],CorrBasisQVarsq}];
dataCorrBasisPVarsq=Transpose[{10 Log10[squeezing],CorrBasisPVarsq}];
dataCorrBasisIneq=Transpose[{10 Log10[squeezing],Max[0,#]&/@((CorrBasisPVarsq^(-1))- CorrBasisQVarsq)}];
Export["./steering_Reid_crit_Ph_subtracted_ps_mode_corr_basis_vs_sq.dat",dataCorrBasisIneq];

(**Steering from the complementary mode**)
CorrBasisQVarsqCompl=CondQVar[{x,y},VBS[0.5,DiagonalMatrix[{1/#,#,#,1/#}]],{{1,Normalize[{1,0,0,0}]}},{{1,0,0,0}},{0,1,0,0},q,2,1]&/@Table[s,{s,1.01,4,0.2}];
CorrBasisPVarsqCompl=CondQVar[{x,y},VBS[0.5,DiagonalMatrix[{1/#,#,#,1/#}]],{{1,Normalize[{1,0,0,0}]}},{{0,0,1,0}},{0,0,0,1},q,2,1]&/@Table[s,{s,1.01,4,0.2}];

(*Linearized version of the Reid criteria*)

dataCorrBasisQVarsqCompl=Transpose[{10 Log10[ squeezing],CorrBasisQVarsqCompl}];
dataCorrBasisPVarsqCompl=Transpose[{10 Log10[ squeezing],CorrBasisPVarsqCompl}];
dataCorrBasisMultVarCompl=Transpose[{10 Log10[ squeezing], ((CorrBasisQVarsqCompl^(-1))-  CorrBasisPVarsqCompl)}];
dataCorrBasisIneqCompl=Transpose[{10 Log10[ squeezing],Max[0,#]&/@((CorrBasisPVarsqCompl^(-1))- CorrBasisQVarsqCompl)}];
Export["./steering_Reid_crit_Ph_subtracted_compl_mode_corr_basis_vs_sq.dat",dataCorrBasisIneq];

Computation of the entropic inequality

From the photon subtracted mode

CorrBasisConHq=Table[CondEntropy[{x,y},VBS[0.5,DiagonalMatrix[{1/s1,s1,s1,1/s1}]],{{1,Normalize[{0,1,0,0}]}},{{1,0,0,0}},Normalize[{0,1,0,0}],q,2,1],{s1,1.01,4,0.2}];
 CorrBasisCondHp=Table[CondEntropy[{x,y},VBS[0.5,DiagonalMatrix[{1/s1,s1,s1,1/s1}]],{{1,Normalize[{0,1,0,0}]}},{{0,0,1,0}},Normalize[{0,0,0,1}],q,2,1],{s1,1.01,4,0.2}];

squeezing=Table[s,{s,1.01,4,0.2}];
dataCorrBasisCondHq=Transpose[{10 Log10[squeezing],CorrBasisConHq}];
dataCorrBasisCondHp=Transpose[{10 Log10[squeezing],CorrBasisCondHp}];
CorrBasisEntropIneq=Transpose[{10 Log10[squeezing],Table[Sh[Abs[CorrBasisConHq[[i]]],Abs[CorrBasisCondHp[[i]]]],{i,Length[CorrBasisConHq]}]}];
Export["./steering_Entropic_crit_Ph_subtracted_PS_mode_corr_basis_vs_sq.dat",CorrBasisEntropIneq];

From the complementary mode

CorrBasisConHqCompl=Table[CondEntropy[{x,y},VBS[0.5,DiagonalMatrix[{1/s1,s1,s1,1/s1}]],{{1,Normalize[{0,1,0,0}]}},{{0,1,0,0}},Normalize[{1,0,0,0}],q,2,1],{s1,1.01,4,0.2}];
 CorrBasisCondHpCompl=Table[CondEntropy[{x,y},VBS[0.5,DiagonalMatrix[{1/s1,s1,s1,1/s1}]],{{1,Normalize[{0,1,0,0}]}},{{0,0,0,1}},Normalize[{0,0,1,0}],q,2,1],{s1,1.01,4,0.2}];

squeezing=Table[s,{s,1.01,4,0.2}];
dataCorrBasisCondHqCompl=Transpose[{10 Log10[squeezing],CorrBasisConHqCompl}];
dataCorrBasisCondHpCompl=Transpose[{10 Log10[squeezing],CorrBasisCondHpCompl}];
(*datagaussCondHppHq=Transpose[{10 Log10[squeezing],gausscondhP+gausscondhX}]*)
CorrBasisEntropIneqCompl=Transpose[{10 Log10[squeezing],Table[Sh[CorrBasisConHqCompl[[i]],CorrBasisCondHpCompl[[i]]],{i,Length[CorrBasisConHqCompl]}]}];
Export["./steering_Entropic_crit_Ph_subtracted_compl_mode_corr_basis_vs_sq.dat",CorrBasisEntropIneq];

(*To generate Fig.5*)

ListLinePlot[{steeringMetrCompl,steeringMetrPSmodedB,CorrBasisEntropIneqCompl,CorrBasisEntropIneq,dataCorrBasisIneqCompl[[;;-2]],dataCorrBasisIneq[[;;-2]]},Frame->True,Ticks->{{{1,2,3,4},Automatic},{{1,2,3,4,5},Automatic}},
PlotLegends->Placed[{Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{SubscriptBox["𝒜","A→B"]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{SubscriptBox["𝒜","B→A"]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubscriptBox["S","H"],"[",RowBox[{SubscriptBox["𝒜","A→B"]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubscriptBox["S","H"],"[",RowBox[{SubscriptBox["𝒜","B→A"]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubsuperscriptBox["S","R","hom"],"[",RowBox[{SubscriptBox["𝒜","A→B"]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubsuperscriptBox["S","R","hom"],"[",RowBox[{SubscriptBox["𝒜","B→A"]}],"]"}]],FontFamily->"Times",FontSize->20]},{{0.45,0.7}}],FrameLabel->{"s(dB)"},FrameStyle->{Directive[Black,19],Directive[Black,19]},BaseStyle->{Large,FontFamily->"Times"},PlotStyle->{Directive[Continue,ColorData[97,"ColorList"][[;;3]][[1]]],Directive[Dashed,ColorData[97,"ColorList"][[;;3]][[1]]],Directive[Continue,ColorData[97,"ColorList"][[;;3]][[2]]],Directive[Dashed,ColorData[97,"ColorList"][[;;3]][[2]]],Directive[Continue,ColorData[97,"ColorList"][[;;4]][[3]]],Directive[Dashed,ColorData[97,"ColorList"][[;;4]][[3]]]},RotateLabel->False,PlotRange->All]

Under the effect of losses (still correlated basis. For Fig.6) (at 5 dB initial squeezing)

Metrological inequality

From the mode where the photon was subtracted

qfLosses=Table[CondF[{x,y},AddLosses[VBS[0.5,DiagonalMatrix[{1/Sqrt[10],Sqrt[10],Sqrt[10],1/Sqrt[10]}]],n1],{{1,Normalize[{0,1,0,0}]}},{{1,0,0,0}},{0,1,0,0},q,2,1,θ],{n1,0,0.5,0.02}]

cVarLosses=Table[CondQVar[{x,y},AddLosses[VBS[0.5,DiagonalMatrix[{1/Sqrt[10],Sqrt[10],Sqrt[10],1/Sqrt[10]}]],n1],{{1,Normalize[{0,1,0,0}]}},{{0,0,1,0}},{0,0,0,1},q,2,1],{n1,0,0.5,0.02}]

lossesAdd=Table[s,{s,0,0.5,0.02}];
qfdataLosses=Transpose[{lossesAdd,qfLosses}];
qVardataLosses=Transpose[{lossesAdd, cVarLosses}];
steeringMetrIneqLossesPSmode=Transpose[{lossesAdd,Max[#,0]&/@(qfLosses- cVarLosses)}]
Export["./steering_Metr_Ineq_losses_Corr_basis_PS_mode.dat",steeringMetrIneqLossesPSmode]

In the complementary direction

qfLossesCompl=Table[CondF[{x,y},AddLosses[VBS[0.5,DiagonalMatrix[{ 1/Sqrt[10],Sqrt[10],Sqrt[10],1/Sqrt[10]}]],n1],{{1,Normalize[{1,0,0,0}]}},{{1,0,0,0}},{0,1,0,0},q,2,1,θ],{n1,0,0.5,0.02}]
cVarLossesCompl=Table[CondQVar[{x,y},AddLosses[VBS[0.5,DiagonalMatrix[{1/Sqrt[10],Sqrt[10],Sqrt[10],1/Sqrt[10]}]],n1],{{1,Normalize[{1,0,0,0}]}},{{0,0,1,0}},{0,0,0,1},q,2,1],{n1,0,0.5,0.02}]

lossesAdd=Table[s,{s,0,0.5,0.02}];
qfdataLossesCompl=Transpose[{lossesAdd,qfLossesCompl}];
qVardataLossesCompl=Transpose[{lossesAdd, cVarLossesCompl}];
steeringMetrIneqLossesCompl=Transpose[{lossesAdd,Max[#,0]&/@(qfLossesCompl- cVarLossesCompl)}]
Export["./steering_Metr_Ineq_losses_Corr_basis_Compl_mode.dat",steeringMetrIneqLossesCompl]

Using the entropic criteria

(*5dB*)
CorrBasisConHqLoss=Table[CondEntropy[{x,y},AddLosses[VBS[0.5,DiagonalMatrix[{ 1/Sqrt[10],Sqrt[10],Sqrt[10],1/Sqrt[10]}]],n1],{{1,Normalize[{0,1,0,0}]}},{{1,0,0,0}},{0,1,0,0},q,2,1],{n1,0,0.15,0.02}]
CorrBasisConHpLoss=Table[CondEntropy[{x,y},AddLosses[VBS[0.5,DiagonalMatrix[{  1/Sqrt[10],Sqrt[10],Sqrt[10],1/Sqrt[10]}]],n1],{{1,Normalize[{0,1,0,0}]}},{{0,0,1,0}},{0,0,0,1},q,2,1],{n1,0,0.15,0.02}]

losses=Table[s,{s,0,0.15,0.02}];
dataCorrBasisCondHqlosses=Transpose[{losses,CorrBasisConHqLoss}];
dataCorrBasisCondHplosses=Transpose[{losses,CorrBasisConHpLoss}];
CorrBasisEntropIneqLosses=Transpose[{losses,Table[Sh[Abs[CorrBasisConHqLoss[[i]]],Abs[CorrBasisConHpLoss[[i]]]],{i,Length[CorrBasisConHqLoss]}]}];
Export["steering_Entr_criteria_losses_Corr_Basis_Ps_mode.dat",CorrBasisEntropIneqLosses]

In the complementary direction

CorrBasisConHqLossCompl=Table[CondEntropy[{x,y},AddLosses[VBS[0.5,DiagonalMatrix[{ 1/Sqrt[10],Sqrt[10],Sqrt[10],1/Sqrt[10]}]],n1],{{1,Normalize[{0,1,0,0}]}},{{0,1,0,0}},{1,0,0,0},q,2,1],{n1,0,0.55,0.05}]
CorrBasisConHpLossCompl=Table[CondEntropy[{x,y},AddLosses[VBS[0.5,DiagonalMatrix[{  1/Sqrt[10],Sqrt[10],Sqrt[10],1/Sqrt[10]}]],n1],{{1,Normalize[{0,1,0,0}]}},{{0,0,0,1}},{0,0,1,0},q,2,1],{n1,0,0.55,0.05}]

losses=Table[s,{s,0,0.55,0.05}];
dataCorrBasisCondHqlossesCompl=Transpose[{losses,CorrBasisConHqLossCompl}];
dataCorrBasisCondHplossesCompl=Transpose[{losses,CorrBasisConHpLossCompl}];
CorrBasisEntropIneqLossesCompl=Transpose[{losses,Table[Sh[Abs[CorrBasisConHqLossCompl[[i]]],Abs[CorrBasisConHpLossCompl[[i]]]],{i,Length[CorrBasisConHqLossCompl]}]}];
Export["steering_Entr_criteria_losses_Corr_Basis_Compl_mode.dat",CorrBasisEntropIneqLosses]

Code to generate Fig.6

ListLinePlot[{steeringMetrIneqLossesCompl,steeringMetrIneqLossesPSmode,CorrBasisEntropIneqLossesCompl,CorrBasisEntropIneqLosses},Frame->True,FrameTicks->{{{0,1,2,3},None},{Automatic,None}},
PlotLegends->Placed[{Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{SubscriptBox["𝒜","A→B"]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{SubscriptBox["𝒜","B→A"]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubscriptBox["S","H"],"[",RowBox[{SubscriptBox["𝒜","A→B"]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubscriptBox["S","H"],"[",RowBox[{SubscriptBox["𝒜","B→A"]}],"]"}]],FontFamily->"Times",FontSize->20]},{{0.3,0.7}}],FrameLabel->{Style["η",FontFamily->"Times",FontSize->19]},FrameStyle->{Directive[Black,19,FontFamily->"Times"],Directive[Black,19,FontFamily->"Times"]},BaseStyle->{Large,FontFamily->"Times"},PlotStyle->{Directive[Continue,ColorData[97,"ColorList"][[;;2]][[1]]],Directive[Dashed,ColorData[97,"ColorList"][[;;2]][[1]]],Directive[Continue,ColorData[97,"ColorList"][[;;2]][[2]]],Directive[Dashed,ColorData[97,"ColorList"][[;;2]][[2]]]},RotateLabel->False]

Steering in the uncorrelated basis (delocalized photon subtraction) Fig. 7,8 and 9

Using the metrological criteria

qfSqueezUncModes=Table[CondF[{x,y},DiagonalMatrix[{1/s1,s1,s1,1/s1}],{{1,Normalize[{1,1,0,0}]}},{{1,0,0,0}},Normalize[{0,1,0,0}],q,2,1,θ],{s1,1.01,4,0.2}]

cVarSqueezUncModes=Table[CondQVar[{x,y},DiagonalMatrix[{1/s1,s1,s1,1/s1}],{{1,Normalize[{1,1,0,0}]}},{{0,0,1,0}},Normalize[{0,0,0,1}],q,2,1],{s1,1.01,4,0.2}]

squeezing=Table[s,{s,1.01,4,0.2}];
qfUncModesdata=Transpose[{10Log10[squeezing],qfSqueezUncModes}];
Export["./QFI_SubtractedPh_State_Unc_basis_vs_sq.dat",qfUncModesdata];
qVarUncModesdata=Transpose[{10Log10[squeezing],cVarSqueezUncModes}];
Export["./QVar_SubtractedPh_State_Unc_basis_vs_sq.dat",qVarUncModesdata];
steeringUncModesdBMetrIneq=Transpose[{10Log10[squeezing],Max[#,0]&/@(qfSqueezUncModes-cVarSqueezUncModes)}]
Export["./steering_Metr_ineq_Uncorrelated_basis_vs_sq.dat",steeringUncModesdBMetrIneq]

(*This generates a plot with the information above (not in the paper)*)
ListLinePlot[{steeringUncModesdBMetrIneq,qfUncModesdata,qVarUncModesdata},Frame->True,FrameTicks->Automatic,PlotRange->{0,8.9},FrameLabel->{"s(dB)"},FrameStyle->{Directive[Black,19],Directive[Black,19]},PlotLegends->Placed[{Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{"𝒜"}],"]"}]],FontFamily->"Times",FontSize->19],Style[ToExpression[RowBox[{SubsuperscriptBox["F","hom","B|A"],"[",RowBox[{"𝒜",",",OverscriptBox["p/2","^"]}],"]"}]],FontFamily->"Times",FontSize->19],Style[ToExpression[RowBox[{SubsuperscriptBox["Var","hom","B|A"],"[",RowBox[{"𝒜",",",OverscriptBox["p","^"]}],"]"}]],FontFamily->"Times",FontSize->19]},{{0.3,0.75}}],BaseStyle->{Large,FontFamily->"Times"},PlotLabels->Placed["(b)",{{5,1.6}}],PlotStyle->{Directive[Continue],Directive[Dashed],Directive[Dotted]}]

Using Reid’s criteria

UncBasisQVarsq=CondQVar[{x,y},DiagonalMatrix[{1/#,#,#,1/#}],{{1,Normalize[{1,1,0,0}]}},{{1,0,0,0}},{0,1,0,0},q,2,1]&/@Table[s,{s,1.01,4,0.2}];
UncBasisPVarsq=CondQVar[{x,y},DiagonalMatrix[{1/#,#,#,1/#}],{{1,Normalize[{1,1,0,0}]}},{{0,0,1,0}},{0,0,0,1},q,2,1]&/@Table[s,{s,1.01,4,0.2}];

dataUncBasisQVarsq=Transpose[{10 Log10[squeezing],(UncBasisQVarsq^(-1))}];
Export["Var_Q_for_Reid_Uncorrelated_basis_vs_sq.dat",dataUncBasisQVarsq]
dataUncBasisPVarsq=Transpose[{10 Log10[squeezing],UncBasisPVarsq}];
Export["Var_P_for_Reid_Uncorrelated_basis_vs_sq.dat",dataUncBasisPVarsq]
dataUncBasisIneqReid=Transpose[{10 Log10[squeezing],Max[0,#]&/@( (UncBasisQVarsq^(-1))-  UncBasisPVarsq)}];
Export["steering_Reid_crit_Uncorrelated_basis_vs_sq.dat",dataUncBasisIneqReid]

ListLinePlot[{dataUncBasisIneqReid,dataUncBasisQVarsq,dataUncBasisPVarsq},FrameLabel->{"s(dB)"},Frame->True,Axes->False,FrameTicks->{{{0,1.5,3.0,4.5},Automatic},{{0,1,2,3,4,5},Automatic}},PlotLegends->Placed[{Style[ToExpression[RowBox[{SubscriptBox["S","R"],"[",RowBox[{"𝒜"}],"]"}]],FontFamily->"Times",FontSize->16],Style[4"/"ToExpression[RowBox[{SubsuperscriptBox["Var","Q","B|A"],"[",RowBox[{"𝒜",",",OverscriptBox["1/2 q","^"]}],"]"}]],FontFamily->"Times",FontSize->16],Style[4 ToExpression[RowBox[{SubsuperscriptBox["Var","Q","B|A"],"[",RowBox[{"𝒜",",",OverscriptBox["p/2","^"]}],"]"}]],FontFamily->"Times",FontSize->16]},{{0.35,0.79},Center}],FrameStyle->{Directive[Black,19],Directive[Black,19]},PlotLabels->Placed[Style["(a)",Large],{{5,1},Center}],BaseStyle->{Large,FontFamily->"Times"},PlotStyle->{Directive[Black],Directive[Dashed,Black],Directive[Dotted,Black]},PlotRange->{Automatic,{-0.01,5.1}}]

Using the entropic inequalities

UncBasisConHq=Table[CondEntropy[{x,y},DiagonalMatrix[{1/s1,s1,s1,1/s1}],{{1,Normalize[{1,1,0,0}]}},{{1,0,0,0}},Normalize[{0,1,0,0}],q,2,1],{s1,1.01,4,0.2}];
UncBasisCondHp=Table[CondEntropy[{x,y},DiagonalMatrix[{1/s1,s1,s1,1/s1}],{{1,Normalize[{1,1,0,0}]}},{{0,0,1,0}},Normalize[{0,0,0,1}],q,2,1],{s1,1.01,4,0.2}];

squeezing=Table[s,{s,1.01,4,0.2}];
dataUncBasisCondHq=Transpose[{10 Log10[squeezing],UncBasisConHq}];
dataUncBasisCondHp=Transpose[{10 Log10[squeezing],UncBasisCondHp}];
UncBasisEntropIneq=Transpose[{10 Log10[squeezing],Table[Sh[UncBasisConHq[[i]],UncBasisCondHp[[i]]],{i,Length[UncBasisConHq]}]}];
Export["steering_Entr_criteria_Uncorrelated_basis_vs_sq.dat",UncBasisEntropIneq]

Generating Fig.7

ListLinePlot[{steeringUncModesdBMetrIneq,UncBasisEntropIneq,dataUncBasisIneqReid},Frame->True,FrameTicks->{{{0,1,2},None},{Automatic,None}},PlotRange->All,FrameLabel->{"s(dB)"},FrameStyle->{Directive[Black,19,FontFamily->"Times"],Directive[Black,19,FontFamily->"Times"]},PlotLegends->Placed[{Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{SubscriptBox["𝒜",RowBox[{"θ","≡",FractionBox["π","4"]}]]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubscriptBox["S","H"],"[",RowBox[{SubscriptBox["𝒜",RowBox[{"θ","≡",FractionBox["π","4"]}]]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubsuperscriptBox["S","R","hom"],"[",RowBox[{SubscriptBox["𝒜",RowBox[{"θ","≡",FractionBox["π","4"]}]]}],"]"}]],FontFamily->"Times",FontSize->20]},{{0.3,0.75}}],BaseStyle->{Large,FontFamily->"Times"},PlotStyle->{Directive[Continue],Directive[Continue],Directive[Continue]}]

Considering the effect of losses Fig. 8

Using the metrological criteria (run over many values of the squeezing)

qfUncModesLosses=Table[Table[CondF[{x,y},AddLosses[DiagonalMatrix[{ 1/s ,s, s, 1/s}],n1],{{1,Normalize[{1,1,0,0}]}},{{1,0,0,0}},{0,1,0,0},q,2,1,θ],{n1,0,0.1,0.0025}],{s,1.01,7.51,0.5}]

cVarUncModesLosses=Table[Table[CondQVar[{x,y},AddLosses[DiagonalMatrix[{ 1/s ,s, s, 1/s}],n1],{{1,Normalize[{1,1,0,0}]}},{{0,0,1,0}},{0,0,0,1},q,2,1],{n1,0,0.1,0.0025}],{s,1.01,7.6,0.5}]

lossesAdd3=Table[s,{s,0,0.1,0.0025}];
qfUncModesdataLosses=Transpose[{lossesAdd3,#}]&/@qfUncModesLosses;
Export["./QFI_SubtractedPh_State_Unc_basis_vs_losses.dat",qfUncModesdataLosses]
qVarUncModesdataLosses=Transpose[{lossesAdd3,#}]&/@cVarUncModesLosses;
Export["./QVar_SubtractedPh_State_Unc_basis_vs_losses.dat",qVarUncModesdataLosses]
steeringUncModesLossesMetrIneq=Transpose[{lossesAdd3,#}]&/@Map[Max[#,0]&,(qfUncModesLosses- cVarUncModesLosses),{2}]

Generating Fig. 9

Show[Table[ListLinePlot[steeringUncModesLossesMetrIneq[[i]],Frame->True,PlotRange->All,PlotStyle->{RGBColor[{i/Length[steeringUncModesLossesMetrIneq[[1;;-2]]],0,0}]},FrameStyle->{Directive[Black,19,Times],Directive[Black,19,Times]},FrameLabel->{Style["η",FontSize->19],Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{"𝒜"}],"]"}]],FontFamily->"Times",FontSize->19]},RotateLabel->False,FrameTicks->{{{0,2,4},Automatic},{{0,0.03,0.06,0.09},Automatic}},PlotLabels->Placed[{Style["s(dB)",FontSize->19,FontFamily->"Times"],Style["",ImageSizeMultipliers->{0.75, 0.75}]},{{0.065,4.},{0.08,0.14}}]],{i,1,Length[steeringUncModesLossesMetrIneq]-1}],PlotRange->{Automatic,{0,5}},BaseStyle->{Large,FontFamily->"Times"}]

Using the entropic inequalities (run over many values of the squeezing)

UncBasisConHqLoss=Table[Table[CondEntropy[{x,y},AddLosses[DiagonalMatrix[{ 1/s,s,s,1/s}],n1],{{1,Normalize[{1,1,0,0}]}},{{1,0,0,0}},{0,1,0,0},q,2,1],{n1,0,0.1,0.002}],{s,1.01,7.6,0.5}]

UncBasisConHpLoss=Table[Table[CondEntropy[{x,y},AddLosses[DiagonalMatrix[{ 1/s,s,s,1/s}],n1],{{1,Normalize[{1,1,0,0}]}},{{0,0,1,0}},{0,0,0,1},q,2,1],{n1,0,0.1,0.002}],{s,1.01,7.6,0.5}]

losses=Table[s,{s,0,0.1,0.002}];
dataUnCorrBasisCondHqlosses=Transpose[{losses,#}]&/@UncBasisConHqLoss;
dataUnCorrBasisCondHplosses=Transpose[{losses,#}]&/@UncBasisConHpLoss;
UncBasisEntropIneqLosses=Table[Transpose[{losses,Table[Sh[Abs[UncBasisConHqLoss[[i,j]]],Abs[UncBasisConHpLoss[[i,j]]]],{j,Length[UncBasisConHqLoss[[i]]]}]}],{i,Length[UncBasisConHqLoss]}];
Export["./steering_Entr_criteria_Ph_State_Unc_basis_vs_losses.dat",UncBasisEntropIneqLosses]

Show[Table[ListLinePlot[UncBasisEntropIneqLosses[[i]],Frame->True,PlotRange->All,PlotStyle->{RGBColor[{i/Length[UncBasisEntropIneqLosses[[1;;-2]]],0,0}]},FrameStyle->{Directive[Black,19],Directive[Black,19]},FrameLabel->{Style["η",FontSize->24],"(Subscript[S, H])"},FrameTicks->{{{0,0.4,0.8,1.2},None},Automatic},PlotLabels->Placed[{Style["s(dB)",FontSize->19,FontFamily->"Times"],Style["",ImageSizeMultipliers->{0.05, 0.05}]},{{0.065,0.8},{0.08,0.06}}]],{i,1,Length[UncBasisEntropIneqLosses]-1}],PlotRange->All,BaseStyle->{Large,FontFamily->"Times"}]

Generating Fig. 8

(*Code to plot Fig. 8*) (*We pick squeezing of around 4.7 dB*)
ListLinePlot[{steeringUncModesLossesMetrIneq[[5]],UncBasisEntropIneqLosses[[5]]},Frame->True,FrameTicks->{{{0,1,2},None},{Automatic,None}},PlotRange->{{0,0.06},All},FrameLabel->{"η"},FrameStyle->{Directive[Black,19,FontFamily->"Times"],Directive[Black,19,FontFamily->"Times"]},PlotLegends->Placed[{Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{SubscriptBox["𝒜",RowBox[{"θ","≡",FractionBox["π","4"]}]]}],"]"}]],FontFamily->"Times",FontSize->20],Style[ToExpression[RowBox[{SubscriptBox["S","H"],"[",RowBox[{SubscriptBox["𝒜",RowBox[{"θ","≡",FractionBox["π","4"]}]]}],"]"}]],FontFamily->"Times",FontSize->20]},{{0.65,0.75}}],BaseStyle->{Large,FontFamily->"Times"},PlotStyle->{Directive[Continue],Directive[Dashed]}]

Code for Fig.10 (comparison of the metrological steering in different scenarios)

ListLinePlot[{dataGaussSteeringsq,steeringMetrCompl,steeringMetrPSmodedB,steeringUncModesdBMetrIneq},PlotRange->All,Frame->True,FrameTicks->Automatic,PlotLegends->Placed[{Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{SubscriptBox["𝒜","G"]}],"]"}]],FontFamily->"Times",FontSize->19],Style[ ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",SubsuperscriptBox["𝒜",RowBox[{"θ","≡","0"}],RowBox[{"(A","→","B)"}]],"]"}]],FontFamily->"Times",FontSize->19],Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",SubsuperscriptBox["𝒜",RowBox[{"θ","≡","0"}],RowBox[{"(B","→","A)"}]],"]"}]],FontFamily->"Times",FontSize->19],Style[ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{SubscriptBox["𝒜",RowBox[{"θ","≡",FractionBox["π","4"]}]]}],"]"}]],FontFamily->"Times",FontSize->19]},{{0.25,0.69},Center}],AxesLabel->{"s(dB)",""},FrameStyle->{Directive[Black,19],Directive[Black,19]},BaseStyle->{Large,FontFamily->"Times"},PlotStyle->{Directive[Continue],Directive[Dashed],Directive[Dotted],Directive[DotDashed]}]

Binning the sample space (Fig. 11 and 12)

5 bins

CondqfBinned5bins=CondFBinned[{x,y},DiagonalMatrix[{1/#,#,#,1/#}],{{1,{1,1,0,0}}},{{1,0,0,0}},{0,1,0,0},q,{{-10,-1.4},{-1.4,-0.3},{-0.3,0.3},{0.3,1.4},{1.4,10}},2,1,θ]&/@Table[sq,{sq,1.01,4.01,0.5}]
CondqvarBinned5bins=CondQVarBinned[{x,y},DiagonalMatrix[{1/#,#,#,1/#}],{{1,{1,1,0,0}}},{{0,0,1,0}},{0,0,0,1},q,{{-10,-1.4},{-1.4,-0.3},{-0.3,0.3},{0.3,1.4},{1.4,10}},2,1]&/@Table[sq,{sq,1.01,4.01,0.5}]

squeezing=Table[s,{s,1.01,4.01,0.5}];
qfUncModesdataBinned5bins=Transpose[{10 Log10[squeezing],#[[1]]&/@CondqfBinned5bins}];
Export["./QFI_SubtractedPh_State_Unc_basis_vs_sq_Binned_Opt_5bins.dat",qfUncModesdataBinned5bins];
qVarUncModesdataBinned5bins=Transpose[{10 Log10[squeezing],CondqvarBinned5bins}];
Export["./QVar_SubtractedPh_State_Unc_basis_vs_sq_Binned_opt_5bins.dat",qVarUncModesdataBinned5bins];
steeringUncModesBinned5bins=Transpose[{10 Log10[squeezing],Max[#,0]&/@(CondqfBinnedOpt5bins-CondqvarBinned5bins)}];

11 bins

CondqfBinned11bins=CondFBinned[{x,y},DiagonalMatrix[{1/#,#,#,1/#}],{{1,{1,1,0,0}}},{{1,0,0,0}},{0,1,0,0},q,{{-Infinity,-6},{-6,-3.5},{-3.5,-2.5},{-2.5,-1},{-1,-0.3},{-0.3,0.3},{0.3,1},{1,2.5},{2.5,3.5},{3.5,6},{6,Infinity}},2,1,θ]&/@Table[sq,{sq,1.01,4.01,0.5}]

CondqvarBinned11bins=CondQVarBinned[{x,y},DiagonalMatrix[{1/#,#,#,1/#}],{{1,{1,1,0,0}}},{{0,0,1,0}},{0,0,0,1},q,{{-Infinity,-6},{-6,-3.5},{-3.5,-2.5},{-2.5,-1},{-1,-0.3},{-0.3,0.3},{0.3,1},{1,2.5},{2.5,3.5},{3.5,6},{6,Infinity}},2,1]&/@Table[sq,{sq,1.01,4.01,0.5}]

squeezing=Table[s,{s,1.01,4.01,0.5}];
qfUncModesdataBinned11bins=Transpose[{10 Log10[squeezing],#[[1]]&/@CondqfBinned11bins}];
Export["./QFI_SubtractedPh_State_Unc_basis_vs_sq_Binned_11_bins.dat",qfUncModesdataBinedOpt11];
qVarUncModesdataBinned11bins=Transpose[{10 Log10[squeezing],CondqvarBinned11bins}];
Export["./QVar_SubtractedPh_State_Unc_basis_vs_sq_Binned_11_bins.dat",qVarUncModesdataBinned11];
steeringUncModesBinned11=Transpose[{10 Log10[squeezing],Max[#,0]&/@(CondqfBined11bins-CondqvarBined11bins)}];

13 bins

CondFSqBinned13Bins=CondFBinned[{x,y},DiagonalMatrix[{1/#,#,#,1/#}],{{1,{1,1,0,0}}},{{1,0,0,0}},{0,1,0,0},q,{{-10,-3.5},{-3.5,-2},{-2,-1.5},{-1.5,-1},{-1,-0.5},{-0.5,-0.1},{-0.1,0.1},{0.1,0.5},{0.5,1},{1,1.5},{1.5,2},{2,3.5},{3.5,10}},2,1,θ]&/@Table[sq,{sq,1.01,4.01,0.2}]

CondVarSqBinned13Bins=CondQVarBinned[{x,y},DiagonalMatrix[{1/#,#,#,1/#}],{{1,{1,1,0,0}}},{{0,0,1,0}},{0,0,0,1},q,{{-10,-3.5},{-3.5,-2},{-2,-1.5},{-1.5,-1},{-1,-0.5},{-0.5,-0.1},{-0.1,0.1},{0.1,0.5},{0.5,1},{1,1.5},{1.5,2},{2,3.5},{3.5,10}},2,1]&/@Table[sq,{sq,1.01,4.01,0.2}]

squeezing=Table[s,{s,1.01,4.01,0.2}];
CondqfUncModesdataBinned13bins=Transpose[{10 Log10[squeezing],CondFSqBinned13Bins}];
Export["./QFI_SubtractedPh_State_Unc_basis_vs_sq_Binned_13bins.dat",CondqfUncModesdataBinned13bins];
CondqVarUncModesdataBinned13bins=Transpose[{10 Log10[squeezing],CondVarSqBinned13Bins}];
Export["./QVar_SubtractedPh_State_Unc_basis_vs_sq_Binned_13bins.dat",CondqVarUncModesdataBinned13bins];
steeringUncModesBinned13bins=Transpose[{10 Log10[squeezing],Max[#,0]&/@(CondqfSqBinned13bins-CondqVarSqBinned13bins)}];

Code to generate Fig.11

ListLinePlot[{steeringUncModesCont,steeringUncModesBinned13bins,steeringUncModesBinned11,steeringUncModesBinned5bins},Frame->True,FrameTicks->{{{0.5,1.5,2.5},Automatic},Automatic},PlotLegends->Placed[{Style["Continuous",FontFamily->"Times",FontSize->19],Style["13 bins",FontFamily->"Times",FontSize->19],Style["11 bins",FontFamily->"Times",FontSize->19],Style["5 bins",FontFamily->"Times",FontSize->19]},{{0.3,0.77},Center}],FrameLabel->{"s(dB)",ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{"𝒜"}],"]"}]],FontFamily->"Times",FontSize->19},FrameStyle->{Directive[Black,19],Directive[Black,19]},AxesOrigin->{0,0},BaseStyle->{Large,FontFamily->"Times"},PlotStyle->{Directive[Continue],Directive[Dashed],Directive[Dotted],Directive[DotDashed]},RotateLabel->False]

Considering the effect of losses (Binned)

5 bins

CondqfBinnedLosses5bins=CondFBinned[{x,y},AddLosses[DiagonalMatrix[{ 1/2.5, 2.5, 2.5, 1/2.5}],#],{{1,{1,1,0,0}}},{{1,0,0,0}},{0,1,0,0},q,{{-10,-1.4},{-1.4,-0.3},{-0.3,0.3},{0.3,1.4},{1.4,10}},2,1,θ]&/@Table[sq,{sq,0,0.1,0.0025}]
CondcvarBinnedLosses5bins=CondQVarBinned[{x,y},AddLosses[DiagonalMatrix[{ 1/2.5, 2.5, 2.5, 1/2.5}],#],{{1,{1,1,0,0}}},{{0,0,1,0}},{0,0,0,1},q,{{-10,-1.4},{-1.4,-0.3},{-0.3,0.3},{0.3,1.4},{1.4,10}},2,1]&/@Table[sq,{sq,0,0.1,0.0025}]

losses=Table[sq,{sq,0,0.1,0.0025}];
qfUncModesdataBinnedLosses5bins=Transpose[{losses,#[[1]]&/@CondqfBinnedLosses5bins}];
Export["./QFI_SubtractedPh_State_Unc_basis_vs_sq_Binned_losses_5_bins.dat",qfUncModesdataBinnedLosses5bins];
qVarUncModesdataBinnedLosses5=Transpose[{losses,CondcvarBinnedLosses5bins}];
Export["/Users/lopetegui/Documents/Stage/data/QVar_SubtractedPh_State_Unc_basis_vs_sq_Binned_losses_5_bins.dat",qVarUncModesdataBinnedLosses5bins];
steeringUncModesBinnedLosses5bins=Transpose[{losses,Max[#,0]&/@(qfUncModesdataBinnedLosses5bins-qVarUncModesdataBinnedLosses5bins)}];

11 bins

CondqfBinnedLosses11bins=CondFBinned[{x,y},AddLosses[DiagonalMatrix[{1/2.5,2.5,2.5,1/2.5}],#],{{1,{1,1,0,0}}},{{1,0,0,0}},{0,1,0,0},q,{{-10,-6},{-6,-3.5},{-3.5,-2.5},{-2.5,-1},{-1,-0.3},{-0.3,0.3},{0.3,1},{1,2.5},{2.5,3.5},{3.5,6},{6,10}},2,1,θ]&/@Table[sq,{sq,0,0.1,0.0025}]
CondcvarBinnedLosses11bins=CondQVarBinned[{x,y},AddLosses[DiagonalMatrix[{1/2.5,2.5,2.5,1/2.5}],#],{{1,{1,1,0,0}}},{{0,0,1,0}},{0,0,0,1},q,{{-10,-6},{-6,-3.5},{-3.5,-2.5},{-2.5,-1},{-1,-0.3},{-0.3,0.3},{0.3,1},{1,2.5},{2.5,3.5},{3.5,6},{6,10}},2,1]&/@Table[sq,{sq,0,0.1,0.0025}]

losses=Table[sq,{sq,0,0.1,0.0025}];
qfUncModesdataBinnedLosses11bins=Transpose[{losses,#[[1]]&/@CondqfBinnedLosses11bins}];
Export["./QFI_SubtractedPh_State_Unc_basis_vs_sq_Binned_Losses_11_bins.dat",qfUncModesdataBinnedLosses11bins];
qVarUncModesdataBinnedLosses11bins=Transpose[{losses,CondcvarBinnedLosses11bins}];
Export["./QVar_SubtractedPh_State_Unc_basis_vs_sq_Binned_Losses_11_bins.dat",qVarUncModesdataBinnedLosses11bins];
steeringUncModesBinnedLosses11bins=Transpose[{losses,Max[#,0]&/@(qfBinedLosses11bins-cvarBinnedLosses11bins)}];

13 bins

CondqfBinnedLosses13bins=CondFBinned[{x,y},AddLosses[DiagonalMatrix[{ 1/2.5, 2.5, 2.5, 1/2.5}],#],{{1,{1,1,0,0}}},{{1,0,0,0}},{0,1,0,0},q,{{-10,-3.5},{-3.5,-2},{-2,-1.5},{-1.5,-1},{-1,-0.5},{-0.5,-0.1},{-0.1,0.1},{0.1,0.5},{0.5,1},{1,1.5},{1.5,2},{2,3.5},{3.5,10}},2,1,θ]&/@Table[sq,{sq,0,0.1,0.0025}]
CondcvarBinnedLosses13bins=CondQVarBinned[{x,y},AddLosses[DiagonalMatrix[{ 1/2.5, 2.5, 2.5, 1/2.5}],#],{{1,{1,1,0,0}}},{{0,0,1,0}},{0,0,0,1},q,{{-10,-3.5},{-3.5,-2},{-2,-1.5},{-1.5,-1},{-1,-0.5},{-0.5,-0.1},{-0.1,0.1},{0.1,0.5},{0.5,1},{1,1.5},{1.5,2},{2,3.5},{3.5,10}},2,1]&/@Table[sq,{sq,0,0.1,0.0025}]

losses=Table[sq,{sq,0,0.1,0.0025}];
qfUncModesdataBinnedLosses13bins=Transpose[{losses,#[[1]]&/@CondqfBinnedLosses13bins}];
Export["./QFI_SubtractedPh_State_Unc_basis_vs_sq_Binned_losses_13_bins.dat",qfUncModesdataBinnedLosses13bins];
qVarUncModesdataBinnedLosses13bins=Transpose[{losses,CondcvarBinnedLosses13bins}];
Export["./QVar_SubtractedPh_State_Unc_basis_vs_sq_Binned_losses_13_bins.dat",qVarUncModesdataBinnedLosses13bins];
steeringUncModesBinnedLosses13bins=Transpose[{losses,Max[#,0]&/@(qfUncModesdataBinnedLosses13bins-qVarUncModesdataBinnedOptLosses13)}];

Code to generate Fig.12

ListLinePlot[{steeringUncModesLossesMetrIneq[[4]][[;;-10]],steeringUncModesBinnedLosses13bins[[;;-10]],steeringUncModesBinnedLosses11bins[[;;-10]],steeringUncModesBinnedLosses5bins[[;;-10]]},Frame->True,FrameTicks->Automatic,PlotLegends->Placed[{Style["Continuous",FontFamily->"Times",FontSize->19],Style["13 bins",FontFamily->"Times",FontSize->19],Style["11 bins",FontFamily->"Times",FontSize->19],Style["5 bins",FontFamily->"Times",FontSize->19]},{Right,Center}],FrameStyle->{Directive[Black,19],Directive[Black,19]},FrameLabel->{Style["η",FontFamily->"Times",FontSize->19],ToExpression[RowBox[{SubsuperscriptBox["S","max","hom"],"[",RowBox[{"𝒜"}],"]"}]],FontFamily->"Times",FontSize->19},PlotRange->All,BaseStyle->{Large,FontFamily->"Times"},PlotStyle->{Directive[Continue],Directive[Dashed],Directive[Dotted],Directive[DotDashed]},RotateLabel->False]