Simulating Ising-Higgs lattice gauge theory on (ℤ^2)

Parameters

L=80;  (* box size *)
n=400; (* number of updates *)
β=2.2;
κ=0.5;

Selects colors to be used in plots

(* sets up colors to be used for the plots *)
colorscheme = "CandyColors"; (* BrightBands, CandyColors, RustTones also work well *)
edgecolor1 = Blend[{ColorData[colorscheme,0.65],White},0.4]; 
edgecolor2 = ColorData[colorscheme,0.8];
edgecolor3 =ColorData[colorscheme,0.5]; 
edgecolor4 = ColorData[colorscheme,0.0];
plaquettecolor = ColorData[colorscheme,0.0];
vertexcolor = ColorData[colorscheme,0.15];
edgecolor5= Blend[{ColorData[colorscheme,0.75],White},0.4];

Defines useful functions

A function which returns the plaquette configuration corresponding to a given edge configuration

vortices[conf_] := With[{L = Length[conf]},Table[ Mod[conf[[i,j,1]]+conf[[i,j+1,2]]+conf[[i+1,j,1]]+conf[[i,j,2]],2] ,{i,1,L-1},{j,1,L-1}]];

A few technical functions which will later be used for the Gibbs samples

hvortices[vort_] :=With[{L = Length[vort]+1},Transpose[Append[Transpose[Prepend[vort,ConstantArray[0,L-1]]+Append[vort,ConstantArray[0,L-1]]],ConstantArray[0,L]]]];

vvortices[vort_]  := With[{L = Length[vort]+1},Append[Transpose[Append[Transpose[vort],ConstantArray[0,L-1]]]+
Transpose[Prepend[Transpose[vort],ConstantArray[0,L-1]]],ConstantArray[0,L]]]; 

