Majorana operators in spin language

This is a simplified version of a small subset of calculations for our paper:
https://arxiv.org/abs/2412.06377

Here I implement the Majorana operators (SuperDagger[η]), η as (2^L) x (2^L) matrices, following e.g. P. Coleman’s book, or the review by Mbeng, Russomanno, and Santoro 
( https://arxiv.org/abs/2009.09208 ). Of course, using big (2^L) x (2^L) matrices is not efficient for long chains. But here I focus on checking basic statements using minimal assumptions.

Load the Quantum Framework:

Remove["Global`*"]//Quiet
PacletInstall["Wolfram/QuantumFramework"];
Needs["Wolfram`QuantumFramework`"];
Needs["Notation`"] //Quiet
Notation[(NotationTemplateTag[(Subscript[J, 0])]) ⟺ (NotationTemplateTag[J0])] //Quiet;

Consider 4 sites.
(This file runs quickly also for L = 5, but already for L = 6 it starts taking a while. This is because I check e.g. commutation relation of (2^6) x (2^6)  = 64 x 64 matrices, but do it (L^3) = 216 times. So for longer chains either only do spot checks, or better, don’t use the spin representation. )

L=4;

The Hamiltonian operator in bra-ket language is:

-J QuantumOperator["X",{1,2}]-J0 QuantumOperator["X",{2,3}]-h QuantumOperator["X",{3,4}]-h QuantumOperator["Z",{2}]-J0 QuantumOperator["Z",{3}]-J QuantumOperator["Z",{4}]

(It is easy to automate this, but it’s easier to read if typed in.)

And here is H as a matrix:

H = %["Matrix"]//Normal

The Ising operator Q is the tensor product of Z operator on all sites:

QuantumOperator["Z",{1,2,3,4}]
Q = %["Matrix"]//Normal

The operator Q commutes with the Hamiltonian:

Commutator[A_,B_] := A.B-B.A;
ZeroMatrix = ConstantArray[0,{2^L,2^L}];

Commutator[H,Q]==ZeroMatrix

Id = QuantumOperator["I"[Table[2,{j,1,L}]]];
δH = %["Matrix"]//Normal;

(“δH” stands for “identity operator in Hilbert space”, to avoid confusion with two other identity operators introduced below.)

The single-site hardcore bosons (SuperDagger[b]), b (se references above, but not necessary, these are just steps on the way to the fermions):

bdSite = 1/2(QuantumOperator["X"] -I QuantumOperator["Y"])
bSite = 1/2(QuantumOperator["X"] +I QuantumOperator["Y"])

The hardcore boson number operator:

nSite = QuantumOperator[Normal[bdSite["Matrix"]].Normal[bSite["Matrix"]]]

Assemble to length-L arrays (Subscript[b, HB]) and (Subscript[b, HB])† (with “HB” for Hardcore Boson,to avoid conflict with the index b below: these (Subscript[b, HB]) are not used again later), and convert to ordinary matrices:

bHBd = Table[Normal[(0Id +QuantumOperator[bdSite,{j}])["Matrix"]],{j,1,L}];
bHB =Table[Normal[(0Id +QuantumOperator[bSite ,{i}])["Matrix"]],{i,1,L}];
n =Table[Normal[(0Id +QuantumOperator[nSite,{i}])["Matrix"]],{i,1,L}];
Z =Table[Normal[(0Id +QuantumOperator["Z",{i}])["Matrix"]],{i,1,L}];

Define K-strings using Z  = (e^(π i (Subscript[n, j]))):

Table[Z[[j]]==MatrixExp[π I n[[j]]],{j,1,L}]

Here we can use Product for matrix product, since Z is diagonal (will not work below):

KString[j_]:=If[j>1,Product[Z[[l]],{l,1,j-1}],δH]

The string goes to the left. 
(In the draft we distinguish operators on spins and on Majorana, here they are the same.)
We make the (SuperDagger[c]), c fermions from (SuperDagger[b]), b Hardcore Bosons “HB” (the b’s won’t be used below):

Table[bHBd[[j]].KString[j],{j,2,L}];
cd = Join[{bHBd[[1]]},%];
Table[bHB[[j]].KString[j],{j,2,L}];
c  = Join[{bHB[[1]]},%];

Note that b and the K-string commute.
Here’s one example, (SuperDagger[c])4 :

cd[[4]]

The c-fermions satisfy the algebra { (SuperDagger[(Subscript[c, j])]), (Subscript[c, k]) } = (Subscript[δ, jk]):

δSite =  IdentityMatrix[L];
Table[cd[[j]].c[[k]]+ c[[k]].cd[[j]]==δSite[[j,k]] δH,{j,1,L},{k,1,L}]

(In particular we do not allow the flipped sign “–(Subscript[δ, jk])” on the RHS in the algebra.)

As expected, they are not Hermitian:

Table[ConjugateTranspose[c[[j]]]== c[[j]],{j,1,L}]

Construct the Majorana fermions (Subscript[η, odd]) = (SuperDagger[c])+ c, (Subscript[η, even]) = i((SuperDagger[c])– c), assembled into a single array indexed by a = 1 .. 2L (note that the MRS review does it differently, they introduce a two-index object instead of a single array):

η=Flatten[Table[{cd[[j]]+c[[j]],I(cd[[j]]-c[[j]])},{j,1,L}],1];

Check the fermion algebra { (Subscript[η, a]), (Subscript[η, b]) } = 2(Subscript[δ, ab]), where again a runs from 1 to 2L:

δM = IdentityMatrix[2L];Table[η[[a]].η[[b]]+η[[b]].η[[a]]==2δM[[a,b]]δH,{a,1,2L},{b,1,2L}]

The Majorana (Subscript[η, a]) are Hermitian:

Table[ConjugateTranspose[η[[a]]]==η[[a]],{a,1,2L}]

Check that ((Subscript[η, a])^2) = 1:

Table[η[[a]].η[[a]]==δH,{a,1,2L}]

Check the spin-Majorana mapping:
1. the L-1 Ising couplings  (Subscript[X, j]) (Subscript[X, j+1]) = –i (Subscript[η, 2j])  (Subscript[η, 2j+1]):

Table[Normal[(0 Id + QuantumOperator["X",{j,j+1}])["Matrix"]]==-I η[[2j]].η[[2j+1]],{j,1,L-1}]

Jvec={J,J0,h};
HJ =  Sum[Jvec[[j]] I η[[2j]].η[[2j+1]],{j,1,L-1}];

For example for L = 4, (Subscript[H, J]) = J (i (Subscript[η, 2])(Subscript[η, 3]) + i (Subscript[η, 4])(Subscript[η, 5]) + i (Subscript[η, 6])(Subscript[η, 7])).
The patterns are:

Table[ArrayPlot[-I η[[2j]].η[[2j+1]]],{j,1,L-1}]

2. the L magnetic field operators (Subscript[Z, j]) = –i (Subscript[η, 2j -1])(Subscript[η, 2j ]):

Table[Normal[(0 Id + QuantumOperator["Z",{j}])["Matrix"]]==-I η[[2j-1]].η[[2j]],{j,1,L}]

For Z, I need to make –1 a different color (here blue) to see the difference:

Table[ArrayPlot[Normal[(0 Id + QuantumOperator["Z",{j}])["Matrix"]],ColorRules->{-1->Blue}],{j,1,L}]

hvec = {h,J0,J};
Hh =  Sum[hvec[[j-1]]I η[[2j-1]].η[[2j]],{j,2,L}];

For example for L = 4, we have (Subscript[H, h]) = h (i (Subscript[η, 3])(Subscript[η, 4]) + i (Subscript[η, 5])(Subscript[η, 6]) + i (Subscript[η, 7])(Subscript[η, 8])).

HJ + Hh == H

Define the Majorana swap gate, and its dagger:

swap = Table[1/Sqrt[2](δH +η[[a]].η[[b]]),{a,1,2L},{b,1,2L}];
swapd = Table[1/Sqrt[2](δH -η[[a]].η[[b]]),{a,1,2L},{b,1,2L}];

Check the action of these swap gates on a single η :

Table[If[a≠b,swap[[a,b]].η[[c]].swapd[[a,b]]==-δM[[a,c]] η[[b]]+δM[[b,c]] η[[a]]+(1-δM[[a,c]])(1-δM[[b,c]])η[[c]],{}],{c,1,2L},{b,1,2L},{a,1,2L}]

The a ≠ b restriction (where the “{ }” appears above) is in the paper.

The Kramers-Wannier operator (Subscript[U, KW]) 
First, note that the Mathematica command Product does not do matrix product in the sense of matrix multiplication, as you might have thought:

MyMatrix[1]=Array[a,{2,2}];
MyMatrix[2]=Array[b,{2,2}];
Product[MyMatrix[i],{i,1,2}]

One way to make an actual matrix product is to construct a list of matrices, and Fold them together using the actual matrix multiplication “Dot”, to make (Subscript[U, KW]):

UKW =Exp[-2 π I L/8] Fold[Dot,Table[swap[[j,j+1]],{j,2L-1,1,-1}]]
UKWd = Exp[2 π I L/8] Fold[Dot,Table[swapd[[j,j+1]],{j,1,2L-1}]];

It’s not so easy to read the above, so check by hand for L = 4:

{UKW ==Exp[-2 π I L/8]swap[[7,8]].swap[[6,7]].swap[[5,6]].swap[[4,5]].swap[[3,4]].swap[[2,3]].swap[[1,2]],
UKWd ==Exp[2 π I L/8]swapd[[1,2]].swapd[[2,3]].swapd[[3,4]].swapd[[4,5]].swapd[[5,6]].swapd[[6,7]].swapd[[7,8]]}

OK.

The action of (Subscript[U, KW]) on η for a > 1 is:

Table[UKW.η[[a]].UKWd == η[[a-1]],{a,2,2L}]

... and the action of (Subscript[U, KW]) on η for the exceptional case a = 1 is:

UKW.η[[1]].UKWd ==-η[[2L]]

OK.

We can double check the action on the spin operators in this language (“double check”, because here it is actually equivalent to the above):
1. (Subscript[Z, j]) maps to (Subscript[X, j-1]) (Subscript[X, j]) under Kramers-Wannier:

Table[UKW.η[[2j-1]].η[[2j]].UKWd == η[[2j-2]].η[[2j-1]],{j,2,L}]

1. (Subscript[X, j]) (Subscript[X, j+1]) maps to (Subscript[Z, j]) under Kramers-Wannier:

Table[UKW.η[[2j]].η[[2j+1]].UKWd == η[[2j-1]].η[[2j]],{j,1,L-1}]

OK. 

We already checked the action on the single (Subscript[η, 1]) above, that in the spin language includes an Ising charge Q. (In the draft we distinguish the notation for the action on spin operators and on the Majorana operators, the latter is curly, but here I do not.)

The reflection operator R 
The claim is that it acts as i , Q and site reflection (note the index on η runs from 1 to 2L).
It acts on spins as SWAP gates (which by definition squares to identity):

swapSite[i_,j_] := Normal[(0 Id+ QuantumOperator["SWAP",{i,j}])["Matrix"]];

Like for (Subscript[U, KW]) above, I construct both a generic R by Fold/Dot and check it explicitly for L = 4:

R =Fold[Dot,Table[swapSite[j,L-j+1],{j,1,L/2}]];
R== swapSite[1,4].swapSite[2,3]

I checked above that fermion parity P is given by Q in the spin language:

Q = Product[Z[[l]],{l,1,L}];

Now I can act with R on the Majoranas (Subscript[η, a]), separately for odd and even a:

hat[a_] := 2L+1-a;Table[R.η[[2j-1]].Inverse[R]==I Q.η[[hat[2j-1]]],{j,1,L}]

Table[R.η[[2j]].Inverse[R]==-I Q.η[[hat[2j]]],{j,1,L}]

R squares to +1 on the Majoranas:

Table[R.R.η[[a]].Inverse[R].Inverse[R]==η[[a]],{a,1,2L}]

R obviously does nothing to the algebra { (Subscript[η, a]), (Subscript[η, b]) } = 2(Subscript[δ, ab]) because of the identity operator.

1. R acting on the spin operator (Subscript[X, j]) in the η representation, i.e. after Jordan-Wigner:

Table[R.KString[j].η[[2j-1]].Inverse[R]==KString[L-j+1].η[[2(L+1-j)-1]],{j,1,L}]

1. R acting on the spin operator (Subscript[Z, j]) in the η representation, i.e. after Jordan-Wigner:

Table[(-I)R.η[[2j-1]].η[[2j]].Inverse[R]==(-I)η[[2(L+1-j)-1]].η[[2(L+1-j)]],{j,1,L}]

The composite operator S = R (Subscript[U, KW])
Finally our operator:

S = R.UKW;

It acts on η as, again dividing into odd and even cases,
and at first excluding a = 1 (which is j = 1):

t[a_] := 2L+2-a;
Table[S.η[[2j-1]].Inverse[S]==-I Q.η[[t[2j-1]]],{j,2,L}]

Table[S.η[[2j]].Inverse[S]==I Q.η[[t[2j]]],{j,1,L}]

We can represent the action of S more compactly if we set:

η = Join[η[[1;;2L]],{-η[[1]]}];
η[[2L+1]] == -η[[1]]

Then we can extend the S action on odd (Subscript[η, a]) above to include also a = 1 ( j = 1) :

With this, double check the known action of S on η, that reflects Majoranas to the tilded ones:

Table[S.η[[2j-1]].Inverse[S]==-I Q.η[[t[2j-1]]],{j,1,L}]

I already check (S^2) = Q in another file, but let’s do it again for completeness:

S.S==Q

OK. The (Subscript[ℤ, 4]) symmetry S squares to the  (Subscript[ℤ, 2]) Ising symmetry Q.
(Whether S merits the name symmetry depends on your definition, another name is “quadrality”.)

Here is the left-right combination plus:

ξp =1/Sqrt[2]Table[η[[a]]+(-1)^a η[[2L+2-a]],{a,1,2L}];

𝒩 = J^2(1-(h/J)^L)/(J^2-h^2)

ηSZM = 1/𝒩 Sum[(h/J)^(L/2-j)ξp[[2j]],{j,1,L/2}];

The commutator of H and the strong zero mode is zero:

Commutator[H,ηSZM]==ZeroMatrix//Simplify

So [H,η] =0.