Oscillator semigroup, quantum measurements, and fractals on squeezed States
by Arkadiusz Jadczyk

Hessian matrix for h(q,p)=(q^2)

h1=D[q^2,{{q,p}},{{q,p}}]

Hessian matrix for h(q,p)=(p^2)

h2=D[p^2,{{q,p}},{{q,p}}]

Hessian matrix for h(q,p)=(q+((p))^2)/2

h3=D[(q+p)^2/2,{{q,p}},{{q,p}}]

Hessian matrix for h(q,p)=(q-((p))^2)/2

h4=D[(q-p)^2/2,{{q,p}},{{q,p}}]

The corresponding generators of one parameter groups of symplectic transformations:

J={{0,1},{-1,0}};
H1=J.h1
H2=J.h2
H3=J.h3
H4=J.h4

Cayley transformation to complex su(1,1) generators

c={{1,I},{1,-I}}/Sqrt[2];
ic=Inverse[c]
cH1=c.H1.ic
cH2=c.H2.ic
cH3=c.H3.ic
cH4=c.H4.ic

One parameter groups of SU(1,1) transformations

u1[t_]=MatrixExp[t*cH1]
u2[t_]=MatrixExp[t*cH2]
u3[t_]=MatrixExp[t*cH3]
u4[t_]=MatrixExp[t*cH4]

One parameter semigroups of contractions of the unit disk

s1[ϵ_]=u1[I*ϵ]
s2[ϵ_]=u2[I*ϵ]
s3[ϵ_]=u3[I*ϵ]
s4[ϵ_]=u4[I*ϵ]

Deformation drawing function:

pl[a_]:=Module[{plr,plf,plc,p0,p10},
p10:=ParametricPlot[{Re[(a[[2]][[2]]*r+a[[2]][[1]])/(a[[1]][[2]]*r+a[[1]][[1]])],Im[(a[[2]][[2]]*r+a[[2]][[1]])/(a[[1]][[2]]*r+a[[1]][[1]])]},{r,-1,1},PlotStyle->{Thickness[0.01],Black}];
plr:=ParametricPlot[Table[{Re[(a[[2]][[2]]*(r*Cos[i*Pi/18]+I*r*Sin[i*Pi/18])+a[[2]][[1]])/(a[[1]][[2]]*(r*Cos[i*Pi/18]+I*r*Sin[i*Pi/18])+a[[1]][[1]])],Im[(a[[2]][[2]]*(r*Cos[i*Pi/18]+I*r*Sin[i*Pi/18])+a[[2]][[1]])/(a[[1]][[2]]*(r*Cos[i*Pi/18]+I*r*Sin[i*Pi/18])+a[[1]][[1]])]},{i,0,35}],{r,0,1},PlotRange->All];
plf:=ParametricPlot[Table[{Re[(a[[2]][[2]]*(i*Cos[ϕ]/10+I*i*Sin[ϕ]/10)+a[[2]][[1]])/(a[[1]][[2]]*(i*Cos[ϕ]/10+I*i*Sin[ϕ]/10)+a[[1]][[1]])],Im[(a[[2]][[2]]*(i*Cos[ϕ]/10+I*i*Sin[ϕ]/10)+a[[2]][[1]])/(a[[1]][[2]]*(i*Cos[ϕ]/10+I*i*Sin[ϕ]/10)+a[[1]][[1]])]},{i,0,10}],{ϕ,0,2Pi}];
plc:=Graphics[Circle[]];
Show[{plr,plf,plc,p10}]]

Disk deformation, compression regions, and fixed points for each of the trenaformations, for ϵ=0.5.  The thick line is the image of the interval  [-1,1] on the real axis.

pl[s1[0.5]]
pl[s2[0.5]]
pl[s3[0.5]]
pl[s4[0.5]]

Fractal generation

The function f below selects randomly, with equal probabilities 1/4, one of the four contractive Mobius transformations corresponding to hamiltonian functions (q^2),(p^2),(((q+p))^2)/2,(((q-p))^2)/2, and applies it to the current point - the second argument of the function. The first argument is the nonnegative value determining the amount of contraction.

