QutipQPU with Bosonic Hamiltonians

If we want to simulate the evolution of a number (or a Fock) state in a bosonic field, we can do so in the usual manner by encoding a Hamiltonian with a Schedule. Just in this case, to represent the creation and annihilation operators in an Observable, we will use the letters A and a for $a^\dagger$ and $a$, respectively. Let us take a look at an example with the Simple Harmonic Oscillator.

Simple harmonic oscillator (SHO)

The Hamiltonian of the SHO takes the form:

$$H = \hbar\omega \left(a^\dagger a + \frac12 \right)$$

where $a^\dagger$ and $a$ are the creation and annihilation operators, and $\omega$ is the oscillator frequency.

If we apply this Hamiltonian to a number state $\left|n\right>$ we would expect to get out exactly the same state with a value $\hbar\omega \left(n + \frac12 \right) $:

$$H\left|n\right> = \hbar\omega \left(n + \frac12 \right)\left|n\right> $$

Let's see how we could obtain these results with tools from the QLM.

Solving the SHO on the QLM

To evolve the initial state $\left|n\right>$ and examine the expected result, we can use the option SAMPLE in the to_job() method of the schedule. This will give us the probabilities of each possible state, assuming infinitely many measurements.

To get the value of the expected state $\left|n\right>$ we can use the OBS option by passing the Hamiltonian as the measured observable.

Specifying the Hilbert space

The QLM is a powerful computer, yet not a real infinite-level quantum system. So, to deal analytically with bosons we need to bound their excitation levels. This is done by the bosonic_levels argument of QutipQPU. It is a list of tuples for which each tuple consists of two numbers - the number of modes and the number of their respective excitations.

For a single SHO, we could say that we have one mode with at most 6 excitations. This can be described as bosonic_levels = [(1, 6)]. (If we were to work with a system containing, say, 5 modes where each has 10 excitations, and 8 modes each with 20 excitations, the respective representation would be bosonic_levels = [(5, 10), (8, 20)].)

Let us further specify that in this space our initial state will be $\left|n\right> = \left|2\right>$.

import numpy as np
from qat.core import Observable, Term, Schedule
from qat.qpus import QutipQPU

# Encode the SHO Hamiltonian in a Schedule
omega = 1.0
hamiltonian = Observable(1, pauli_terms=[Term(1, "Aa", [0, 0])], constant_coeff=1 / 2)
drive_sho = 1 * omega * hamiltonian  # taking h bar as 1
schedule = Schedule(drive=drive_sho,
                    tmax=1.0)

# Define the bosonic space structure
n_modes_SHO = 1
n_excitations_SHO = 6
bosonic_levels = [(n_modes_SHO, n_excitations_SHO)]

# Define the inital state of the boson
n_state = 2
psi_0_array = np.zeros(n_excitations_SHO) # the number state of the boson
psi_0_array[n_state] = 1

# Two different possible jobs depending on whether we want
# the final states or to measure an observable
job_sample = schedule.to_job(job_type="SAMPLE",  # "SAMPLE" can actually be omitted as it is the default option
                             psi_0=psi_0_array)
job_obs = schedule.to_job(job_type="OBS",
                          observable=hamiltonian,
                          psi_0=str(n_state))  # another way to specify the initial state - by giving it as a string

# Call the QPU by supplying with the Hilbert spaces
qpu = QutipQPU(bosonic_levels=bosonic_levels)

# Submit the jobs
res_state = qpu.submit(job_sample)
res_value = qpu.submit(job_obs)

We would expect the system to stay in the state $\left|2\right>$ and that the value of the Observable is 2.5. The amplitude of the state $\left|2\right>$ could be easily derived as the solution of this time-independedent Schrödinger equation is

$$\left|n(t)\right> = e^{-iHt / \hbar} \left|n(0)\right> = e^{-i \omega (n + 1/2) t} \left|n(0)\right>$$

Hence for $\left|n(0)\right> = \left|2\right>$ we should get $\left|n(t)\right> = e^{-i \omega (5/2) t} \left|2\right> \approx (-0.8011 - 0.5985i) \left|2\right>$ for $\omega = 1$ and $t = 1$.

from scipy.linalg import expm

# Constructing H in a matrix form, obeying the number of levels
H = np.zeros((n_excitations_SHO, n_excitations_SHO), dtype=complex)
for diag_index in range(n_excitations_SHO):
    H[diag_index, diag_index] = -1j * (diag_index + 0.5)
exp_H = expm(H)

print("Expected state and amplitude: |" + str(n_state) + ">", np.dot(exp_H, psi_0_array)[n_state])
print("Computed state and amplitude:", res_state.raw_data[0].state, res_state.raw_data[0].amplitude)

