# Numerically finding the energy diagram of the hamiltonian

I'm looking at a collection of three two level systems (qubits) coupled to each other (with known bare state energies and couplings). The hamiltonian is given by $$\mathcal{H}=\sum_{i=1}^3{\omega_ia^\dagger_i a_i} + \sum_{i>j}{J_{ij}(a^\dagger_i a_j + a_i a^\dagger_j)}$$

Now, what I'm interested in is finding the energy levels of this system in terms of the states $|000\rangle$, $|001\rangle$, ..., $|111\rangle$, and I'll be trying to do so using Mathematica. I'm already not sure how to start.

In my mind, what I should do is write the qubit states as kronecker products so that the state vector of the system is given by a vector of length 8. Then I should subsequently express the hamiltonian as an eight by eight matrix, but I'm not sure how to do so. I can probably introduce Pauli matrices somehow, but how does that work in such a three qubit case? I'm sure this is a textbook example of numerically analysing a hamiltonian, but I haven't been able to find the recipe.

• When you say you've got three qubits, but then represent them by bosonic creation and annihilation operators... yeah, that doesn't make that much sense. – Emilio Pisanty Mar 8 '16 at 11:34
• Hm okay, well lets see where I'm being imprecise. They are essentially three coupled harmonic oscillators (very small anharmonicity) of which we only use the first two levels. Omega is the transition frequency of the ground state to the first excited state. – user129412 Mar 8 '16 at 11:41

You have 3 steps:
1. Build the hermitian coupling matrix $C=C^{\dagger}$-
$$C=\begin{pmatrix} \omega_1, J_{1,2},J_{1,3} \\ J_{2,1},\omega_2,J_{2,3} \\ J_{3,1},J_{3,1}, \omega_3 \end{pmatrix}$$

1. Build 3 anihilation matrices $a_n$ (n=1,2,3, acts on n'th mode) (each is $8\times8$ since you have 8 states). Choose an ordered basis, for example -
$$\begin{pmatrix} |000\rangle \\ |001\rangle \\ |010\rangle \\ |011\rangle \\ |100\rangle \\ |101\rangle \\ |110\rangle \\ |111\rangle\end{pmatrix}$$ and build each element of each matrix by mapping of the $a_n$ operator - for example $a_{2}|a,b,c\rangle=\sqrt{b}|a,b-1,c\rangle$

2. $$H=\sum_{i=1} ^3\sum_{j=1} ^3C_{i,j}a^{\dagger}_{i}a_{j}$$

Now you got your Hamiltonian and can find it's eigenvalues and eigenstates (both symbolically and numerically). The procedure is easily extendible to system of any size andnumber of qubits.

• This makes sense. I guess then the main point is how to do step two. Do I just take the two qubit annihilation and creation and Kronecker product them with identity matrices? – user129412 Mar 8 '16 at 11:49
• It can be implemented in many ways, choose the easiest for yourself. The trivial implementation is setting matrix of zeros and building the non zero elements (maximum one per column) by "for" loop, checking the mapping of the operator on each basis vector. – Alexander Mar 8 '16 at 11:55