Spin squeezing

The following is a Mathematica notebook working with the spin squeezing theory from M. Kitagawa and M. Ueda paper in Phys. Rev. A 47, 5138 (1993) http://dx.doi.org/10.1103/PhysRevA.47.5138

Symmetric states

We are going to study the squeezing properties of spin states made from N spin s=1/2 subsystems. We can do several simplifications. The most obvious one is that we will start and thus stay with the Hilbert space spanned by symmetric states, ((Ket[ϕ])^(⊗N)). This basis is made from states with maximal angular momentum (Ket[(N/2),m]) for m=-(N/2),...,(N/2). We start by writing the representation of angular momentum operators with total spin S=N/2.

Clear[Jz];
Jz[Nspins_]:=Module[{j=Nspins/2,nstates=Nspins+1,n0=Nspins/2+1},
SparseArray[Table[{m+n0,m+n0}->N[m],{m,-j,j}],{nstates,nstates}]];

Clear[Jplus];
Jplus[Nspins_]:=Module[{j=Nspins/2,nstates=Nspins+1,n0=Nspins/2+1},
SparseArray[Table[{m+1+n0,m+n0}->N[(Sqrt[j(j+1)-m(m+1)])],{m,-j,j-1}],{nstates,nstates}]];

Clear[Jminus];
Jminus[Nspins_]:=Module[{j=Nspins/2,nstates=Nspins+1,n0=Nspins/2+1},
SparseArray[Table[{m-1+n0,m+n0}->N[(Sqrt[j(j+1)-m(m-1)])],{m,-j+1,j}],{nstates,nstates}]];

Clear[Jx];
Jx[Nspins_]:=0.5(Jplus[Nspins]+Jminus[Nspins]);
Clear[Jy];
Jy[Nspins_]:=0.5I(Jplus[Nspins]-Jminus[Nspins]);

Jx[1]//MatrixForm

Jz[1]//MatrixForm

Variances

The following are the normalized variances. The problem with this function is that it works with the “canonical” variables and does not give the directions of optimal squeezing.

Var[ψ_]:=
Module[{Nspins=Length[ψ]-1},
Module[{Sxψ=Jx[Nspins].ψ,Syψ=Jy[Nspins].ψ,Szψ=Jz[Nspins].ψ},
(* Angular momentum averages ⟨(Subscript[S, x])⟩,⟨(Subscript[S, y])⟩,⟨(Subscript[S, z])⟩ *)
Module[{meanX=ψ\[Conjugate].Sxψ,meanY=ψ\[Conjugate].Syψ,meanZ=ψ\[Conjugate].Szψ},
(* Variances ((((Subscript[ΔS, α])))^2)=⟨(Subscript[S, α])(Subscript[S, α])⟩-((⟨(Subscript[S, α])⟩)^2) *)
Module[{varX=Sxψ\[Conjugate].Sxψ-(meanX^2),varY=Syψ\[Conjugate].Syψ-(meanY^2),varZ=Szψ\[Conjugate].Szψ-(meanZ^2)},
Re[{meanX,meanY,meanZ,varX,varY,varZ}]]]]]

We solve this by computing the full covariance matrix

Cov[ψ_]:=
(* Compute the average value of the angular momentum operators and covariance matrix *)
Module[{Nspins=Length[ψ]-1},
Module[{Sxψ=Jx[Nspins].ψ,Syψ=Jy[Nspins].ψ,Szψ=Jz[Nspins].ψ},
(* Angular momentum averages ⟨(Subscript[S, x])⟩,⟨(Subscript[S, y])⟩,⟨(Subscript[S, z])⟩ *)
Module[{meanX=ψ\[Conjugate].Sxψ,meanY=ψ\[Conjugate].Syψ,meanZ=ψ\[Conjugate].Szψ},
(* Covariances (Subscript[C, αβ])=⟨(Subscript[S, α])(Subscript[S, β])⟩-⟨(Subscript[S, α])⟩⟨(Subscript[S, β])⟩ *)
Module[{
cov=Re[(({{Sxψ\[Conjugate].Sxψ, Sxψ\[Conjugate].Syψ, Sxψ\[Conjugate].Szψ}, {Syψ\[Conjugate].Sxψ, Syψ\[Conjugate].Syψ, Syψ\[Conjugate].Szψ}, {Szψ\[Conjugate].Sxψ, Szψ\[Conjugate].Syψ, Szψ\[Conjugate].Szψ}}))-(({{meanX}, {meanY}, {meanZ}})).(({{meanX, meanY, meanZ}}))]},
Re[{meanX,meanY,meanZ,cov}]]]]]