print("\nExpected observable value:", n_state + 0.5)
print("Computed observable value:", res_value.value)

Expected state and amplitude: |2> (-0.8011436155469337-0.5984721441039565j)
Computed state and amplitude: |2> (-0.8011461677843632-0.598468727540904j)

Expected observable value: 2.5
Computed observable value: 2.5

A qubit and a SHO

Let us now add a qubit to the system. With no coupling, we would expect the two systems to evolve independently. The total Hamiltonian could look like:

$$H = \alpha \, \sigma_x + \hbar\omega \left(a^\dagger a + \frac12 \right)$$

where $\alpha$ is a constant term and $\sigma_{x}$ is the usual Pauli matrix. With this operator, we would get a superposition for the qubit if it starts from an initial state $\left|0\right>$ and is evolved for time $\frac\pi4$ with $\alpha = 1$. If the boson starts from the state $\left|3\right>$ it should stay in that state.

A remark on bosonic_levels

We need to mention that when we have both qubits and bosons, the way we represent the space of the system is by still providing only the bosonic space structure in bosonic_levels. That is, if we have 3 qubits and (like before) 5 modes with 10 excitations, and 8 modes each having 20 excitations, the respective QutipQPU argument would still be bosonic_levels = [(5, 10), (8, 20)].

A remark on specifying psi_0

If one wants to specify the initial state of the total system (as a string), then declaring the states of the qubits (either 0 or 1) will be needed. But what should be the order of the qubits with respect to the modes?
The way to go is by declaring the state of the 2-level systems first (i.e. the qubits), then all modes with 3 excitations (if any), then all with 4 excitations (if any), etc. For example, if we want to start from 3 qubits in states $\left|1\right>$, $\left|1\right>$ and $\left|0\right>$, a 2 modes having 5 excitations each in states $\left|3\right>$ and $\left|2\right>$ and 1 mode with 9 excitations in state $\left|7\right>$, the string corresponding to the compound initial state would be 110327 and the respective space structure would be expressed as bosonic_levels = [(2, 5), (1, 9)].
The string representation wouldn't allow one to specify the state $\left|21\right>$ of a mode with 30 excitations, because of the possible ambiguity with the initial state of any smaller systems present. Even for such extreme cases though, one may still specify the initial state of the total system - via using an array representation of the desired psi_0 (which can be obtained by doing a Kronecker product between the states of the individual modes in the respective order).

With these remarks in mind, let us now evolve the above $H$:

# Hamiltonian terms
alpha = 1.0
omega = 1.0
evol_time_qubit_sho = 1 / 4 * np.pi

# Constructing the Hamiltonian of the system in a drive for the schedule
drive_qubit_sho = [(alpha, Observable(2, pauli_terms=[Term(1, "X", [0])])),
                   (1.0 * omega, Observable(2,
                                            pauli_terms=[Term(1, "Aa", [1, 1])],
                                            constant_coeff=1 / 2))]
# Create a Schedule and a Job
schedule = Schedule(drive=drive_qubit_sho, tmax=evol_time_qubit_sho)
psi_0 = '13' # start from |1> for the qubit and |3> for the boson
job = schedule.to_job(psi_0=psi_0)
                   
# Call the QPU
qpu = QutipQPU(bosonic_levels = [(1, 6)]) # one mode with 6 excitations

# Submit the job and get the result
result = qpu.submit(job)

sample_list = result.raw_data
for sample in sample_list:
    print(sample.state, sample.probability)

|0>|3> 0.5000004772383915
|1>|3> 0.4999995227616084

Adding a weak coupling to the qubit and the SHO

The SHO does not change the state of the boson. We can have the same behaviour for the qubit, namely with a $\sigma_z$ instead of a $\sigma_x$. Now, we can apply a coupling between the two quantum systems and see how the compound state would change. This could be achieved by the Hamiltonian:

$$H = \beta * \sigma_z + \hbar\omega \left(a^\dagger a + \frac12 \right) + g(\sigma^{+} a + \sigma^{-} a^\dagger)$$

where $\sigma^{+}$ and $\sigma^{-}$ are the spin-raising and lowering operators and $g$ is the coupling constant.

If we start again from the initial state $\left|1\right>\left|3\right>$ we would expect that for very weak coupling, i.e. very low $g$, the system will be in the state $\left|0\right>\left|2\right> = \sigma^{+} a \left|1\right>\left|3\right>$ only with a very low probability. However, this probability will get higher if $g$ is increased.

