Definitions

The representation of the lattice and Hilbert space

The basis of the Hilbert space

nϕ=(2^3);
ϕMaxval=(Sqrt[(π/2)])((nϕ-1)/(Sqrt[nϕ]))//N;
δϕ[ϕMax_]:=((2ϕMax)/(nϕ-1));
ϕval[ϕMax_,iϕ_]:=-ϕMax+iϕ δϕ[ϕMax]
basis[iϕ_]:=Table[If[i==iϕ,1,0],{i,0,nϕ-1}]

The Fourier transform and canonical momenta

This transformation goes between the k states and the ϕ states, and is the symmetric Fourier transform that Martin mentions in the digitization paper

FTop=(1/(Sqrt[nϕ]))Table[Exp[((2 I π )/(nϕ))j k],{j,-((nϕ-1)/2),((nϕ-1)/2)},{k,((-nϕ)/2),((nϕ)/2)-1}]//N;
FTInvop=ConjugateTranspose[FTop];

πMax[ϕMax_]:=(π/(δϕ[ϕMax]));
δπ[ϕMax_]:=((2π)/(δϕ[ϕMax] nϕ));
πval[ϕMax_,βk_]:=-πMax[ϕMax]+(βk+(1/2)) δπ[ϕMax]

ListPlot[Table[{ϕval[ϕMaxval,iϕ],1},{iϕ,0,nϕ-1}],PlotRange->{{-4,4},{-.8,1.2}}]
ListPlot[Table[{πval[ϕMaxval,iϕ],1},{iϕ,0,nϕ-1}],PlotRange->{{-4,4},{-.8,1.2}}]

ListPlot[Table[{ϕval[1.2ϕMaxval,iϕ],1},{iϕ,0,nϕ-1}],PlotRange->{{-4,4},{-.8,1.2}}]
ListPlot[Table[{πval[1.2ϕMaxval,iϕ],1},{iϕ,0,nϕ-1}],PlotRange->{{-4,4},{-.8,1.2}}]

The operators

The operators involving powers of ϕ

ϕop[ϕMax_]:=DiagonalMatrix[Table[ϕval[ϕMax,iϕ],{iϕ,0,nϕ-1}]];
ϕ2op[ϕMax_]:=DiagonalMatrix[Table[((ϕval[ϕMax,iϕ])^2),{iϕ,0,nϕ-1}]];

The operators involving powers of π

Using finite differences (without the Fourier transform)

Define the derivatives

ival[iϕ_]=Mod[iϕ,nϕ];
deriv1[i_,iϕ_]=Piecewise[{{(1/2),i==ival[iϕ+1]},{((-1)/2),i==ival[iϕ-1]}},0];
deriv2[i_,iϕ_]=Piecewise[{{2,i==ival[iϕ]},{-1,i==ival[iϕ-1]},{-1,i==ival[iϕ+1]}},0];

Table[deriv1[iϕ,i],{i,0,nϕ-1},{iϕ,0,nϕ-1}]//MatrixForm
Table[deriv2[iϕ,i],{i,0,nϕ-1},{iϕ,0,nϕ-1}]//MatrixForm

The operators involving π

iπop[1][ϕMax_]:=(1/(δϕ[ϕMax]))Table[deriv1[iϕ,i],{i,0,nϕ-1},{iϕ,0,nϕ-1}];
π2op[1][ϕMax_]:=(1/((δϕ[ϕMax])^2))Table[deriv2[i,iϕ],{i,0,nϕ-1},{iϕ,0,nϕ-1}];

Using the Fourier transform of the finite difference operators

iπop[2][ϕMax_]:=FTInvop.DiagonalMatrix[Table[-(I/(δϕ[ϕMax]))Sin[πval[ϕMax,βk]δϕ[ϕMax]],{βk,0,nϕ-1}]].FTop;
π2op[2][ϕMax_]:=FTInvop.DiagonalMatrix[Table[(4/((δϕ[ϕMax])^2))((Sin[((πval[ϕMax,βk]δϕ[ϕMax])/2)])^2),{βk,0,nϕ-1}]].FTop;

((δϕ[ϕMaxval])^2)π2op[2][ϕMaxval]//Chop//MatrixForm

Using the Fourier transform using exact momentum operator

iπop[3][ϕMax_]:=FTInvop.DiagonalMatrix[Table[-I πval[ϕMax,βk],{βk,0,nϕ-1}]].FTop;
π2op[3][ϕMax_]:=FTInvop.DiagonalMatrix[Table[((πval[ϕMax,βk])^2),{βk,0,nϕ-1}]].FTop;

((δϕ[ϕMaxval])^2)π2op[3][ϕMaxval]//Chop//MatrixForm

The derived operators needed in the rest of the code

The Hamiltonian

HopV[1][ϕMax_]:=(1/2)(ϕ2op[ϕMax]+π2op[1][ϕMax]);
HopV[2][ϕMax_]:=(1/2)(ϕ2op[ϕMax]+π2op[2][ϕMax]);
HopV[3][ϕMax_]:=(1/2)(ϕ2op[ϕMax]+π2op[3][ϕMax]);

Creation and annihilation operators

adagopV[1][ϕMax_]:=(1/(Sqrt[2]))(ϕop[ϕMax]-iπop[1][ϕMax]);
aopV[1][ϕMax_]:=(1/(Sqrt[2]))(ϕop[ϕMax]+iπop[1][ϕMax]);
adagopV[2][ϕMax_]:=(1/(Sqrt[2]))(ϕop[ϕMax]-iπop[2][ϕMax]);
aopV[2][ϕMax_]:=(1/(Sqrt[2]))(ϕop[ϕMax]+iπop[2][ϕMax]);
adagopV[3][ϕMax_]:=(1/(Sqrt[2]))(ϕop[ϕMax]-iπop[3][ϕMax]);
aopV[3][ϕMax_]:=(1/(Sqrt[2]))(ϕop[ϕMax]+iπop[3][ϕMax]);

