IMPORTANT NOTE: All models in this notebook have been transformed into CSS models. Therefore, instead of being cluster states, each model has stabilizers consisting of either Pauli X or Pauli Z exclusively. CZ gates in a cluster state becomes CNOT in a CSS model. These are completely equivalent ways to write the models, however. To convert between them, just apply Hadamard to each qubit on one of the sublattices.
NOTE 2: All cells should be run in order.

1.Tools

General ERG Tools

(*Most of these adapted from the attached mathematica notebook in Phys.Rev.B 89,075119 (arXiv:1310.4507)*)

Clear[x,y,z,w];translationVariables={x,y,z};

(*displays entries mod 2*)
display2[mat_]:=MatrixForm[Expand[PolynomialMod[mat,2]]]
expandEntry[poly_,x_,n_]:=Module[{f,p,c,xp,col,monos,exps,gen},
f[i_,j_]:=If[i-1==j,1,If[i-1+n==j,x,0]];
gen=Table[f[i,j],{i,1,n},{j,1,n}];
p=Expand[poly];
If[p===0,
(*return zero matrix if the input is zero*)
Table[0,{n},{n}],
(*Else, compute properly*)
col=Collect[p,x];
monos=If[Head[col]===Plus,Apply[List,col],{col}];
exps=Map[Exponent[#,x]&,monos];
c=Power[x,-exps]*monos;
xp=Map[MatrixPower[gen,#]&,exps];
Expand[c.xp]
(*end-if*)]
]
(*coarse grains a matrix along x by n lattice spacings*)
coarseGrain[mat_,{x_,n_}]:=
Apply[Join,Table[
Join[Table[
Flatten[
Map[expandEntry[#,x,n][[nn]]&,mat[[row]]]
],{nn,n}]
],{row,Length[mat]}]
]

(*takes eg x -> 1/x, corresponds to reflection of lattice coords*)
monomialAntipode[f_]:=
Module[{e,v,c,mono},
v=Intersection[Variables[f],translationVariables];
e=Map[Exponent[f,#]&,v];
mono=Apply[Times,Power[v,e]];
c=If[mono===1,f,Coefficient[f,mono]];
c*Apply[Times,Power[v,-e]]]
antipode[f_]:=
Module[{ff,fff},
ff=Expand[f];
If[Head[ff]===Plus,fff=ff,fff={ff}];
Apply[Plus,Map[monomialAntipode,Apply[List,Expand[fff]]]]
]
monomialAntipode[f_[x_]]:=f[1/x]
monomialAntipode[f_[x_]^n_]:=f[1/x]^n
monomialAntipode[mono_*f_[x_]]:=monomialAntipode[mono]*monomialAntipode[f[x]]
monomialAntipode[mono_*f_[x_]^n_]:=monomialAntipode[mono]*monomialAntipode[f[x]]^n

(*apply translationally invariant cnot from src to tgt with a displacement of f. If f is a polynomial apply multiple cnots*)
cNotMatrix[q_,tgt_,src_,f_]:=Module[{m},
m=IdentityMatrix[2q];
m[[tgt,src]]=f;
m[[src+q,tgt+q]]=-antipode[f];
m
]
cNot[mat_,{tgt_,src_,f_}]:=cNotMatrix[Length[mat]/2,tgt,src,f].mat
ApplyCNOTSequence[mat_,srctgts_]:=Fold[(cNotMatrix[Dimensions[#1][[1]]/2,#2/.List->Sequence].#1)&,mat,srctgts]

(*drop rows or columns*)
factorOut[mat_,tgt_]:=
(*`tgt' can either be a number or a list of numbers*)
Fold[Delete,mat,Reverse[Sort[Flatten[{{tgt},{tgt+Length[mat]/2}}]]]]
dropColumn[mat_,cols_]:=
(*`cols' can either be a number or a list of numbers*)
Transpose[Fold[Delete,Transpose[mat],Reverse[Sort[Flatten[{cols}]]]]]

Twist Phase Tools

(*calculate points where two operators anticommute*)
SymMat[n_]:=KroneckerProduct[{{0,1},{1,0}},IdentityMatrix[n]]
ErrorSyndrome[vecs1_,vecs2_]:=Transpose[vecs1].SymMat[Dimensions[vecs1][[1]]/2].Map[antipode,vecs2,{2}]

(*tools for simplifying and truncating fractal symmetries*)
TruncateYs[poly_,maxN_]:=ExpandAll[poly]/.{y^n_->If[n>maxN||n<-1*maxN,0,y^n]}
ImposePBC2[poly_,n_]:=ExpandAll[poly]/.{g[1/x]^c_ /; c≥(2*n)->g[1/x]^(c-2*n)*(h[1/x]^(2*n)/x^n+1),g[x]^c_ /; c≥(2*n)->g[x]^(c-2*n)*(h[x]^(2*n)*x^n+1)}
SimplifyRule1[expr_,f_,l_]:=ExpandAll[ExpandAll[expr]/.{S->Sum[y^i*(f)^i,{i,0,l-1}],Q->(Sum[(y)^(Mod[i,l])*(antipode[f])^(Mod[l-i,l]),{i,0,l-1}])}]/.{y^c_->y^(Mod[c,l])}
SimplifyRule3[expr_]:=ExpandAll[expr]/.{S*y*g[x]^2->S*(y*x*h[x]^2+1),S*y^(c_)*g[x]^(k_)/; BitAnd[c,c-1]== 0 &&k==2*c->S*(1+(x*y)^c*h[x]^(k)),Q*g[1/x]^2/y->Q*(h[1/x]^2/(x*y)+1),Q*g[1/x]^(k_)*(1/y)^c_/;BitAnd[c,c-1] == 0 &&k==2*c->Q*(1+h[1/x]^(k)*(1/(x*y))^c)}

(*tools for simplifying and truncating linear symmetries*)
SimplifyRuleLin[expr_,l_]:=expr/.{Xs->Sum[x^i,{i,0,l-1}],Ys->Sum[y^i,{i,0,l-1}]}
SimplifyRuleLin2[expr_]:=ExpandAll[expr]/.{Xs*x^n_->Xs,Xs*x->Xs,Ys*y^n_->Ys,Ys*y->Ys}
SimplifyRuleLin3[expr_,l_]:=ExpandAll[expr]/.{c_^n_/;n==l-1->1/c}
TruncateDir[poly_,var_,maxN_]:=ExpandAll[poly]/.{var^n_/;(n>maxN||n<-1*maxN)->0}
PBCDir[poly_,var_,l_]:=ExpandAll[poly]/.{var^n_/;(n<0)->var^(l+n)}/.{var^n_/;(n==l)->1}

(*given a set of error syndromes of a truncated symmetry (points where it anticommutes), calculate an equivalent operator that only acts on one edge*)
GetUVecs[syndromes_,nSyms_,nStabs_]:=Join[Join[SymMat[nStabs/2].Transpose[syndromes][[;;,1;;nSyms/2]],ConstantArray[0,{nStabs,nSyms/2}]],
Join[ConstantArray[0,{nStabs,nSyms/2}],SymMat[nStabs/2].Transpose[syndromes][[;;,nSyms/2+1;;nSyms]]],2]

(*takes a 2nxn matrix with 1s on diagonal and adds a copy of some other model to srcqubits and tgtqubits. Used for self-bifurcation*)
AddModel[model_,srcQubits_,tgtQubits_,mat_]:=With[{h=Dimensions[mat][[1]],w=Dimensions[mat][[2]],modelQ=Dimensions[model][[1]]/4,arr=Table[{srcQubits[[i]],tgtQubits[[j]],model[[j,i]]},{i,1,Length[srcQubits]},{j,1,Length[srcQubits]}]},
ReplacePart[mat,
Flatten[Map[With[{yq=#[[2]],xq=#[[1]]-h/4,elem=#[[3]]},
{{yq,xq}->elem,{3*h/4+xq,w/2+yq}->antipode[elem]}]&,arr,{2}]]
]
]

Local Gates Tools

(*Make local gates for bifurcating the B model. These correspond to cnots joining two rectangles in each sublattice.*)
GetQubitsSquare[w1_,h1_,xd1_,yd1_,w2_,h2_,xd2_,yd2_]:=Module[{mat1,matrows1,mat2,matrows2,nonzeros1,nonzeros2},
mat1=Transpose[(x^xd1*y^yd1)*{{1,0,0,0},{x^w1,0,0,0},{y^h1,0,0,0},{x^w1*y^h1,0,0,0}}]~coarseGrain~{x,2}~coarseGrain~{y,2}~coarseGrain~{x,2}~coarseGrain~{y,2};
mat2=Transpose[(x^xd2*y^yd2)*{{0,1,0,0},{0,x^w2,0,0},{0,y^h2,0,0},{0,x^w2*y^h2,0,0}}]~coarseGrain~{x,2}~coarseGrain~{y,2}~coarseGrain~{x,2}~coarseGrain~{y,2};
matrows1=mat1[[;;,{1,17,33,49}]];
matrows2=mat2[[;;,{1,17,33,49}]];
nonzeros1 = Position[matrows1,c_/;!(c===0),{2},Heads->False];
nonzeros2 = Position[matrows2,c_/;!(c===0),{2},Heads->False];
Join[
{matrows1[[#[[1]],#[[2]]]],#[[1]]}&/@nonzeros1,
{matrows2[[#[[1]],#[[2]]]],#[[1]]}&/@nonzeros2
]]
SquareGatesList[w1_,h1_,xd1_,yd1_,w2_,h2_,xd2_,yd2_]:=With[{
qubits =GetQubitsSquare[w1,h1,xd1,yd1,w2,h2,xd2,yd2]},
Flatten[Table[With[{tgt=qubits[[i]],src=qubits[[4+j]]},{tgt[[2]],src[[2]],tgt[[1]]/src[[1]]}],{i,1,4},{j,1,4}],1]
]
SquareGatesMatrix[w1_,h1_,xd1_,yd1_,w2_,h2_,xd2_,yd2_]:=Module[{gatesList,gatesListDropped},
gatesList=SquareGatesList[w1,h1,xd1,yd1,w2,h2,xd2,yd2];
gatesListDropped = gatesList/.{c_/;(c>8 && c≤24) ->c-4,c_/;(c>28)->c-8};
If[MemberQ[gatesList,c_/;((c>4&&c≤8)||(c>24&&c≤28)),2],
{},ApplyCNOTSequence[IdentityMatrix[48],gatesListDropped]
]]
ApplySquareGateSequence[mat_,squares_]:=Fold[(SquareGatesMatrix[#2/.List->Sequence].#1)&,mat,squares]

(*make local gates for the bifurcation of linear cluster state. Similar to above but on a smaller lattice and only 1x1 squares.*)
GetQubitsSquareSmall[xd1_,yd1_,xd2_,yd2_]:=Module[{mat1,matrows1,mat2,matrows2,nonzeros1,nonzeros2},
mat1=Transpose[(x^xd1*y^yd1)*{{1,0,0,0},{x,0,0,0},{y,0,0,0},{x*y,0,0,0}}]~coarseGrain~{x,2}~coarseGrain~{y,2};
mat2=Transpose[(x^xd2*y^yd2)*{{0,1,0,0},{0,x,0,0},{0,y,0,0},{0,x*y,0,0}}]~coarseGrain~{x,2}~coarseGrain~{y,2};
matrows1=mat1[[;;,{1,5,9,13}]];
matrows2=mat2[[;;,{1,5,9,13}]];
nonzeros1 = Position[matrows1,c_/;!(c===0),{2},Heads->False];
nonzeros2 = Position[matrows2,c_/;!(c===0),{2},Heads->False];
Join[
{matrows1[[#[[1]],#[[2]]]],#[[1]]}&/@nonzeros1,
{matrows2[[#[[1]],#[[2]]]],#[[1]]}&/@nonzeros2
]]
SquareGatesListSmall[xd1_,yd1_,xd2_,yd2_]:=With[{
qubits =GetQubitsSquareSmall[xd1,yd1,xd2,yd2]},
Flatten[Table[With[{tgt=qubits[[i]],src=qubits[[4+j]]},{tgt[[2]],src[[2]],tgt[[1]]/src[[1]]}],{i,1,4},{j,1,4}],1]
]
SquareGatesMatrixSmall[xd1_,yd1_,xd2_,yd2_]:=ApplyCNOTSequence[IdentityMatrix[16],SquareGatesListSmall[xd1,yd1,xd2,yd2]]
ApplySquareGateSequenceSmall[mat_,params_]:=Fold[((SquareGatesMatrixSmall@@#2).#1)&,mat,params]

(*make local gates for the Sierpinski fractal f(x)=1+x+1/x. If you want to change to a different fractal, simply change localshape*)
localShape={1,x,x^2,x*y};
GetQubitsShape[xd1_,yd1_,xd2_,yd2_]:=
Module[{mat1,mat2,matrows1,matrows2,nonzeros1,nonzeros2},
mat1=Transpose[(x^xd1*y^yd1*{antipode[#],0,0,0})&/@localShape]~coarseGrain~{y,2}~coarseGrain~{x,2}~coarseGrain~{y,2}~coarseGrain~{x,2};
mat2 = Transpose[(x^xd2*y^yd2*{0,#,0,0})&/@localShape]~coarseGrain~{y,2}~coarseGrain~{x,2}~coarseGrain~{y,2}~coarseGrain~{x,2};
matrows1=mat1[[;;,{1,17,33,49}]];
matrows2=mat2[[;;,{1,17,33,49}]];
nonzeros1 = Position[matrows1,c_/;!(c===0),{2},Heads->False];
nonzeros2 = Position[matrows2,c_/;!(c===0),{2},Heads->False];
Join[
{matrows1[[#[[1]],#[[2]]]],#[[1]]}&/@nonzeros1,
{matrows2[[#[[1]],#[[2]]]],#[[1]]}&/@nonzeros2
]
]
ShapeGatesList[xd1_,yd1_,xd2_,yd2_]:=With[{qubits=GetQubitsShape[xd1,yd1,xd2,yd2]},
With[{src=#[[2]],tgt=#[[1]]},{tgt[[2]],src[[2]],tgt[[1]]/src[[1]]}]&/@Tuples[{qubits[[1;;4]],qubits[[5;;8]]}]
]
ShapeGatesMatrix[xd1_,yd1_,xd2_,yd2_]:=ApplyCNOTSequence[IdentityMatrix[64],ShapeGatesList[xd1,yd1,xd2,yd2]]
ApplyShapeGateSequence[mat_,params_]:=Fold[((ShapeGatesMatrix@@#2).#1)&,mat,params]

2.Twist Phases for Fractals

l=0 General RG

(*coarse graining. We write f(x)=g(x^2)+x*h(x^2)*)
stabsCG=coarseGrain[coarseGrain[{{(u[x]+x*v[x])*(1+y*(g[x]+x*h[x])),0},{1,0},{0,1},{0,(u[1/x]+v[1/x]/x)*(1+(g[1/x]+h[1/x]/x)/y)}},{y,2}],{x,2}]

(*Symmetries. S=F(x,y) and Q=F(1/x,1/y).*)
symsCG = coarseGrain[dropColumn[coarseGrain[{{0,0},{S,0},{0,Q},{0,0}}/.{S->S*(1+y*(g[x]+x*h[x])),Q->Q*(1+(1/y)(g[1/x]+h[1/x]/x))},{y,2}],{2,4}],{x,2}]

(*calculating twist phases*)
truncatedSyms=SimplifyRule1[symsCG,g[x]^2+x*h[x]^2,8];
syndromes = ImposePBC2[TruncateYs[ErrorSyndrome[truncatedSyms,stabsCG],4],8];
SimplifyRule3[ErrorSyndrome[symsCG,GetUVecs[syndromes,4,8]]]//display2

(*creating the bifurcating models. Adding 4 models with w[x]=u,u,v,xv*)
GeneralStabs[w_,f_]:={{w*(1+y*f),0},{1,0},{0,1},{0,antipode[w]*(1+antipode[f]/y)}}
baseMat=Join[ConstantArray[0,{4,8}],IdentityMatrix[8],ConstantArray[0,{4,8}]];
bifurcatedStabs=AddModel[GeneralStabs[x*v[x],g[x]^2+x*h[x]^2],{8},{3},
AddModel[GeneralStabs[v[x],g[x]^2+x*h[x]^2],{7},{4},
AddModel[GeneralStabs[u[x],g[x]^2+x*h[x]^2],{6},{2},
AddModel[GeneralStabs[u[x],g[x]^2+x*h[x]^2],{5},{1},
baseMat]]]]

(*twist phases of bifurcating models*)

syndromes2 = ImposePBC2[TruncateYs[ErrorSyndrome[truncatedSyms,bifurcatedStabs],4],8];
SimplifyRule3[ErrorSyndrome[symsCG,GetUVecs[syndromes2,4,8]]]//display2

l=1 RG For (Subscript[H, A])

HA={{(1+y*(g[x]^2+x*h[x]^2)),0},{1,0},{0,1},{0,(1+(g[1/x]^2+h[1/x]^2/x)/y)}};
Hp={{x*(1+y*(g[x]^2+x*h[x]^2)),0},{1,0},{0,1},{0,(1/x)*(1+(g[1/x]^2+h[1/x]^2/x)/y)}};
Hm={{(1/x)*(1+y*(g[x]^2+x*h[x]^2)),0},{1,0},{0,1},{0,x*(1+(g[1/x]^2+h[1/x]^2/x)/y)}};

(*HA but with l=1*)
stabsHAyCG=coarseGrain[coarseGrain[{{y*(1+y*(g[x]+x*h[x])),0},{1,0},{0,1},{0,(1+(g[1/x]+h[1/x]/x)/y)/y}},{y,2}],{x,2}]

(*using a different basis for the symmetries is more convenient*)

symsCG2 = coarseGrain[dropColumn[coarseGrain[{{0,0},{S,0},{0,Q*y},{0,0}}/.{S->S*(1+y*(g[x]+x*h[x])),Q->Q*(1+(1/y)(g[1/x]+h[1/x]/x))},{y,2}],{2,4}],{x,2}]

(*calculating twist phases*)

truncatedSyms2=SimplifyRule1[symsCG2,g[x]^2+x*h[x]^2,8];
syndromes3 = ImposePBC2[TruncateYs[ErrorSyndrome[truncatedSyms2,stabsHAyCG],4],8];
SimplifyRule3[ErrorSyndrome[symsCG2,GetUVecs[syndromes3,4,8]]]//display2

(*making the bifurcating copy with 2 copies of HA (but on different qubits than we did in the first example)*)
bifurcatedStabs2=
AddModel[HA,{5},{3},
AddModel[HA,{6},{4},
baseMat]]

(*twist phases of bifurcating model*)

syndromes4 = ImposePBC2[TruncateYs[ErrorSyndrome[truncatedSyms2,bifurcatedStabs2],4],8];
SimplifyRule3[ErrorSyndrome[symsCG2,GetUVecs[syndromes4,4,8]]]//display2

l=1 RG For (SubPlus[H])

(*(SubPlus[H]) but with l=1*)
stabsHpyCG=coarseGrain[coarseGrain[{{x*y*(1+y*(g[x]+x*h[x])),0},{1,0},{0,1},{0,(1/(x*y))*(1+(g[1/x]+h[1/x]/x)/y)}},{y,2}],{x,2}]

(*calculating twist phases*)
syndromes5 = ImposePBC2[TruncateYs[ErrorSyndrome[truncatedSyms2,stabsHpyCG],4],8];
SimplifyRule3[ErrorSyndrome[symsCG2,GetUVecs[syndromes5,4,8]]]//display2

(*creating bifurcating model with (Subscript[H, A]) and (SubPlus[H])*)
bifurcatedStabs3=AddModel[Hp,{6},{3},
AddModel[HA,{5},{4},
baseMat]]

(*calculating twist phases*)
syndromes6 = ImposePBC2[TruncateYs[ErrorSyndrome[truncatedSyms2,bifurcatedStabs3],4],8];
SimplifyRule3[ErrorSyndrome[symsCG2,GetUVecs[syndromes6,4,8]]]//display2

l=-1 RG For (Subscript[H, A])

(*HA but with l=-1*)
stabsHAyCG2=coarseGrain[coarseGrain[{{(1+y*(g[x]+x*h[x]))/y,0},{1,0},{0,1},{0,y*(1+(g[1/x]+h[1/x]/x)/y)}},{y,2}],{x,2}]

(*using a different basis for the symmetries is more convenient*)

symsCG3 = coarseGrain[dropColumn[coarseGrain[{{0,0},{S,0},{0,Q/y},{0,0}}/.{S->S*(1+y*(g[x]+x*h[x])),Q->Q*(1+(1/y)(g[1/x]+h[1/x]/x))},{y,2}],{2,4}],{x,2}]

(*calculating twist phases*)

truncatedSyms3=SimplifyRule1[symsCG3,g[x]^2+x*h[x]^2,8];
syndromes7 = ImposePBC2[TruncateYs[ErrorSyndrome[truncatedSyms3,stabsHAyCG2],4],8];
SimplifyRule3[ErrorSyndrome[symsCG3,GetUVecs[syndromes7,4,8]]]//display2

(*making the bifurcating copy with 2 copies of HA *)
bifurcatedStabs4=
AddModel[HA,{7},{1},
AddModel[HA,{8},{2},
baseMat]]

(*twist phases of bifurcating model*)

syndromes8 = ImposePBC2[TruncateYs[ErrorSyndrome[truncatedSyms3,bifurcatedStabs4],4],8];
SimplifyRule3[ErrorSyndrome[symsCG3,GetUVecs[syndromes8,4,8]]]//display2

l=-1 RG For (SubPlus[H])

(*(SubPlus[H]) but with l=1*)
stabsHpyCG2=coarseGrain[coarseGrain[{{x*(1+y*(g[x]+x*h[x]))/y,0},{1,0},{0,1},{0,(y/x)*(1+(g[1/x]+h[1/x]/x)/y)}},{y,2}],{x,2}]

(*calculating twist phases*)
syndromes9 = ImposePBC2[TruncateYs[ErrorSyndrome[truncatedSyms3,stabsHpyCG2],4],8];
SimplifyRule3[ErrorSyndrome[symsCG3,GetUVecs[syndromes9,4,8]]]//display2

(*creating bifurcating model with (Subscript[H, A]) and (SubPlus[H])*)
bifurcatedStabs5=AddModel[Hp,{8},{1},
AddModel[HA,{7},{2},
baseMat]]

(*calculating twist phases*)
syndromes10 = ImposePBC2[TruncateYs[ErrorSyndrome[truncatedSyms3,bifurcatedStabs5],4],8];
SimplifyRule3[ErrorSyndrome[symsCG3,GetUVecs[syndromes10,4,8]]]//display2

3.Twist Phases for Linear Model

Bifurcation of (Subscript[H, lin])

linStabs={{1+x+1/y+x/y,0},{1,0},{0,1},{0,1+y+1/x+y/x}}

bStabs=Transpose[{{0,1+1/x,1+1/x,1,0,0,0,0,0,0,0,0},{1+y,1/x+y,1+1/x,0,1,0,0,0,0,0,0,0},{1+y,1+y,0,0,0,1,0,0,0,0,0,0},{0,0,0,0,0,0,1,0,0,0,1+1/y,1+1/y},{0,0,0,0,0,0,0,1,0,1+x,x+1/y,1+1/y},{0,0,0,0,0,0,0,0,1,1+x,1+x,0}}]

(*coarse graining original model*)
linStabsCG =coarseGrain[coarseGrain[linStabs,{x,2}],{y,2}]

(*coarse graining symmetries. Ys=1+y+y^2+... Xs=1+x+x^2+...*)
symsLinCG=dropColumn[coarseGrain[coarseGrain[{{0,0,0,0},{Xs,Ys,0,0},{0,0,Xs,Ys},{0,0,0,0}}/.{Xs->(1+x)*Xs,Ys->(1+y)*Ys},{x,2}],{y,2}],{3,4,6,8,11,12,14,16}]

(*calculating twist phases*)
symsLinTruncated = SimplifyRuleLin[symsLinCG,8];
syndromes11=ErrorSyndrome[symsLinTruncated,linStabsCG];
syndromes11Y=PBCDir[TruncateDir[syndromes11,y,4],x,8];
syndromes11X=PBCDir[TruncateDir[syndromes11,x,4],y,8];

(*y twist phase*)
SimplifyRuleLin3[SimplifyRuleLin2[ErrorSyndrome[GetUVecs[syndromes11Y,8,8],symsLinCG]],8]//display2

(*x twist phase*)
SimplifyRuleLin3[SimplifyRuleLin2[ErrorSyndrome[GetUVecs[syndromes11X,8,8],symsLinCG]],8]//display2

(*making the bifurcated model*)
bifurcatedStabsLin=AddModel[linStabs,{7},{2},
AddModel[bStabs,{5,6,8},{1,3,4},
baseMat]]

(*calculating twist phase*)
syndromes12=ErrorSyndrome[symsLinTruncated,bifurcatedStabsLin];
syndromes12Y=PBCDir[TruncateDir[syndromes12,y,4],x,8];
syndromes12X=PBCDir[TruncateDir[syndromes12,x,4],y,8];

(*y twist phase*)
SimplifyRuleLin3[SimplifyRuleLin2[ErrorSyndrome[GetUVecs[syndromes12Y,8,8],symsLinCG]],8]//display2

(*x twist phase*)
SimplifyRuleLin3[SimplifyRuleLin2[ErrorSyndrome[GetUVecs[syndromes12X,8,8],symsLinCG]],8]//display2

Bifurcation of (Subscript[H, B])

(*coarse graining symmetries, but just on the B model qubits*)

symsLinCG2=dropColumn[coarseGrain[coarseGrain[factorOut[symsLinCG,{2,7}]/.{Xs->(1+x)*Xs,Ys->(1+y)*Ys},{x,2}],{y,2}],{3,4,7,8,10,12,14,16,19,20,23,24,26,28,30,32}]

(*coarse graining stabilizers*)
bStabsCG=coarseGrain[coarseGrain[bStabs,{x,2}],{y,2}]

(*calculating twist phases*)
symsLinTruncated2= SimplifyRuleLin[symsLinCG2,8];
syndromes13=ErrorSyndrome[symsLinTruncated2,bStabsCG];
syndromes13Y=PBCDir[TruncateDir[syndromes13,y,4],x,8];
syndromes13X=PBCDir[TruncateDir[syndromes13,x,4],y,8];

(*y twist phases*)
SimplifyRuleLin2[ErrorSyndrome[GetUVecs[syndromes13Y,16,24],symsLinCG2]]//display2

(*x twist phases*)
SimplifyRuleLin2[ErrorSyndrome[GetUVecs[syndromes13X,16,24],symsLinCG2]]//display2

(*Now building bifurcated model*)

baseMat2=Join[ConstantArray[0,{12,24}],IdentityMatrix[24],ConstantArray[0,{12,24}]];
bifurcatedStabsLin2=AddModel[bStabs,{15,19,23},{3,7,11},AddModel[bStabs,{14,18,22},{2,6,10},baseMat2]]

(*calculating twist phases*)
syndromes14=ErrorSyndrome[symsLinTruncated2,bifurcatedStabsLin2];
syndromes14Y=PBCDir[TruncateDir[syndromes14,y,4],x,8];
syndromes14X=PBCDir[TruncateDir[syndromes14,x,4],y,8];

(*y twist phases*)
SimplifyRuleLin2[ErrorSyndrome[GetUVecs[syndromes14Y,16,24],symsLinCG2]]//display2

(*x twist phases*)
SimplifyRuleLin2[ErrorSyndrome[GetUVecs[syndromes14X,16,24],symsLinCG2]]//display2

4.Local Gates

In this section, we take the gates used for bifurcation, and use local gates to turn these into identity. Since these gates are made up of commuting CNOTs, this means they must themselves be local.

Bifurcation of (Subscript[H, A])

This is an example for the Sierpinski fractal, f(x)=1+x+1/x. But a similar algorithm may be used for any fractal, just by changing localshape (see tools).

(*coarse graining syms again and substitute f(x)=1+x+1/x*)
symsCG4=dropColumn[coarseGrain[coarseGrain[(symsCG/.{S->S*(1+y*(h[x]+x*g[x])),Q->Q*(1+(h[1/x]+g[1/x]/x)/y)})/.{g[c_]->1,h[c_]->1+1/c},{y,2}],{x,2}],{2,4,6,8,10,12,14,16}];
symsCG4//display2

(*same for stabilizers*)
stabsCG2=(HA/.{g[c_]->1,h[c_]->1+1/c})~coarseGrain~{y,2}~coarseGrain~{x,2}~coarseGrain~{y,2}~coarseGrain~{x,2};
stabsCG2//display2

(*coarse graining bifurcated stabilizers and substituting 1+x+1/x*)

bifurcatedStabsCG=coarseGrain[coarseGrain[AddModel[HA,{6},{2},AddModel[HA,{5},{1},baseMat]]/.{g[c_]->1,h[c_]->1+1/c},{y,2}],{x,2}];
bifurcatedStabsCG//display2

(*gates needed to take original model to bifurcated model*)

padLeft = Join[IdentityMatrix[16],ConstantArray[0,{48,16}]];
padRight =  Join[ConstantArray[0,{48,16}],IdentityMatrix[16]];
totalGates=Join[padLeft,bifurcatedStabsCG,padRight,2].Join[padLeft,stabsCG2,padRight,2];
totalGates//display2

(*Shows that 8 sets of local gates make this identity*)

gateSequence={{2,0,0,-2},{3,0,1,-2},{4,0,2,-2},{5,0,3,-2},{2,2,0,0},{3,2,1,0},{4,2,2,0},{5,2,3,0}}

(*Each row is one local unitary operator. In the language of appendix C, the first two entries in each gate are the position of PA, the second two entries are position of PB. We can write this using c(x,y) and d(x,y) as in the appendix. For example, in the first gate, c=1, d=y^-2, in the second gate, c=x, d=x*y^-2, in the third, c=x^2, d=x^2*y^-2, etc.*)

ApplyShapeGateSequence[totalGates,gateSequence]//display2

Bifurcation of (Subscript[H, B])

(*Gates required to take HB to two copies of itself*)
padLeft = Join[IdentityMatrix[12],ConstantArray[0,{36,12}]];
padRight =  Join[ConstantArray[0,{36,12}],IdentityMatrix[12]];
totalGates2=Join[padLeft,bStabsCG,padRight,2].Join[padLeft,bifurcatedStabsLin2,padRight,2];
totalGates2//display2

(*Show that 4 sets of local gates make this identity*)

gateSequence2={{-2,1,1,0,2,1,0,0},{2,1,1,2,2,1,2,2},{1,2,0,0,1,2,0,-1},{1,2,2,2,1,2,2,1}}

ApplySquareGateSequence[totalGates2,gateSequence2]//display2

Bifurcation of (Subscript[H, lin])

(*Gates required to take Hlin to a copy of HB and Hlin*)
padLeft = Join[IdentityMatrix[4],ConstantArray[0,{12,4}]];
padRight =  Join[ConstantArray[0,{12,4}],IdentityMatrix[4]];
totalGates3=Join[padLeft,linStabsCG,padRight,2].Join[padLeft,bifurcatedStabsLin,padRight,2];
totalGates3//display2

(*4 copies of a 1x1 square gate make this identity*)
gateSequence3={{0,0,0,0},{1,0,1,0},{0,1,0,1}}

ApplySquareGateSequenceSmall[totalGates3,gateSequence3]//display2