This file contains supplementary data for Marina Cuzminschi, Alexei Zubarev and Aurelian Isar, “Extractable quantum work from a two-mode Gaussian state in a noisy channel,” Scientific Reports, 11(1) 2021 24286. https://doi.org/10.1038/s41598-021-03752-4.

Extractable Quantum Work from a Two-Mode Gaussian State In a Noisy Channel (Resonant Case)

Alexei Zubarev

k=1;
m=1; 
(*We set Boltzmann constant and mass to one.*)
λ=0.1;
(*λ - dissipation parameter*)

a[r_,n1_,n2_]=n1 ((Cosh[r])^2)+n2 ((Sinh[r])^2)+(1/2)Cosh[2r];
b[r_,n1_,n2_]=n1 ((Sinh[r])^2)+n2 ((Cosh[r])^2)+(1/2)Cosh[2r];
c[r_,n1_,n2_]=(1/2)(n1+n2+1)Sinh[2r];
σ0[r_,n1_,n2_]=(({{a[r,n1,n2], 0, c[r,n1,n2], 0}, {0, a[r,n1,n2], 0, -c[r,n1,n2]}, {c[r,n1,n2], 0, b[r,n1,n2], 0}, {0, -c[r,n1,n2], 0, b[r,n1,n2]}}));
(*Here we defined the covariance matrix for initial time in terms of squeezing parameter r, average numbers of thermal photons of the first and second mode correspondingly ((Subscript[n, 1]) and (Subscript[n, 2])). The initial state is a squeezed thermal state.*)
nth[T_]=(1/2)(Coth[(1/(2T))]-1); 
(*(Subscript[n, th]) is the average number of thermal photons of the resource state. We consider a resonant case first (The frequencies of the modes are (Subscript[ω, 1])=(Subscript[ω, 2])=1).*)

Nn[T_,Rr_]=nth[T](((Cosh[Rr])^2)+((Sinh[Rr])^2))+((Sinh[Rr])^2);
(*Rr is the squeezing parameter of the noisy channel.*)

MRe[T_,Rr_,ϕ_]=-(2nth[T]+1)Cosh[Rr]Sinh[Rr]Cos[ϕ];
MIm[T_,Rr_,ϕ_]=-(2nth[T]+1)Cosh[Rr]Sinh[Rr]Sin[ϕ];
(*Noisy channel is described in terms of it squeezing parameter Rr, temperature T and the phase ϕ.*)

σ∞[T_,Rr_,ϕ_]=(({{(1/2)+(Nn[T,Rr]+MRe[T,Rr,ϕ]), MIm[T,Rr,ϕ], 0, 0}, {MIm[T,Rr,ϕ], (1/2)+(Nn[T,Rr]-MRe[T,Rr,ϕ]), 0, 0}, {0, 0, (1/2)+(Nn[T,Rr]+MRe[T,Rr,ϕ]), MIm[T,Rr,ϕ]}, {0, 0, MIm[T,Rr,ϕ], (1/2)+(Nn[T,Rr]-MRe[T,Rr,ϕ])}}));
(*Here we define the asymptotic covariance matrix for time t=∞*)
σa∞[T_,Rr_,ϕ_]=(({{σ∞[T,Rr,ϕ][[1,1]], σ∞[T,Rr,ϕ][[1,2]]}, {σ∞[T,Rr,ϕ][[2,1]], σ∞[T,Rr,ϕ][[2,2]]}}));

σb∞[T_,Rr_,ϕ_]=(({{σ∞[T,Rr,ϕ][[3,3]], σ∞[T,Rr,ϕ][[3,4]]}, {σ∞[T,Rr,ϕ][[4,3]], σ∞[T,Rr,ϕ][[4,4]]}}));
σc∞[T_,Rr_,ϕ_]=(({{σ∞[T,Rr,ϕ][[1,3]], σ∞[T,Rr,ϕ][[1,4]]}, {σ∞[T,Rr,ϕ][[2,3]], σ∞[T,Rr,ϕ][[2,4]]}}));
(*Here are given the the asymptotic covariance matrices for mode 1 ((Subscript[σ, a])(∞)), the second mode ((Subscript[σ, b])(∞)), and correlations between them ((Subscript[σ, c])(∞)).*)

σ[t_,r_,n1_,n2_,T_,Rr_,ϕ_]=Exp[-λ t]σ0[r,n1,n2]+(1-Exp[-λ t])σ∞[T,Rr,ϕ];
(*The time evolution of the two-mode state in a noisy channel is presented above.*)