With this we can compute the minimal covariance matrices

Clear[MinCov];
MinCov[ψ_]:=
(* Compute the covariance matrix and the average spin *)
Module[{s=Cov[ψ]},
 (* We diagonalize the covariance matrix and add it to the list *)
Join[s[[{1,2,3}]],Sort[Re[Eigenvalues[s[[4]]]]],{s[[4]]}]]

We can also do this along directions orthoginal to the spin orientation

Clear[MinCov2D];
MinCov2D[ψ_]:=
(* We start with the spin covariance matrix *)
Module[{s=Cov[ψ]},
(* Average orientation of spin *)
Module[{v1=Normalize[s[[{1,2,3}]]]},
(* Projector onto space orthogonal to v1 *)
Module[{Port=IdentityMatrix[3]-{v1}\[ConjugateTranspose].{v1}},(* Extract the part of the cov. matrix in the orthogonal space *)
Module[{c=Port.s[[4]].Port},
Module[{λ=Sort[Re[Eigenvalues[c]]]},
Join[s[[{1,2,3}]],λ[[{2,3}]],{s[[4]]}]]]]]]

Coherent spin states and uncertainties

We also need states with well defined angular momentum

Clear[ψ];
ψ[N_,m_]:=Table[If[r==m,1.0,0.0],{r,-N/2,N/2}];
ψ[N_]:=ψ[N,-N/2];

Note that this state is aligned along the Z axis, pointing towards ⟨(Subscript[S, z])⟩=N/2. Note also that the variance of (Subscript[S, z]) is zero, but the other ones are N/4

Var[ψ[30]]

We will also need spin coherent states

Clear[ψcss];
ψcss[Nspins_,θ_,ϕ_]:=
Table[
Which[n==0,((Sin[(θ/2)])^Nspins)Exp[I (Nspins)ϕ],
n==Nspins,((Cos[(θ/2)])^n)Exp[I (Nspins-n)ϕ],
True,(Sqrt[Binomial[Nspins,n]])((Sin[(θ/2)])^(Nspins-n))((Cos[(θ/2)])^n)Exp[I (Nspins-n)ϕ]],{n,0,Nspins}];

Coherent spin states are oriented along various directions. Hence, the variances of the standard angular momentum operators are not the smallest ones

Var[ψcss[30,0.5,1]]

Instead we can use the covariance matrix to show that one angular momentum direction is perfectly well defined and the other two have the standard Heisenberg-limited variance

MinCov[ψcss[30,0.5,1]]

With this we can build the Q representation

Q[θ_,ϕ_,ψ_]:=
Module[{Nspins=Length[ψ]-1},
(((Nspins+1))/(4π))((Abs[ψcss[Nspins,θ,ϕ]\[Conjugate].ψ])^2)];

And finally we have the possibility of plotting variances and uncertainties using this Q function:

Clear[Linspace];
Linspace[a_,b_,Npoints_]:=Table[a+(b-a)/(Npoints-1)i,{i,0,Npoints-1}];

