second order correlation function for a coupled atom cavity system as a function of detuning

[jakeb...@googlemail.com]

I need to calculate g2_0 for the following setup of a coupled atom cavity system with leaks to the environment over a range of different input frequencies w.

so far i have the following:

from qutip import*

%matplotlib inline

import matplotlib.pyplot as plt

import numpy as np

from cmath import *

# In units where h bar = 1

wc = 1.0  * 2 * np.pi                      # cavity frequency

wa = 1.0  * 2 * np.pi                      # atom frequency

g  = 1 * 2 * np.pi                            # coupling strength

kappa = 1                                      # cavity dissipation rate

gamma = 1                                    # atom dissipation rate

N = 20                                           # number of cavity fock states

a_in = 2                                         # Input power from laser

taus = np.linspace(0, 25.0, 200)    # calculation times

i=0  

for w in np.arange(-6 * np.pi, 6 * np.pi, np.pi):    # loop over different values of input frequency

    # initial state

    rho0 = tensor(basis(N,1), basis(2,0))              # start with atom in ground state

    # operators

    a  = tensor(destroy(N), qeye(2))

    sm = tensor(qeye(N), destroy(2))

    # Hamiltonian

    H_driven = (wc - w) * a.dag() * a + (wa - w) * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag()) + sqrt(kappa) * a_in * (a.dag() + a)

    # collapse operator

    c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

    # first calculate the occupation number as a function of time

    n = mesolve(H_driven, rho0, taus, c_ops, [a.dag() * a, sm.dag() * sm]).expect[0]

    # calculate the correlation function G2 and normalize with n(0)n(t) to obtain g2

    G2 = correlation_3op_1t(H_driven, rho0, taus, c_ops, a.dag(), a.dag() * a, a) 

    g2 = G2 / (n[0] * n)

    

    g2_0[i] = full(g2[0])

    

    i = i + 1

However I get the following error message:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-78-0b551972326d> in <module>()
     10     # Hamiltonian
     11 
---> 12     H_driven = (wc - w) * a.dag() * a + (wa - w) * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag()) + sqrt(kappa) * a_in * (a.dag() + a)
     13 
     14     # collapse operator

//anaconda/lib/python3.4/site-packages/qutip/qobj.py in __rmul__(self, other)
    466         if isinstance(other, (list, np.ndarray)):
    467             # if other is a list, do element-wise multiplication
--> 468             return np.array([item * self for item in other],
    469                             dtype=object)
    470 

TypeError: iteration over a 0-d array

Is there a way to do this without using a loops? Im not to sure about the end part where i try and turn the Qobj g2[0] into a numpy array so that i can index it and store g2[0] for each loop ready to plot against w. Im aware that you can't index Qobj so is there a better way to do this?

Any advise would be much appreciated as Im a beginner at computing in general as well as Qutip.

Thanks!