σa[t_,r_,n1_,n2_,T_,Rr_,ϕ_]=(({{σ[t,r,n1,n2,T,Rr,ϕ][[1,1]], σ[t,r,n1,n2,T,Rr,ϕ][[1,2]]}, {σ[t,r,n1,n2,T,Rr,ϕ][[2,1]], σ[t,r,n1,n2,T,Rr,ϕ][[2,2]]}}));
σb[t_,r_,n1_,n2_,T_,Rr_,ϕ_]=(({{σ[t,r,n1,n2,T,Rr,ϕ][[3,3]], σ[t,r,n1,n2,T,Rr,ϕ][[3,4]]}, {σ[t,r,n1,n2,T,Rr,ϕ][[4,3]], σ[t,r,n1,n2,T,Rr,ϕ][[4,4]]}}));
σc[t_,r_,n1_,n2_,T_,Rr_,ϕ_]=(({{σ[t,r,n1,n2,T,Rr,ϕ][[1,3]], σ[t,r,n1,n2,T,Rr,ϕ][[1,4]]}, {σ[t,r,n1,n2,T,Rr,ϕ][[2,3]], σ[t,r,n1,n2,T,Rr,ϕ][[2,4]]}}));
(*The matrix of the (1^st) and (2^nd) modes are defined, together with the matrix of their correlations.*)

Ref[θ_]=Cos[θ]IdentityMatrix[2]-I Sin[θ](({{0, -I}, {I, 0}}));
(*θ is the measurement angle.*)
γπb[θ_,μ_]=Ref[θ] .(({{(μ/2), 0}, {0, (1/(2μ))}})).Transpose[Ref[θ]];  (*(γ^((Subscript[π, b]) ))is the seed of measurement matrix. μ=0 Homodyne detection, μ=1 Heterodyne detection*)

σaπb[t_,r_,n1_,n2_,T_,Rr_,ϕ_,θ_,μ_]=σa[t,r,n1,n2,T,Rr,ϕ]-σc[t,r,n1,n2,T,Rr,ϕ].Inverse[σb[t,r,n1,n2,T,Rr,ϕ]+γπb[θ,μ]].Transpose[σc[t,r,n1,n2,T,Rr,ϕ]];
(*(Subscript[σ, (a^(Subscript[π, b]))]) is the covariance matrix of the state of the first mode reached after the measurement of the (2^nd) mode has been performed. The state changes due to back-reaction.*)

W[t_,r_,n1_,n2_,T_,Rr_,ϕ_,θ_,μ_]=(T/2)Log[((Det[σa[t,r,n1,n2,T,Rr,ϕ]])/(Det[σaπb[t,r,n1,n2,T,Rr,ϕ,θ,μ]]))];
(*The quentum work extracted when the state of the mode 1 reaches the new equilibrium with the noisy channel in terms of Renyi entropy.*)
Δ[t_,r_,n1_,n2_,T_,Rr_,ϕ_]=4Det[σa[t,r,n1,n2,T,Rr,ϕ]]+4Det[σb[t,r,n1,n2,T,Rr,ϕ]]+8Det[σc[t,r,n1,n2,T,Rr,ϕ]];
(*Seralian of the two-mode states.*)
(*(ν^2)-Δν+Det[σ]=0*)
ν1[t_,r_,n1_,n2_,T_,Rr_,ϕ_]=Sqrt[((Δ[t,r,n1,n2,T,Rr,ϕ]-Sqrt[((Δ[t,r,n1,n2,T,Rr,ϕ])^2)-64Det[σ[t,r,n1,n2,T,Rr,ϕ]]])/2)];
ν2[t_,r_,n1_,n2_,T_,Rr_,ϕ_]=Sqrt[((Δ[t,r,n1,n2,T,Rr,ϕ]+Sqrt[((Δ[t,r,n1,n2,T,Rr,ϕ])^2)-64Det[σ[t,r,n1,n2,T,Rr,ϕ]]])/2)];
(*Simplectic eigenvalues of the two-mode states.*)
S[t_,r_,n1_,n2_,T_,Rr_,ϕ_]=((ν1[t,r,n1,n2,T,Rr,ϕ]+1)/2)Log2[((ν1[t,r,n1,n2,T,Rr,ϕ]+1)/2)]-((ν1[t,r,n1,n2,T,Rr,ϕ]-1)/2)Log2[((Abs[ν1[t,r,n1,n2,T,Rr,ϕ]-1])/2)]+((ν2[t,r,n1,n2,T,Rr,ϕ]+1)/2)Log2[((ν2[t,r,n1,n2,T,Rr,ϕ]+1)/2)]-((ν2[t,r,n1,n2,T,Rr,ϕ]-1)/2)Log2[((Abs[ν2[t,r,n1,n2,T,Rr,ϕ]-1])/2)];
(*Von Neumann Entropy of a two-mode state.*)
Werasnou[t_,r_,n1_,n2_,T_,Rr_,ϕ_]=T S[t,r,n1,n2,T,Rr,ϕ];
Wmnou[r_,n1_,n2_,T_,Rr_,ϕ_]:=FindMaximum[{Werasnou[t,r,n1,n2,T,Rr,ϕ],0<20<100},{t,1}][[1]]
(*We take the maximal entropy in time.*)
ηnou[t_,r_,n1_,n2_,T_,Rr_,ϕ_,θ_,μ_]:=((W[t,r,n1,n2,T,Rr,ϕ,θ,μ])/(Wmnou[r,n1,n2,T,Rr,ϕ]));
(*The entropy is obtained*)