PlotCSS[ψ_]:=
Module[{Nspins=Length[ψ]-1},
Module[{maxoverlap=Quiet[FindMaximum[(Sqrt[Q[θ,ϕ,ψ]]),{θ,0.13},{ϕ,-0.5}]][[1]]},
Module[{myf=Function[{θ,ϕ},(Sqrt[Q[θ,ϕ,ψ]/maxoverlap])]},
SphericalPlot3D[1,{θ,0,π},{ϕ,0,2 π},PlotPoints->30,
ColorFunction->(Darker[White,myf[#4,#5]]&),
PlotRange->Full,ColorFunctionScaling->False,ViewVector->5*Normalize[Var[ψ][[{1,2,3}]]],
AxesLabel->{"x","y","z"},Mesh->3,MeshStyle->Dashed]]]];

We do this for a spin coherent state. Note that PlotCSS points the camera to the point of largest overlap in the sphere of CSS orientations.

PlotCSS[ψcss[10,0,0]]

We change the CSS orientation and the camera is moved:

PlotCSS[ψcss[40,π/2,0.0]]

Squeezing (Jz^2) Hamiltonian

We work with the squeezing Hamiltonian

H=χ (Subscript[J, z]^2)

which stretches a state that is initially polarized along X.

ψsq[Nspins_,μ_]:=
Module[{J=Jz[Nspins]},
Module[{H=(J^2)},
MatrixExp[(-I μ/2) H, ψcss[Nspins,π/2.0,0.0]]]]

PlotCSS[ψsq[40,0.075]]

The following are plots from Ueda’s paper and show a coherent spin state at various degrees of squeezing

Quiet[{PlotCSS[ψsq[40,0.0]],PlotCSS[ψsq[40,0.199]],PlotCSS[ψsq[40,0.399]]}]

Note that ordinary variances do not seem squeezed

Plot[{Var[ψsq[30,μ]][[4]],Var[ψsq[30,μ]][[5]],Var[ψsq[30,μ]][[6]]},{μ,0,2},PlotLegends->{"((Subscript[ΔS, x])^2)","((Subscript[ΔS, y])^2)","((Subscript[ΔS, z])^2)"}]

This is so because the squeezing happens along some direction on the ZY plane which rotates with the Hamiltonian H. One possibility is to compute the covariance matrix and diagonalize it to find the optimal squeezing direction

Plot[{MinCov[ψsq[30,μ]][[4]],MinCov[ψsq[30,μ]][[5]],MinCov[ψsq[30,μ]][[6]]},{μ,0,2}]

This procedure has the problem that it shows zero uncertainty along the direction along which the spin is well defined, which is not the relevant plane. It is better to focus just on the part that happens on the directions which are orthogonal to (⟨(Subscript[S, x])⟩, ⟨(Subscript[S, y])⟩, ⟨(Subscript[S, z])⟩). This is what MinCov2D does:

Plot[{MinCov2D[ψsq[8,μ]][[4]],MinCov2D[ψsq[16,μ]][[4]],MinCov2D[ψsq[32,μ]][[4]]},{μ,0,2}]

We can compute the maximum squeezing for different number of particles

aux=ParallelTable[{Nspins,Quiet[FindMinimum[MinCov2D[ψsq[Nspins,μ]][[4]],{μ,0.3}]][[1]]},
{Nspins,{4,8,16,32,64,128}}]

ListLogLogPlot[aux,Joined->True,Frame->True,FrameLabel->{"N","((ΔS)^2)"},BaseStyle->FontSize->16]

The uncertainty slightly increases but because the value of the total spin increases faster, this results in a total squeezing (N^(-3/2)) as shown in Ueda’s paper:

Quiet[
Show[
ListLogLogPlot[Map[{#1[[1]],#1[[2]]/(((#1[[1]])^2)/4)}&,aux],Joined->True,Frame->True,FrameLabel->{"N","((ΔS)^2)/(S^2)"},BaseStyle->FontSize->16],
LogLogPlot[1/(N^(3/2)),{N,4,256},PlotStyle->Dashed]]
]

Squeezing ((SubPlus[J])^2)-((SubMinus[J])^2) Hamiltonian

We now work with the squeezing Hamiltonian

H=((iχ)/2) ((Subscript[J, +]^2)-(Subscript[J, -]^2))

This Hamiltonian squeezes the (Subscript[J, x]) observable for an initial state polarized along the Z direction

ψsq2[Nspins_,μ_]:=
Module[{Sp=Jplus[Nspins],Sm=Jminus[Nspins]},
Module[{H=(I/2)(Sp.Sp-Sm.Sm)},
MatrixExp[(-I μ/4) H, ψcss[Nspins,0,0.0]]]]

Quiet[{PlotCSS[ψsq2[40,0.0]],PlotCSS[ψsq2[40,0.203]],PlotCSS[ψsq2[40,0.248]]}]

Note that now we can use ordinary variances as there is no “rotation”. The optimally squeezed direction is (Subscript[S, x])

Plot[{Var[ψsq2[30,μ]][[4]],Var[ψsq2[30,μ]][[5]],Var[ψsq2[30,μ]][[6]]},{μ,0,0.5},PlotLegends->{"((Subscript[ΔS, x])^2)","((Subscript[ΔS, y])^2)","((Subscript[ΔS, z])^2)"}]

We can compute the maximum squeezing for different number of particles

aux=ParallelTable[{Nspins,Quiet[FindMinimum[{Var[ψsq2[Nspins,μ]][[4]],0≤μ≤0.5},{μ,0.3}]][[1]]},
{Nspins,{4,8,16,32,64,128}}]

ListLogLogPlot[aux,Joined->True,Frame->True,FrameLabel->{"N","((ΔS)^2)"},BaseStyle->FontSize->16]

The uncertainty tends towards 1/2, so that when we take into account the size of the spin, squeezing decreases as  1/(N^2)

Quiet[
Show[
ListLogLogPlot[Map[{#1[[1]],#1[[2]]/(((#1[[1]])^2)/4)}&,aux],Joined->True,Frame->True,FrameLabel->{"N","((ΔS)^2)/(S^2)"},BaseStyle->FontSize->16],
LogLogPlot[1/(N^2),{N,4,256},PlotStyle->Dashed]]
]