In the description of $H$ we will make use of the definitions of the raising and lowering operators as functions of the Pauli $\sigma_x$ and $\sigma_y$ matrices, namely $\sigma^{+} = \frac{\sigma_x + i\sigma_y}{2}$ and $\sigma^{-} = \frac{\sigma_x - i\sigma_y}{2}$. One should bear in mind that $\sigma^{+}\left|1\right> = \left|0\right>$ and $\sigma^{+}\left|0\right> = 0$. The inverse holds for $\sigma^{-}$.

# Hamiltonian terms
beta = 1.0
omega = 1.0
g_coupling = 0.01
evol_time_w_coupling = 1 / 4 * np.pi

# Constructing the Hamiltonian of the system in a drive for the schedule
drive_w_coupling = [(beta, Observable(2, pauli_terms=[Term(1, "Z", [0])])),
                   (1.0 * omega, Observable(2,
                                            pauli_terms=[Term(1, "Aa", [1, 1])],
                                            constant_coeff=1 / 2)),
                    # representing the spin-raising and lowering operators
                   (g_coupling, Observable(2, pauli_terms=[Term(1 / 2, "Xa", [0, 1]),
                                                           Term(1j / 2, "Ya", [0, 1]), 
                                                           Term(1 / 2, "XA", [0, 1]),
                                                           Term(-1j / 2, "YA", [0, 1])
                                                          ]))]
# Create a Schedule and a Job
schedule = Schedule(drive=drive_w_coupling, tmax=evol_time_w_coupling)
psi_0 = np.kron(np.array([0, 1]), # start from |1> for the qubit
                np.array([0, 0, 0, 1, 0, 0])) # and |3> for the boson
job = schedule.to_job(psi_0=psi_0)
                   
# Call the QPU
qpu = QutipQPU(bosonic_levels = [(1, 6)])

# Submit the job and get the result
result = qpu.submit(job)

/tmp/ipykernel_5827/714577833.py:22: UserWarning: psi_0 is converted into an array of floats. Prefer setting the ndarray dtype to float of complex to avoid this copy in the future.
  job = schedule.to_job(psi_0=psi_0)

sample_list = result.raw_data
for sample in sample_list:
    print(sample.state, sample.probability)

|0>|2> 0.00017572735922002792
|1>|3> 0.9998242726407799

The qubit in a superposition

We can now combine the two Hamiltonians above to see how even for a qubit in a superposition a similar result will be expected, i.e. some probability for the state $\left|0\right>\left|2\right>$. However, increasing the coupling this time will also cause "leakage" to the other possible states. The compound $H$ will look like:

$$H = \gamma * \sigma_x + \hbar\omega \left(a^\dagger a + \frac12 \right) + g(\sigma^{+} a + \sigma^{-} a^\dagger)$$

# Hamiltonian terms
gamma = 1.0
omega = 1.0
g_coupling = 0.1
evol_time_with_g = 1 / 4 * np.pi

# Constructing the Hamiltonian of the system in a drive for the schedule
drive_with_g = [(gamma, Observable(2, pauli_terms=[Term(1, "X", [0])])),
                (1.0 * omega, Observable(2,
                                         pauli_terms=[Term(1, "Aa", [1, 1])],
                                         constant_coeff=1 / 2)),
                (g_coupling, Observable(2, pauli_terms=[Term(1 / 2, "Xa", [0, 1]),
                                                        Term(1j / 2, "Ya", [0, 1]), 
                                                        Term(1 / 2, "XA", [0, 1]),
                                                        Term(-1j / 2, "YA", [0, 1])
                                                       ]))]
# Create a Schedule and a Job
schedule = Schedule(drive=drive_with_g,
                    tmax=evol_time_with_g)
job = schedule.to_job(psi_0='13')
                   
# Call the QPU
qpu = QutipQPU(bosonic_levels = [(1, 6)])

# Submit the job and get the result
result = qpu.submit(job)

# Get the Samples and sort them for nice printing
sample_list = result.raw_data
state_prob_list = [(sample.state, sample.probability) for sample in sample_list]
state_prob_list.sort(key=lambda x: x[1], reverse=True)
print("The states, sorted by probability:")
for state, probability in state_prob_list:
    # no constraint on the size of the probabilities would output them all
    if probability > 1e-5:
        print(state, probability)

The states, sorted by probability:
|0>|3> 0.4929770050456866
|1>|3> 0.4902990284801713
|0>|2> 0.011307946899560486
|1>|4> 0.0029635073478640663
|1>|2> 0.002227179916392618
|0>|4> 0.00022173355529079544

Further possibilities

One can build and simulate even more complicated Hamiltonians, for example:

• with time-dependent coefficients (like a qubit in a time-varying magnetic field)
• by adding a HardwareModel with
  • stochastic noise terms
  • Lindblad jump operators - see this example notebook for reference.