Plot[{(1/π)NIntegrate[W[0,1.8,0,0,T,0.2,π/4,θ,0],{θ,0,π}],(1/π)NIntegrate[W[1,1.8,0,0,T,0.2,π/4,θ,0],{θ,0,π}],(1/π)NIntegrate[W[5,1.8,0,0,T,0.2,π/4,θ,0],{θ,0,π}]},{T,0,5},Frame->True,FrameLabel->{"T","W"},FrameStyle->26,ImageSize->1000,PlotStyle->{Red,Green,Blue},PlotLegends->Placed[{Style["t=0",FontSize->26],Style["t=1",FontSize->26],Style["t=5",FontSize->26]},Above]]

Plot[{(1/π)NIntegrate[W[t,1.8,0,0,0.1,0.2,π/4,π/4,0],{θ,0,π}],(1/π)NIntegrate[W[t,1.8,0,0,1,0.2,π/4,π/4,0],{θ,0,π}],(1/π)NIntegrate[W[t,1.8,0,0,2,0.2,π/4,π/4,0],{θ,0,π}]},{t,0,5},Frame->True,FrameLabel->{"t","W"},FrameStyle->36,ImageSize->1000,PlotStyle->{Red,Green,Blue},PlotLegends->Placed[{Style["T=0.1",FontSize->30],Style["T=1",FontSize->30],Style["T=2",FontSize->30]},Above]]

Plot[{W[0,1.8,0,0,0.1,Rr,π/4,π/4,0],W[0,1.8,0,0,1,Rr,π/2,π/4,0],W[0,1.8,0,0,2,Rr,π/4,π/4,0]},{Rr,0,5},Frame->True,FrameLabel->{"R","W"},FrameStyle->22,ImageSize->1000,PlotStyle->{Red,Green,Blue},PlotLegends->Placed[{Style["t=0",FontSize->26],Style["t=1",FontSize->26],Style["t=5",FontSize->26]},Above]]

(*t,r,n1,n2,T,Rr,ϕ,θ,μ*)

Plot[{(1/π)NIntegrate[W[0,r,0,0,0.1,0.2,π/4,θ,0],{θ,0,π}],(1/π)NIntegrate[W[1,r,0,0,0.1,0.2,π/4,θ,0],{θ,0,π}],(1/π)NIntegrate[W[2,r,0,0,0.1,0.2,π/4,θ,0],{θ,0,π}]},{r,0,1.8},Frame->True,FrameLabel->{"r","W"},FrameStyle->36,ImageSize->1000,PlotStyle->{Red,Green,Blue},PlotLegends->Placed[{Style["t=0",FontSize->30],Style["t=1",FontSize->30],Style["t=2",FontSize->30]},Above]]

(*W[t,r,n1,n2,T,Rr,ϕ,θ,μ]*)

Plot3D[(1/π)NIntegrate[W[t,1.8,0,0,T,0.2,π/4,θ,0],{θ,0,π}],{t,0,5},{T,0,5},PlotRange->All,AxesLabel->{"t","T","W"},AxesStyle->36,PlotPoints->75,ImageSize->1000,ViewPoint->{4, -10, 5}]

