Untitled

simpl = 0 # No simplified SYK Hamiltonian 

def make_Majorana(N):

    # Make 'N' Majorana fermions. Set of N Hermitian matrices psi_i, i=1,..N
    # obeying anti-commutation relations {psi_i,psi_j} = δ_{ij}
    
    X = np.array([[0, 1], [1, 0]])
    Y = np.array([[0, -1j], [1j, 0]])
    Z = np.array([[1, 0], [0, -1]])
    I = np.array([[1, 0], [0, 1]])

    psi = dict()

    if simpl == 0: 
        for i in range(1, N+1):

            if (i % 2) == 1:
                matlist = [Z] * int((i-1)/2)
                matlist.append(X)
                matlist = matlist + [I] * int((N/2 - (i+1)/2))
                psi[i] = 1/np.sqrt(2)*reduce(np.kron, matlist)
            else:
                matlist = [Z] * int((i - 2) / 2)
                matlist.append(Y)
                matlist = matlist + [I] * int((N/2 - i/2))
                psi[i] = 1/np.sqrt(2)*reduce(np.kron, matlist)

    if simpl == 1: 
        for i in range(1, N+1):

            if (i % 2) == 1:
                matlist = [I] * int((i-1)/2)
                matlist.append(X)
                matlist = matlist + [I] * int((N/2 - (i+1)/2))
                psi[i] = 1/np.sqrt(2)*reduce(np.kron, matlist)
            else:
                matlist = [I] * int((i - 2) / 2)
                matlist.append(Y)
                matlist = matlist + [I] * int((N/2 - i/2))
                psi[i] = 1/np.sqrt(2)*reduce(np.kron, matlist)


    for i in range(1, N+1):
        for j in range(1, N+1):

            if simpl == 0:
            # Not checking this for simplified H yet.  

                if i != j:
                    if np.allclose(psi[i] @ psi[j], -psi[j] @ psi[i]) == False:
                        print ("Does not satisfy algebra")

                if i == j:
                    if np.allclose(psi[i] @ psi[j] + psi[j] @ psi[i], np.eye(int(2**(N/2)))) == False:
                        print ("Does not satisfy algebra for i=j")

    return psi


def make_Hamiltonian(psi, N, instances, J_squared):
    
    # Creates multiple realisations of the SYK Hamiltonian
    # Variance of couplings is given by 'J_squared * 3!/N^3'.

    H = 0
    J = dict()
    sigma_sq = 6.*J_squared/(N**3)
    sigma = math.sqrt(sigma_sq)

    for i in range(1, N+1):
        for j in range(i+1, N+1):
            for k in range(j+1, N+1):
                for l in range(k+1, N+1):
                    J[i, j, k, l] = np.random.normal(loc=0, scale=sigma,size=instances)
                    M = psi[i] @ psi[j] @ psi[k] @ psi[l]
                    H = H + np.array([element * M for element in J[i, j, k, l]])

    return H