How to realize Su-Schrieffer-Heeger model in Qiskit

Question

This is a very specific question, which I try to implement a simple dimerized tight-binding Hamiltonian on qiskit. The model is one dimensional, and defined as

$$ H = \sum_{\langle i,j\rangle} t_{ij} c^\dagger_i c_j $$

where $\langle i,j\rangle$ denotes the nearest neighbor coupling with strength $t_{ij}$. The important aspect of the model is that the coupling is dimerized in the following sense. If we label the sites as 1,2,3,4...., then the nearest neighbor coupling $t_{ij}$ reads, for example,

$$ t_{1,2} = 1, t_{2,3} = 2 , t_{3,4} = 1 , t_{4,5} = 2, etc $$

which is alternating.

I try to pretend this is a molecular Hamiltonian and realize it in the same way with "FermionicOperator". However, I am not sure how to properly index the "one-body" integral here, which is the nearest neighbor coupling $t_{ij}$ here.

Answer 1

I don't think it really matter how you index your $t_{ij}$. As you mentioned you can use FermionicOperator all you need to do is define the one body integral. This of course can be done in many ways. Here is a convenient way.

def ssh_ham(gamma, lamda, n): 
    sigmax = np.array([[0,1],[1,0]], dtype=np.complex_)
    sigmay = np.array([[0,-1j],[1j,0]], dtype=np.complex_)

    op_eye_x = np.eye(n)
    op_cos_x = 1/2*(np.eye(n, k=1) + np.eye(n,k=-1))
    op_sin_x = 1j/2*(np.eye(n, k=1) - np.eye(n,k=-1))

    h = np.kron(gamma*op_eye_x + lamda*op_cos_x, sigmax) + np.kron(lamda*op_sin_x, sigmay)

    return h 

This function would return the one-body integral for you. In this function, the $\sigma$ degree of freedom represents the two atoms in the unit cell, and $n$ represents how many unit cells you have in the SSH chain, and for the dimmerized limit you asked for you can take $\gamma = 1,\ \lambda = 2 $. You can use the output of this function fer_op = FermionicOperator(h1 = ssh_ham(1,2,10))

That being said, I'm not sure what the benefit would be. Trying to solve this Hamiltonian using VQE for example is an over-kill since this Hamiltonian can be readily diagonalized from its one-body Hamiltonian. Putting this on a quantum computer would give you $2n$ qubits, this is a Hilbert space of $2^n$ dimensions, whereas you can diagonalize the the Hamiltonian by diagonalizing the $2n \times 2n$ matrix $t_{ij}$.

Answer 2

This is the code I wrote, and I do feel that the index of the single particle integral is important. Below, h1 is the single particle integral in the basis ABABAB, where A,B are the sublattices (similar to spin), and I use "unitary_matrix" to change it AAA...BBB... which is the basis for h1_qubit. The construction of "unitary_matrix" can be found in the qiskit function "FermionicOperator".

Now, one important checking is that the ground state energy of two Hamiltonians, one is the fermionic single particle Hamiltonian (which is h1), and one is the Hamiltonian in the qubit basis (which is qubitOp), are identical. For the latter, it is a many-body Hamiltonian, and the ground state is the many-body ground state. If we naively compare this to the "lowest energy" of h1, then we will find that they are different, for h1 is the single-particle Hamiltonian. To properly obtain the many-body ground state from h1, we need to fill all the states below the Fermi-level of the spectrum of h1. Once we have done that, we found that the two Hamiltonians give the same many-body ground state energy.

import numpy as np
from scipy.linalg import eig
from qiskit.aqua.algorithms import VQE, NumPyEigensolver
from qiskit.chemistry import FermionicOperator
from qiskit.aqua.operators import Z2Symmetries


t1 = np.random.rand() # intra
t2 = np.random.rand() # inter
t0 = np.random.rand() # onsite energy 


N = 4 # number of unit cell

BC = 1 # BC=1 is PBC, BC=0 is OBC

h1 = np.zeros( (2*N , 2*N ) )


for i in range(0,(2*N-1)):
    if i % 2==1:
        h1[ i , i ] = t0
        h1[ i , i+1 ] = t2       
        h1[ i+1 , i ] = t2        
    else:
        h1[ i , i ] = -t0        
        h1[ i , i+1 ] = t1
        h1[ i+1 , i ] = t1        

h1[2*N-1,2*N-1] = t0

if BC==1:
    h1[ 0 , 2*N-1 ] = t2
    h1[ 2*N-1 , 0 ] = t2    
        
        
print( 'The fermionic Hamiltonian is\n' , h1 ) 

eigh1 = eig(h1)[0]

idx = np.argsort(eigh1)
eigh1 = eigh1[idx]

print( 'The *single* particle spectrum is \n' , eigh1 )

unitary_matrix = np.zeros((h1.shape), h1.dtype)
n = unitary_matrix.shape[0]
j = np.arange(n // 2)
unitary_matrix[j, 2 * j] = 1.0
unitary_matrix[j + n // 2, 2 * j + 1] = 1.0

h1_qubit = unitary_matrix.dot(h1.dot(unitary_matrix.T.conj()))


ferOp = FermionicOperator(h1=h1_qubit )

qubitOp = ferOp.mapping(map_type='parity' , threshold=0.00000001)

numEigenvalue = 1
result = NumPyEigensolver( qubitOp ,  k = numEigenvalue ).run()
eigvals = result['eigenvalues']
eigvecs = result['eigenstates']
print( 'The *many-body* grond state energy is \n ' , eigvals[0] )

print( 'We check the *many-body* ground state energy is the same as the sum of the energy below the fermi level of the *single* particle spectrum')

print( eigvals[0] - sum( eigh1[0:int(len(eigh1)/2)]))