Plot3D[(1/π)NIntegrate[ηnou[t,1.8,0,0,T,0.2,π/4,θ,0],{θ,0,π}],{t,0.1,5},{T,0.1,5},PlotRange->All,AxesLabel->{"t","T","η"},AxesStyle->36,PlotPoints->75,ImageSize->1000,ViewPoint->{4, -10, 5},PlotStyle->Cyan]

Plot3D[(1/π)NIntegrate[ηnew[t,1.8,0,0,T,0.2,π/4,θ,0],{θ,0,π}],{t,0.0001,5},{T,0.0001,5},PlotRange->All,AxesLabel->{"t","T","η"},AxesStyle->36,PlotPoints->75,ImageSize->1000,ViewPoint->{4, -10, 5},PlotStyle->Cyan]

(*W[t,r,n1,n2,T,Rr,ϕ,θ,μ]*)

Plot3D[W[t,1.8,0,0,3,0.2,π/4,θ,0],{t,0,5},{θ,0,2π},PlotRange->All,AxesLabel->{"t","θ","W"},AxesStyle->36,PlotPoints->75,ImageSize->1000,ViewPoint->{4, -10, 5}]

(*W[t,r,n1,n2,T,Rr,ϕ,θ,μ]*)

Plot3D[(1/π)NIntegrate[W[t,1.8,0,0,3,0.2,ϕ,θ,0],{θ,0,π}],{t,0,5},{ϕ,0,6π},PlotRange->All,AxesLabel->{"t","ϕ","W"},AxesStyle->36,ImageSize->1000,ViewPoint->{4, -10, 5}]

Plot3D[(1/π)NIntegrate[W[t,1.8,0,0,3,0.2,0,θ,μ],{θ,0,π}],{t,0,5},{μ,0,6},PlotRange->All,AxesLabel->{"t","μ","W"},AxesStyle->36,ImageSize->1000,ViewPoint->{4, -10, 5}]

Plot3D[(1/π)NIntegrate[ηnew[t,1.8,0,0,3,0.2,0,θ,μ],{θ,0,π}],{t,0,5},{μ,0,6},PlotRange->All,AxesLabel->{"t","μ","η"},AxesStyle->36,ImageSize->1000,ViewPoint->{4, -10, 5},PlotStyle->Cyan]

(*t,r,n1,n2,T,Rr,ϕ,θ,μ*)

Plot3D[(1/π)NIntegrate[W[t,r,0,0,3,0.2,0,θ,0],{θ,0,π}],{t,0,5},{r,0,1.8},PlotRange->All,AxesLabel->{"t","r","W"},AxesStyle->36,ImageSize->1000,ViewPoint->{4, -10, 5}]

Plot3D[(1/π)NIntegrate[ηnew[t,r,0,0,3,0.2,0,θ,0],{θ,0,π}],{t,0,5},{r,0,1.8},PlotRange->All,AxesLabel->{"t","r","η"},AxesStyle->36,ImageSize->1000,ViewPoint->{4, -10, 5},PlotStyle->Cyan]

Plot3D[W[t,1.8,n1,0,3,0.2,0,π/4,1],{t,0,5},{n1,0,3},PlotRange->All,AxesLabel->{"t","(Subscript[n, 1])","W"},AxesStyle->36,PlotPoints->75,ImageSize->1000,ViewPoint->{4, -10, 5}]

Plot3D[ηnew[t,1.8,n1,0,3,0.2,0,π/4,1],{t,0,5},{n1,0,3},PlotRange->All,AxesLabel->{"t","(Subscript[n, 1])","η"},AxesStyle->36,PlotPoints->75,ImageSize->1000,ViewPoint->{4, -10, 5}, PlotStyle->Cyan]

Plot3D[(1/π)NIntegrate[W[t,1.8,0,n2,3,0.2,0,π/4,1],{t,0,5}],{n2,0,3},PlotRange->All,AxesLabel->{"t","n2","W"},AxesStyle->22,PlotPoints->75,ImageSize->1000,ViewPoint->{4, -10, 5}]

(*W[t,r,n1,n2,T,Rr,ϕ,θ,μ]*)

Plot3D[(1/π)NIntegrate[W[t,1.8,0,0,3,Rr,0,θ,0],{θ,0,π}],{t,0,5},{Rr,0,2},PlotRange->All,AxesLabel->{"t","R","W"},AxesStyle->36,ImageSize->1000,ViewPoint->{4, -10, 5}]