f=Compile[{{ϵ,_Real},{xy,_Real,1}},Module[{s,z,z1,r,res},s={{{1+ϵ,ϵ},{-ϵ,1-ϵ}},{{1+ϵ,-ϵ},{ϵ,1-ϵ}},{{1+ϵ,I*ϵ},{I*ϵ,1-ϵ}},{{1+ϵ,-I*ϵ},{-I*ϵ,1-ϵ}}};z=xy[[1]]+I*xy[[2]];r=RandomInteger[{1,4}];
z1=(s[[r]][[2]][[2]]*z+s[[r]][[2]][[1]])/(s[[r]][[1]][[2]]*z+s[[r]][[1]][[1]]);res={Re[z1],Im[z1]}]];

In fractal generation below we use 10000000 points, the staring point is always the same {0.1,0.2}. Instead of displaying the points we divide the unit square into 1000x1000 cells and count the points in each cell. The higher the number of points in a given cell, the brigter is its appearance in the image. The log function is used in order to see more structure. Black means never visited pixel. Of course the outside region of the disk is always black, because our transformations map disk points into the their compression regions inside the disk.

Equal probabilities, ϵ=0.5:

n05=NestList[f[0.5,#] &,{0.1,0.2},10000000];bc05=BinCounts[n05,{-1,1,0.002},{-1,1,0.002}];
bc05l=N[Log[bc05+1]];
lp05=ArrayPlot[bc05l,ColorFunction->GrayLevel,DataRange->{{-1,1},{-1,1}}];
cg=Graphics[{White,Circle[]}];
lp05s=Show[{lp05,cg}]

Equal probabilities, ϵ=0.4:

n04=NestList[f[0.4,#] &,{0.1,0.2},10000000];bc04=BinCounts[n04,{-1,1,0.002},{-1,1,0.002}];
bc04l=N[Log[bc04+1]];
lp04=ArrayPlot[bc04l,ColorFunction->GrayLevel,DataRange->{{-1,1},{-1,1}}];
c=Graphics[{White,Circle[]}];
lp04s=Show[{lp04,c}]

Equal probabilities, ϵ=0.8:

n08=NestList[f[0.8,#] &,{0.1,0.2},10000000];bc08=BinCounts[n08,{-1,1,0.002},{-1,1,0.002}];
bc08l=N[Log[bc08+1]];
lp08=ArrayPlot[bc08l,ColorFunction->GrayLevel,DataRange->{{-1,1},{-1,1}}];
c=Graphics[{White,Circle[]}];
lp08s=Show[{lp08,c}]

There is a natural function that can be used for place dependent probabilities. For a semigroup matrix {{λ,μ},{ν,ρ}} and a point z in th disk, the norm squared of the resulting Hilbert space state is  (1-z*z)^(1/2)/(|μz+λ|^2-|ρz+ν|^2)^(1/2). When using place dependent probabilities generating the fractal takes longer. Here is our fractal for ϵ=0.5 with place dependent probabilities proportional to  (|μz+λ|^2-|ρz+ν|^2)^(-1/2). There is more structure displayed than in the equal probabilities case.

f1=Compile[{{ϵ,_Real},{xy,_Real,1}},Module[{s,z,z1,r,pr,res},s={{{1+ϵ,ϵ},{-ϵ,1-ϵ}},{{1+ϵ,-ϵ},{ϵ,1-ϵ}},{{1+ϵ,I*ϵ},{I*ϵ,1-ϵ}},{{1+ϵ,-I*ϵ},{-I*ϵ,1-ϵ}}};z=xy[[1]]+I*xy[[2]];
pr=Table[(Abs[s[[i]][[1]][[2]]*z+s[[i]][[1]][[1]]]^2-Abs[s[[i]][[2]][[2]]*z+s[[i]][[2]][[1]]]^2)^(-1/2),{i,1,4}];
r=RandomChoice[pr->{1,2,3,4}];
z1=(s[[r]][[2]][[2]]*z+s[[r]][[2]][[1]])/(s[[r]][[1]][[2]]*z+s[[r]][[1]][[1]]);res={Re[z1],Im[z1]}]];

n051=NestList[f1[0.5,#] &,{0.1,0.2},10000000];bc051=BinCounts[n051,{-1,1,0.002},{-1,1,0.002}];
bc051l=N[Log[bc051+1]];
lp051=ArrayPlot[bc051l,ColorFunction->GrayLevel,DataRange->{{-1,1},{-1,1}}];
cg=Graphics[{White,Circle[]}];
lp051s=Show[{lp051,cg}]

Another choice, with probabilities proportional to (|μz+λ|^2-|ρz+ν|^2)^2, the indefinite U(1,1) norm squared of the transformed vector {1,z} changes the measure in the opposite direction - less structure is visible

f2=Compile[{{ϵ,_Real},{xy,_Real,1}},Module[{s,z,z1,r,pr,res},s={{{1+ϵ,ϵ},{-ϵ,1-ϵ}},{{1+ϵ,-ϵ},{ϵ,1-ϵ}},{{1+ϵ,I*ϵ},{I*ϵ,1-ϵ}},{{1+ϵ,-I*ϵ},{-I*ϵ,1-ϵ}}};z=xy[[1]]+I*xy[[2]];
pr=Table[(Abs[s[[i]][[1]][[2]]*z+s[[i]][[1]][[1]]]^2-Abs[s[[i]][[2]][[2]]*z+s[[i]][[2]][[1]]]^2)^2,{i,1,4}];
r=RandomChoice[pr->{1,2,3,4}];
z1=(s[[r]][[2]][[2]]*z+s[[r]][[2]][[1]])/(s[[r]][[1]][[2]]*z+s[[r]][[1]][[1]]);res={Re[z1],Im[z1]}]];

n052=NestList[f2[0.5,#] &,{0.1,0.2},10000000] ; 
bc052=BinCounts[n052,{-1,1,0.002},{-1,1,0.002}];
bc052l=N[Log[bc052+1]];
lp052=ArrayPlot[bc052l,ColorFunction->GrayLevel,DataRange->{{-1,1},{-1,1}}];
cg=Graphics[{White,Circle[]}];
lp052s=Show[{lp052,cg}]

Though the differences are small. Here are the three patterns compared:

gg3=GraphicsGrid[{{lp05,lp051,lp052,cg}},ImageSize->800]

It is also possible to add symplectic transformations to the contractions. Here we add u1[0.5], u2[0.5], u3[0.5], u4[0.5]. Together with s1[0.8], s2[0.8], s3[0.8], s4[0.8], we will have 8 transformations. We will use equal probabilities 1/8.

u8=N[{s1[0.8],s2[0.8],s3[0.8],s4[0.8],u1[0.5],u2[0.5],u3[0.5],u4[0.5]}]
f3=Compile[{{xy,_Real,1}},Module[{s,z,z1,r,res},s=u8;z=xy[[1]]+I*xy[[2]];r=RandomInteger[{1,8}];
z1=(s[[r]][[2]][[2]]*z+s[[r]][[2]][[1]])/(s[[r]][[1]][[2]]*z+s[[r]][[1]][[1]]);res={Re[z1],Im[z1]}]];
n8=NestList[f3[#] &,{0.1,0.2},10000000];
bc8=BinCounts[n8,{-1,1,0.002},{-1,1,0.002}];
bc8l=N[Log[bc8+1]];

lp8s=Show[{ArrayPlot[bc8l,ColorFunction->GrayLevel]}]

Another version, with 0.8 replaced by 2.0, and colored

u81=N[{s1[2.0],s2[2.0],s3[2.0],s4[2.0],u1[0.5],u2[0.5],u3[0.5],u4[0.5]}]
f31=Compile[{{xy,_Real,1}},Module[{s,z,z1,r,res},s=u81;z=xy[[1]]+I*xy[[2]];r=RandomInteger[{1,8}];
z1=(s[[r]][[2]][[2]]*z+s[[r]][[2]][[1]])/(s[[r]][[1]][[2]]*z+s[[r]][[1]][[1]]);res={Re[z1],Im[z1]}]];
n81=NestList[f31[#] &,{0.1,0.2},10000000];
bc81=BinCounts[n81,{-1,1,0.002},{-1,1,0.002}];
bc8l1=N[Log[bc81+1]];
lp81s=Show[{ArrayPlot[bc8l1,ColorFunction->"SunsetColors",ImageSize->1000]}]