edgeweightaroundpoint[conf_]:=With[{we=1-2conf,L=Length[conf]},
Map[Total,we,{2}]+Prepend[we[[1;;L-1,All,2]],ConstantArray[0,L]]+Map[Prepend[#,0]&,we[[All,1;;L-1,1]],{1}]];

totalvertexedgeweight[conf_,κ_] := Exp[κ edgeweightaroundpoint[conf]];

A function which returns a list of all frustrated edges in a configuration, as a list of endpoints of lines, and a matching color assignment function (using the edge clusters)

presentedges[conf_] :=  With[{L=Length[conf]},Join[Select[Flatten[Table[If[conf[[i,j,1]]==1,{{j,i},{j+1,i}},{}],{i,1,L},{j,1,L-1}] ,1],UnsameQ[#,{}]&],
Select[Flatten[Table[If[conf[[i,j,2]]==1,{{j,i},{j,i+1}},{}],{i,1,L-1},{j,1,L}] ,1],UnsameQ[#,{}]&]]];

presentedgescolors[conf_, edgeclusters_] := Map[ColorData[colorscheme,#]&,clusteredgecolorsforpresentedges[conf, edgeclusters]];

A function which assigns an integer for each edge in the configuration, which is the same for edges which belong to the same edge cluster, thus identifying all edge clusters

clusters[conf_]:=(

L = Length[conf];
c = Table[If[conf[[i,j,k]]==1,L^2 i+L j+k,0],{i,1,L},{j,1,L},{k,1,2}];
c0=0*c;

While[ c0≠c,
c0=c;

(* compared edges: = *)
c[[All,All,1]] = Map[Max[#]*If[#[[1]]>0,1,0]&,Table[ {c[[i,j,1]],If[i>1,c[[i-1,j,1]],0],If[i<L,c[[i+1,j,1]],0]},{i,1,L},{j,1,L}],{2}]; 

(*compared edges: |_ *)
c[[All,All,1]] = Map[Max[#]*If[#[[1]]>0,1,0]&,c,{2}];  
c[[All,All,2]] = Map[Max[#]*If[#[[2]]>0,1,0]&,c,{2}];  

(* compared edges: || *)
c[[All,All,2]] = Map[Max[#]*If[#[[1]]>0,1,0]&,Table[ {c[[i,j,2]],If[j>1,c[[i,j-1,2]],0],If[j<L,c[[i,j+1,2]],0]},{i,1,L},{j,1,L}],{2}];

(*compared edges: _| *)
c[[All,All,1]]=Map[Max[#]*If[#[[1]]>0,1,0]&,Table[{c[[i,j,1]],If[j<L,c[[i,j+1,2]],0]},{i,1,L},{j,1,L}],{2}] ;
c[[All,All,2]]=Map[Max[#]*If[#[[1]]>0,1,0]&,Table[{c[[i,j,2]],If[j>1,c[[i,j-1,1]],0]},{i,1,L},{j,1,L}],{2}];

(*compared edges: |- *) 
c[[All,All,1]]=Map[Max[#]*If[#[[1]]>0,1,0]&,Table[{c[[i,j,1]],If[i>1,c[[i-1,j,2]],0]},{i,1,L},{j,1,L}],{2}];
c[[All,All,2]]=Map[Max[#]*If[#[[1]]>0,1,0]&,Table[{c[[i,j,2]],If[i<L,c[[i+1,j,1]],0]},{i,1,L},{j,1,L}],{2}];

(*compared edges: -| *)
c[[All,All,1]]=Map[Max[#]*If[#[[1]]>0,1,0]&,Table[{c[[i,j,1]],If[i>1&&j<L,c[[i-1,j+1,2]],0]},{i,1,L},{j,1,L}],{2}];
c[[All,All,2]]=Map[Max[#]*If[#[[1]]>0,1,0]&,Table[{c[[i,j,2]],If[i<L&&j>1,c[[i+1,j-1,1]] ,0]},{i,1,L},{j,1,L}],{2}];
];

finalcolors = Sort[DeleteDuplicates[Flatten[c]]];
permutedcolors = Prepend[Permute[Range[Length[finalcolors]-1],RandomPermutation[Length[finalcolors]-1]],0];
Do[c=c/.{finalcolors[[i]]->permutedcolors[[i]]},{i,1,Length[finalcolors]}] ;
c);

A function which assign edge cluster numbers to the plaquettes

(* lists which cluster each plaquette belongs to*) 
plaquetteclusters[edgeclusters_] :=With[{L = Length[edgeclusters]},Map[Drop[#,-1]&,Drop[ Table[Max[{edgeclusters[[i,j,1]],edgeclusters[[i,j,2]],If[j<L,edgeclusters[[i,j+1,2]],0],If[i<L,edgeclusters[[i+1,j,1]],0]}],{i,1,L},{j,1,L}],-1]]];

A function which returns a list of the edge clusters which supports vortices

edgeclusterswhichsupportvortices[configuration_, ec_]:=Drop[Sort[DeleteDuplicates[Flatten[plaquetteclusters[ec]*vortices[configuration]]]],1];

A function which returns a coloring of the present edges relative to the edge clusters

clusteredgecolorsforpresentedges[conf_,edgeclusters_] := With[{L=Length[conf],hue = (edgeclusters-1)/(Max[edgeclusters]-1) },Join[Select[Flatten[Table[If[conf[[i,j,1]]==1,hue[[i,j,1]],{}],{i,1,L},{j,1,L-1}] ,1],UnsameQ[#,{}]&],
Select[Flatten[Table[If[conf[[i,j,2]]==1,hue[[i,j,2]],{}],{i,1,L-1},{j,1,L}] ,1],UnsameQ[#,{}]&]]];

Functions which plots the edges in a configuration

(* color is either a single color, or one color for each edge *)edgeplotcolor[conf_,color_] :=  With[{L=Length[conf]},ListPlot[presentedges[conf]
, Joined->True,Axes->False,PlotStyle->color,PlotRange->{{1-0.2,L+0.2},{1-0.2,L+0.2}},AspectRatio->1,Frame->True, FrameTicks->None]];

Functions which plots the plaquettes in a configuration

plaquetteplot[conf_] := With[{L=Length[conf]},ListPlot[Select[Flatten[Table[ If[Mod[conf[[i,j,1]]+conf[[i,j+1,2]]+conf[[i+1,j,1]]+conf[[i,j,2]],2]==1,{j+0.5,i+0.5},{}] ,{i,1,L-1},{j,1,L-1}],1],UnsameQ[#,{}]&]
,PlotMarkers->{Graphics[{plaquettecolor,Rectangle[{0,0}]}],250/L},PlotRange->{{1-0.2,L+0.2},{1-0.2,L+0.2}},AspectRatio->1,Frame->True, FrameTicks->None]];

Functions which converts from unitary gauge to the general measure

verticalhiggscontribution[higgsconfig_]:=Map[Drop[#,-1]&,Drop[higgsconfig,1] + Drop[higgsconfig,-1] ,{1}];
horisontalhiggscontribution[higgsconfig_] := Drop[Map[Drop[#,1]&,higgsconfig,{1}]+Map[Drop[#,-1]&,higgsconfig,{1}],-1];
nonunitaryconfiguration[conf_,higgsconfig_]:=
Table[Mod[conf[[i,j,k]]+If[k==1,(horisontalhiggscontribution[higgsconfig])[[i,j]],(verticalhiggscontribution[higgsconfig])[[i,j]]],2],{i,1,Length[conf]},{j,1,Length[conf]},{k,1,2}];

Functions which plots a Higgs field configuration

higgsplot[higgsconfig_]:=ListPlot[Select[Flatten[Table[If[higgsconfig[[i,j]]==1,{i-1,j-1},{}],{i,1,Length[higgsconfig]},{j,1,Length[higgsconfig]}] ,1],UnsameQ[#,{}]&],PlotStyle->vertexcolor];

The Gibbs sampler

sampleconfiguration[configuration_,β_,κ_,n_]:=(
index=-1; 
c=configuration;
L = Length[configuration];

Do[
index = Mod[index+1,8];
 

conf0 = c;
conf1 = c;
 
 
(* The Markov chain below performs two different types of updates. If "index" is equal to 0,1,2, or 3, then the update corresponds to a gauge transform of the model, and if "index" is equal to 4,5,6, or 7, then the update corresponsds to resampling individual edges. Both types of updates are performed in a chessboard pattern to speed up convergence. *)

If[index==0,

(* proposed updated configuration *)
conf1[[Range[1,L,2],All, 1]] = 1-conf1[[Range[1,L,2],All, 1]] ;
conf1[[All,Range[1,L,2], 2]] = 1-conf1[[All,Range[1,L,2], 2]] ;

(* calculates the total weight around each vertex *)
vertexweight0 = totalvertexedgeweight[conf0,κ];
vertexweight1 = totalvertexedgeweight[conf1,κ];

(* use randomness to choose which configuration we should use around each "active" vertex *) 
rvertexlist= Table[RandomChoice[{vertexweight0[[i,j]],vertexweight1[[i,j]]}->{0,1}], {i,1,L},{j,1,L}];

(* use the previous step to indicate which configuration we should use for each edge *) 
redgelist = 0*conf0;
redgelist[[Range[1,L,2],Range[1,L,2],1]]=rvertexlist[[Range[1,L,2],Range[1,L,2]]];
redgelist[[Range[1,L,2],Range[1,L,2],2]]=rvertexlist[[Range[1,L,2],Range[1,L,2]]];
redgelist[[Range[1,L,2],Range[2,L-1,2],1]]=rvertexlist[[Range[1,L,2],Range[3,L,2]]]; 
redgelist[[Range[2,L-1,2],Range[1,L,2],2]]=rvertexlist[[Range[3,L,2],Range[1,L,2]]]; 

(* updates according to the randomness at sites corr. to step. *) 
c = conf0*(1-redgelist) + conf1*redgelist;
]; 


If[index==1,

(* proposed updated configuration *)
conf1[[Range[2,L,2],All, 1]] = 1-conf1[[Range[2,L,2],All, 1]] ;
conf1[[All,Range[1,L,2], 2]] = 1-conf1[[All,Range[1,L,2], 2]] ;

(* calculates the total weight around each vertex *)
vertexweight0 = totalvertexedgeweight[conf0,κ];
vertexweight1 = totalvertexedgeweight[conf1,κ];

(* use randomness to choose which configuration we should use around each "active" vertex *) 
rvertexlist= Table[RandomChoice[{vertexweight0[[i,j]],vertexweight1[[i,j]]}->{0,1}], {i,1,L},{j,1,L}];

(* use the previous step to indicate which configuration we should use for each edge *) 
redgelist = 0*conf0;
redgelist[[Range[2,L,2],Range[1,L,2],1]]=rvertexlist[[Range[2,L,2],Range[1,L,2]]];
redgelist[[Range[2,L,2],Range[1,L,2],2]]=rvertexlist[[Range[2,L,2],Range[1,L,2]]];
redgelist[[Range[1,L,2],Range[1,L,2],2]]=rvertexlist[[Range[2,L,2],Range[1,L,2]]]; 
redgelist[[Range[2,L,2],Range[2,L-1,2],1]]=rvertexlist[[Range[2,L,2],Range[3,L,2]]];  

(* updates according to the randomness at sites corr. to step. *) 
c = conf0*(1-redgelist) + conf1*redgelist;
]; 


If[index==2,

(* proposed updated configuration *)
conf1[[Range[1,L,2],All, 1]] = 1-conf1[[Range[1,L,2],All, 1]] ;
conf1[[All,Range[2,L,2], 2]] = 1-conf1[[All,Range[2,L,2], 2]] ;

(* calculates the total weight around each vertex *)
vertexweight0 = totalvertexedgeweight[conf0,κ];
vertexweight1 = totalvertexedgeweight[conf1,κ];

(* use randomness to choose which configuration we should use around each "active" vertex *) 
rvertexlist= Table[RandomChoice[{vertexweight0[[i,j]],vertexweight1[[i,j]]}->{0,1}], {i,1,L},{j,1,L}];

(* use the previous step to indicate which configuration we should use for each edge *) 
redgelist = 0*conf0;
redgelist[[Range[1,L,2],Range[2,L,2],1]]=rvertexlist[[Range[1,L,2],Range[2,L,2]]];
redgelist[[Range[1,L,2],Range[2,L,2],2]]=rvertexlist[[Range[1,L,2],Range[2,L,2]]];  
redgelist[[Range[1,L,2],Range[1,L,2],1]]=rvertexlist[[Range[1,L,2],Range[2,L,2]]]; 
redgelist[[Range[2,L-1,2],Range[2,L,2],2]]=rvertexlist[[Range[3,L,2],Range[2,L,2]]];  

(* updates according to the randomness at sites corr. to step. *) 
c = conf0*(1-redgelist) + conf1*redgelist;
]; 


If[index==3,

(* proposed updated configuration *)
conf1[[Range[2,L,2],All, 1]] = 1-conf1[[Range[2,L,2],All, 1]] ;
conf1[[All,Range[2,L,2], 2]] = 1-conf1[[All,Range[2,L,2], 2]] ;

(* calculates the total weight around each vertex *)
vertexweight0 = totalvertexedgeweight[conf0,κ];
vertexweight1 = totalvertexedgeweight[conf1,κ];

(* use randomness to choose which configuration we should use around each "active" vertex *)
rvertexlist= Table[RandomChoice[{vertexweight0[[i,j]],vertexweight1[[i,j]]}->{0,1}], {i,1,L},{j,1,L}];

(* use the previous step to indicate which configuration we should use for each edge *)
redgelist = 0*conf0;
redgelist[[Range[2,L,2],Range[2,L,2],1]]=rvertexlist[[Range[2,L,2],Range[2,L,2]]];
redgelist[[Range[2,L,2],Range[2,L,2],2]]=rvertexlist[[Range[2,L,2],Range[2,L,2]]];
redgelist[[Range[2,L,2],Range[1,L,2],1]]=rvertexlist[[Range[2,L,2],Range[2,L,2]]]; 
redgelist[[Range[1,L,2],Range[2,L,2],2]]=rvertexlist[[Range[2,L,2],Range[2,L,2]]];

(* updates according to the randomness at sites corr. to step. *)
c = conf0*(1-redgelist) + conf1*redgelist;
]; 

 
 If[index==4,
conf0[[Range[1,L,2],All, 1]] = 0*conf0[[Range[1,L,2],All, 1]] ;
conf1[[Range[1,L,2], All, 1]] = 1-0*conf0[[Range[1,L,2],All, 1]] ; 

weights0 = Exp[β (2-2 hvortices[vortices[conf0]]) + κ];
weights1 = Exp[β (2-2 hvortices[vortices[conf1]]) - κ];
 
newspins = Table[RandomChoice[{weights0[[i,j]],weights1[[i,j]]}->{0,1}], {i,1,L},{j,1,L}];

c[[Range[1,L,2],All, 1]] = newspins[[Range[1,L,2],All]]; 
];

  If[index==5,
conf0[[Range[2,L,2],All, 1]] = 0*conf0[[Range[2,L,2],All, 1]] ;
conf1[[Range[2,L,2], All, 1]] = 1-0*conf0[[Range[2,L,2],All, 1]] ; 

weights0 = Exp[β (2-2 hvortices[vortices[conf0]]) + κ];
weights1 = Exp[β (2-2 hvortices[vortices[conf1]]) - κ];
 
newspins = Table[RandomChoice[{weights0[[i,j]],weights1[[i,j]]}->{0,1}], {i,1,L},{j,1,L}];

c[[Range[2,L,2],All, 1]] = newspins[[Range[2,L,2],All]]; 
]; 

If[index==6,
conf0[[All,Range[1,L,2], 2]] = 0*conf0[[All,Range[1,L,2],2]] ;
conf1[[All, Range[1,L,2], 2]] = 1-0*conf0[[All,Range[1,L,2],2]] ; 

weights0 = Exp[β (2-2 vvortices[vortices[conf0]]) + κ];
weights1 = Exp[β (2-2 vvortices[vortices[conf1]]) - κ];
 
newspins = Table[RandomChoice[{weights0[[i,j]],weights1[[i,j]]}->{0,1}], {i,1,L},{j,1,L}];

c[[All, Range[1,L,2],2]] = newspins[[All, Range[1,L,2]]]; 

];

If[index==7,
conf0[[All,Range[2,L,2], 2]] = 0*conf0[[All,Range[2,L,2],2]] ;
conf1[[All, Range[2,L,2], 2]] = 1-0*conf0[[All,Range[2,L,2],2]] ; 

weights0 = Exp[β (2-2 vvortices[vortices[conf0]]) + κ];
weights1 = Exp[β (2-2 vvortices[vortices[conf1]]) - κ];
 
newspins = Table[RandomChoice[{weights0[[i,j]],weights1[[i,j]]}->{0,1}], {i,1,L},{j,1,L}];

c[[All, Range[2,L,2],2]] = newspins[[All, Range[2,L,2]]]; 

];   

,8 n];
c);

Picks an initial distribution

(* Here one could start from either a zero configuration, or a random edge configuration. *)
configuration=0*Table[RandomInteger[],{x,1,L},{y,1,L},{dir,1,2}] ;

Samples configurations

configuration = sampleconfiguration[configuration,β,κ,n]; 
higgsconfiguration = RandomInteger[1,{L+1,L+1}];

Draws the sampled configuration

Draws the sampled configuration

Show[edgeplotcolor[configuration, edgecolor1],plaquetteplot[configuration]]

Draws a corresponding configuration in non-unitary gauge

Show[edgeplotcolor[nonunitaryconfiguration[configuration,higgsconfiguration], edgecolor5],higgsplot[higgsconfiguration]]

Colors the sampled configuration according to edge clusters

edgecolors = presentedgescolors[configuration, clusters[configuration]];
Show[edgeplotcolor[configuration, edgecolors],plaquetteplot[configuration]]

Draws illustrations of the coupling argument

Setup

configurationIsing  = sampleconfiguration[configuration,16,κ,n]; 
jointconfiguration = Table[Max[configuration[[i,j,k]],configurationIsing[[i,j,k]]],{i,1,L},{j,1,L},{k,1,2}];
ec = clusters[jointconfiguration];
cwsv = edgeclusterswhichsupportvortices[configuration, ec];
edgesinvortexclusters = With[{L = Length[jointconfiguration]},Table[
If[MemberQ[cwsv,ec[[i,j,k]]],jointconfiguration[[i,j,k]],0]
,{i,1,L},{j,1,L},{k,1,2}]];

Draws the LGT configuration

Show[edgeplotcolor[configuration, edgecolor2],plaquetteplot[configuration]]

Draws the Ising configuration

Show[edgeplotcolor[configurationIsing, edgecolor3]]

Draws both configurations simultaneously

Show[edgeplotcolor[configuration*configurationIsing,edgecolor4],
edgeplotcolor[configuration*(1-configurationIsing),edgecolor2],
edgeplotcolor[(1-configuration)*configurationIsing,edgecolor3],
plaquetteplot[configuration]
]

Draws the the edge clusters in both configurations that are connected to vortices in the joint graph

Show[edgeplotcolor[configuration*configurationIsing*edgesinvortexclusters,edgecolor4],
edgeplotcolor[configuration*(1-configurationIsing)*edgesinvortexclusters,edgecolor2],
edgeplotcolor[(1-configuration)*configurationIsing*edgesinvortexclusters,edgecolor3],
plaquetteplot[configuration]
]

Draws the edge clusters in the LGT configuration that are not connected to vortices in the joint graph

Show[
edgeplotcolor[configuration*(1-edgesinvortexclusters),edgecolor1]

Draws the edge clusters in the Ising configuration that are not connected to vortices in the joint graph

Show[
edgeplotcolor[configurationIsing*(1-edgesinvortexclusters),edgecolor3]
]

Draws a representation of the final coupling

Show[edgeplotcolor[configuration*configurationIsing*edgesinvortexclusters,edgecolor4],
edgeplotcolor[configuration*(1-configurationIsing)*edgesinvortexclusters,edgecolor2],
edgeplotcolor[(1-configuration)*configurationIsing*edgesinvortexclusters,edgecolor3],
plaquetteplot[configuration],
edgeplotcolor[configuration*(1-edgesinvortexclusters), edgecolor1]]