The time evolution. These are now always using the Fourier transform with the exact momentum

evolveHϕV[t_][ϕMax_]:=DiagonalMatrix[Table[Exp[-I( (((ϕval[ϕMax,i])^2)/2))t],{i,0,nϕ-1}]];
evolveHπV[t_][ϕMax_]:=FTInvop.DiagonalMatrix[Table[Exp[-I (((πval[ϕMax,i])^2)/2)t],{i,0,nϕ-1}]].FTop;

evolveV[t_,nn_][ϕMax_]:=MatrixPower[evolveHϕ[ϕMax][(t/(2nn))].evolveHπ[ϕMax][(t/nn)].evolveHϕ[ϕMax][(t/(2nn))],nn]

Eigenvalues[HopV[1][ϕMaxval]]

Eigenvalues[HopV[2][ϕMaxval]]

Eigenvalues[HopV[3][ϕMaxval]]

Create all the operators with fixed values of ϕMax so rest of code still works

ϕval[iϕ_]:=ϕval[ϕMaxval,iϕ]

The Hamiltonian

Hop[n_][λλ_]:=HopV[n][λλ][ϕMaxval]

Creation and annihilation operators

adagop[n_]:=adagopV[n][ϕMaxval]
aop[n_]:=aopV[n][ϕMaxval]

The time evolution. These are now always using the Fourier transform with the exact momentum

evolveHϕ[t_][λλ_]:=evolveHϕV[t][λλ][ϕMaxval]
evolveHπ[t_]:=evolveHπV[t][λλ][ϕMaxval]

evolve[t_,nn_][λλ_]:=evolveV[t,nn][λλ][ϕMaxval]

Check the commutator of the x,p

Clear[vac,comm]

vac[n_,pre_]:=Eigensystem[HopV[n][0][pre ϕMaxval]][[2,-1]]//Chop;
comm[n_,pre_]:=((ϕop[pre ϕMaxval].iπop[n][pre ϕMaxval]-iπop[n][pre ϕMaxval].ϕop[pre ϕMaxval]).vac[n,pre])[[nϕ/2]]/vac[n,pre][[nϕ/2]]//Abs//Chop

HopV[3][0][pre ϕMaxval][[-1]]/.pre->1

Plot[Abs[comm[3,pre]-1],{pre,.5,2}]
Plot[Eigenvalues[HopV[3][0][pre ϕMaxval]][[-1]],{pre,.5,2}]

Get the Eigenvalues

ListPlot[{{3,.474},{4,.487},{5,.494},{6,.497},{7,.498}}]

Eigenvalues[HopV[1][0][ ϕMaxval]][[-1]]

Plot[Eigenvalues[HopV[1][0][ pre ϕMaxval]][[-1]],{pre,.5,2}]
Plot[Eigenvalues[HopV[2][0][ pre ϕMaxval]][[-1]],{pre,.5,2}]
Plot[Eigenvalues[HopV[3][0][ pre ϕMaxval]][[-1]],{pre,.5,2}]

Create the first few states of the free theory

vac=Eigensystem[Hop[3][0]][[2,-1]]//Chop
Print["Normalization = ",vac.vac//N//Chop]
Print["Energy = ",vac.Hop[3][0].vac//N//Chop]
ListPlot[Table[{ϕval[i],vac[[i+1]]},{i,0,nϕ-1}],Joined->True]
((Hop[3][0].vac)/vac)//N//Chop;
ListPlot[%,PlotRange->{.3,.7}]

first=adagop[3].vac;
Print["Normalization = ",Conjugate[first].first//N//Chop]
Print["Energy = ",Conjugate[first].Hop[3][0].first//N//Chop]
ListPlot[first//Chop,Joined->True]
((Hop[3][0].first)/first)//N//Chop;
ListPlot[%,PlotRange->{1,2}]

second=(Sqrt[(1/(Factorial[2]))])MatrixPower[adagop[3],2].vac;
Print["Normalization = ",Conjugate[second].second//N//Chop]
Print["Energy = ",Conjugate[second].Hop[3][0].second//N//Chop]
ListPlot[second//Chop,Joined->True]
((Hop[3][0].second)/second)//N//Chop;
ListPlot[%,PlotRange->{2,3}]

Perform time evolution in the free theory

vacEvolved[0]=evolve[3,30][0].vac;
Print["Normalization = ",Conjugate[vacEvolved[0]].vacEvolved[0]//N//Chop]
Print["Energy = ",Conjugate[vacEvolved[0]].Hop[3][0].vacEvolved[0]//N//Chop]
ListPlot[Abs[((vacEvolved[0])/vac)],Joined->True]
ListPlot[Arg[vacEvolved[0]],Joined->True]

firstEvolved[0]=evolve[3,30][0].first;
Print["Normalization = ",Conjugate[firstEvolved[0]].firstEvolved[0]//N//Chop]
Print["Energy = ",Conjugate[firstEvolved[0]].Hop[3][0].firstEvolved[0]//N//Chop]
ListPlot[Abs[((firstEvolved[0])/first)],Joined->True]
ListPlot[Arg[((firstEvolved[0])/first)],Joined->True]

secondEvolved[0]=evolve[3,30][0].second;
Print["Normalization = ",Conjugate[secondEvolved[0]].secondEvolved[0]//N//Chop]
Print["Energy = ",Conjugate[secondEvolved[0]].Hop[3][0].secondEvolved[0]//N//Chop]
ListPlot[Abs[((secondEvolved[0])/second)],Joined->True]
ListPlot[Arg[((secondEvolved[0])/second)],Joined->True]