I am trying to implement a physics paper related to Quantum Computing. The part that I am implementing right now is done numerically. In this, I want to evolve the Hamiltonian H1 which is written in the momentum space (i.e in the reciprocal space - not a physical space). The hamiltonian expression is in the image below - Equation S11 is the time evolution operator for H1 Equation above S13 is the effective hamiltonian Heff, which can be simplified using the equation S14 enter image description here U1(epsilon) is the periodized time evolution operator for 1D Topological Matter as given in S16 enter image description here S29 is the expression for winding number which has to be computed enter image description here
The ultimate goal is to implement this part and understand it. I have not got into the part for plotting this, but I can expect what values I should be obtaining for winding number for different values of alpha and beta - winding number will be zero only for alpha and beta, but I am getting it zero for any values that I am inputting into the expression in Mathematica. enter image description here
The errors that I am getting are - Unable to prove the integrating limits are real {0, T/2} enter image description here Another issue that I am facing is that, when I try plotting the winding number as a function of alpha and beta, its running for too long. The range of alpha and beta - -2 to +2.
I tried using Sympy for symbolic computation, but it is taking too long to compute. Here is the Code for Sympy-
import matplotlib.pyplot as plt
import numpy as np
from numpy import random
# import scipy.integrate as integrate
from sympy import symbols, Matrix
import sympy as sp
from sympy import series
tau_x = np.array([[0.0, 1.0], [1.0, 0.0]])
tau_y = np.array([[0.0, complex(0.0, -1.0)], [ complex(0.0, 1.0), 0.0]])
tau_z = np.array([[1.0, 0.0], [0.0, -1.0]])
# alpha, beta, T, kx, e= 0, 0, 2, 0, 0.1
J_0, J_S, kx, e, T, alpha, beta = symbols('J_0, J_S, kx, e, T, alpha, beta')
J_0 = 3*alpha/(4*T)
J_S = 3*beta/(2*T)
#### Define H and Heff - figure out how this is different from Floquet Hamiltonian H
M1 = complex(0, -1)*sp.Matrix([[0, (J_0*T/3)+(J_S*T*np.exp(kx*complex(0,1))/6)],[(J_0*T/3)+(J_S*T*np.exp(kx*complex(0,1))/6),0]])
M2 = (T/2)*complex(0, 1)*sp.Matrix([[2*np.pi/T, -(J_0*T/(3*np.log(e)))-(J_S*T*np.exp(kx*complex(0,1))/(6*np.log(e)))],[-(J_0*T/(3*np.log(e)))-(J_S*T*np.exp(kx*complex(0,1))/(6*np.log(e))),2*np.pi/T]])
Error in the Cell
from sympy import symbols
x, y = symbols('x,y')
a = x + 2*y
print(a*a)
(x + 2*y).subs({x:10, y:3})
Output of the above cell-
(x + 2*y)**2
16
ERROR! Session/line number was not unique in database. History logging moved to new session 223
Also, in numpy, I guess it does not exponentiate the matrix so that the final matrix evaluates close to the one obtained through Taylor Series Expansion, it simply takes exponent for each matrix element. Because exponentiating Identity matrix - np.exp([[1,0],[0,1]]) as [2.71, 1], [1, 2.71]
Is there any way which can be done using in-built function in python?
Here is the Mathematica Code that I wrote-
pauliX[n_Integer] :=
KroneckerProduct[IdentityMatrix[n/2], {{0, 1}, {1, 0}}]
pauliY[n_Integer] :=
KroneckerProduct[IdentityMatrix[n/2], {{0, -I}, {I, 0}}]
pauliZ[n_Integer] :=
KroneckerProduct[IdentityMatrix[n/2], {{1, 0}, {0, -1}}]
(*Print[Dimensions[pauliX[8]]]*)
T=2*Pi
H1[t_,kx_, beta_, alpha_] := Piecewise[{{(3*alpha/(4*T*Pi))*pauliX[2], 0 <= t <= T/3}, {(3*beta/(2*T*Pi))*(Cos[kx]*pauliX[2]-Sin[kx]*pauliY[2]), T/3 <= t <= 2*T/3}, {(3*alpha/(4*T*Pi))*pauliX[2], 2*T/3 <= t <= T}}]
U1[t_, kx_, beta_, alpha_] := Exp[-I*Integrate[H1[t1, kx, beta, alpha], {t1, 0, t}]]
Heff[kx_, epil_, beta_, alpha_] := If[epil==Pi,(I/T)*((2*Pi*I)+(Log[(U1[T/2, kx, beta, alpha])])/Log[(epil)]),(I/T)*(-(2*Pi*I)-(Log[(U1[T/2, kx, beta, alpha])])/Log[(epil)])]
U1epil[t_, kx_, epil_, beta_, alpha_] := U1[t, kx, beta, alpha] . (Exp[(I*t*Heff[kx, epil, beta, alpha])])
K = Tr[pauliZ[2]*Inverse[U1epil[T/2, kx, epil, beta, alpha]]*(D[U1epil[T/2,Pi, Pi, 1.8*Pi, 1.8*Pi]])]
Print[FullSimplify[K]]
kx =.;
Plot[Tr[pauliZ[2]*Inverse[U1epil[T/2, kx, Pi, 1.8*Pi, 1.8*Pi]]*(D[U1epil[T/2, kx, Pi, 1.8*Pi, 1.8*Pi], kx])],{kx,-3.14,3.14}]
v[kx_, epil_, beta_, alpha_] := (I/(4*Pi))*Integrate[Tr[pauliZ[2]*Inverse[U1epil[T/2, kx, epil, beta, alpha]]*(D[U1epil[T/2, kx, epil, beta, alpha], kx])], kx, -Pi, Pi]
Print[FullSimplify[U1epil[T/2, kx, Pi, 1.8*Pi,1.8*Pi]]]
Print[v[kx,Pi, 1.8*Pi, 1.8*Pi]]
Errors-
Integrate: Invalid integration variable or limit(s) in -\[Pi].
General: -3.13987 is not a valid variable.
General: -3.13987 is not a valid variable.
General: -3.13987 is not a valid variable.
General: "Further output of General::ivar will be suppressed during this calculation"
Any suggestions would be appreciated.