The code for semidefinite optimization

QOT[ρ_,ω_,cost_]:=Module[{σ, A,m,C,c,sol,Π},
σ={{{1,0},{0,1}},{{0,1},{1,0}},{{0,-I},{I,0}},{{1,0},{0,-1}}};
A=Flatten[Table[KroneckerProduct[σ[[j]],σ[[k]]],{j,1,4},{k,1,4}],1];
m=Length[cost];
C=Sum[MatrixPower[KroneckerProduct[cost[[l]],IdentityMatrix[2]]-KroneckerProduct[IdentityMatrix[2],Transpose[cost[[l]]]],2],{l,1,m}];
c=Table[Tr[A[[n]].C],{n,1,16}];
sol=SemidefiniteOptimization[c.{1/4,x12,x13,x14,x21,x22,x23,x24,x31,x32,x33,x34,x41,x42,x43,x44},
{1/4 σ[[1]]+x21 σ[[2]]+x31 σ[[3]]+x41 σ[[4]]==ω/2,1/4 σ[[1]]+x12 σ[[2]]+x13 σ[[3]]+x14 σ[[4]]==Transpose[ρ]/2,
 VectorGreaterEqual[{{1/4,x12,x13,x14,x21,x22,x23,x24,x31,x32,x33,x34,x41,x42,x43,x44}.A,0}, {"SemidefiniteCone", 4}] },
{x12,x13,x14,x21,x22,x23,x24,x31,x32,x33,x34,x41,x42,x43,x44}];
Π=Chop[A[[1]]/4+Sum[sol[[n,2]]A[[n+1]],{n,1,15}]];
{Sqrt[Re[Tr[C.Π]]],MatrixForm[Π]}
]

Numerics on the "clock and shift" cost operator

P[x_,y_,z_]:= 1/2{{1+z,x- I y},{x+I y,1-z}}

C13={{{0,1},{1,0}},{{1,0},{0,-1}}};

Case 1: y=+-1/2

T1={Table[QOT[P[x Sqrt[3/8],1/2,z Sqrt[3/8]], P[0,1/2,0],C13][[1]],{x,-1,1,0.1},{z,-1,1,0.1}],
Table[QOT[P[x Sqrt[3/8],-1/2,z Sqrt[3/8]], P[0,1/2,0],C13][[1]],{x,-1,1,0.1},{z,-1,1,0.1}]};

ListPlot3D[
T1,
DataRange ->{{-Sqrt[3/8], Sqrt[3/8]}, {-Sqrt[3/8], Sqrt[3/8]}},
PlotLegends->{"y=1/2","y=-1/2"},
AxesLabel -> {x,z}
]

Case 2: y=+-4/5

T2={Table[QOT[P[x Sqrt[9/2]/5,4/5,z Sqrt[9/2]/5], P[0,1/2,0],C13][[1]],{x,-1,1,0.1},{z,-1,1,0.1}],
Table[QOT[P[x Sqrt[9/2]/5,-4/5,z Sqrt[9/2]/5], P[0,1/2,0],C13][[1]],{x,-1,1,0.1},{z,-1,1,0.1}]};

ListPlot3D[
T2,
DataRange ->{{- Sqrt[9/2]/5, Sqrt[9/2]/5}, {- Sqrt[9/2]/5, Sqrt[9/2]/5}},
PlotLegends->{"y=4/5","y=-4/5"},
AxesLabel -> {x,z}
]

Case 3: y=+-1/9

T3={Table[QOT[P[t Sqrt[40/81],-1/9,t Sqrt[40/81]], P[0,1/2,0],C13][[1]],{t,-1,1,0.1}],
Table[QOT[P[t Sqrt[40/81],1/9,t Sqrt[40/81]], P[0,1/2,0],C13][[1]],{t,-1,1,0.1}]};

ListPlot[
T3,
DataRange ->{- Sqrt[40/81], Sqrt[40/81]},
AxesLabel ->"x=z"
]

ListPlot