Possible bug in time dependent evolution

[M Cotrufo]

Hi, consider the following simple code, which describes an optical cavity pumped by an external classical field with a gaussian temporal shape. The parameter Input_TWidth  fixes the temporal width of the gaussian input. However, I found I very strange things while trying different values for this variable. For, e.g., Input_TWidth = 0.3/2.355 or Input_TWidth = 0.09/2.355, the response of the cavity is a nice quasi-gaussian smooth plot, as expected. For other values, like  Input_TWidth = 0.1/2.355 or  Input_TWidth = 0.2/2.355, the response of the cavity is completely unrealistic, with very low values and a piecewise graph.

Can you try to run this code with the aforementioned values and check the results?

Many thanks

# -*- coding: utf-8 -*-

"""

Created on Thu Dec 17 15:59:44 2015

@author: mcotrufo

"""

import matplotlib.pyplot as plt

import numpy as np

from qutip import *

from math import * 

c_const = 299792458

N = 20             # number of cavity fock states

wc = (c_const/(1300.03))*(2*pi) 

wl = wc;

delta1 = wl-wc

psi0 = tensor(basis(N,0)) 

QFactorCavity = 1*10**5 #QFactor cavity

KappaOUT = wc/QFactorCavity

# Parameters for pumping

KappaIN=1*KappaOUT 

FluxPhotons =1

Multiplicator = 1

E=sqrt(KappaIN)*sqrt(FluxPhotons)

Input_Tstart = 1 #nanoseconds

Input_TWidth = 0.3/2.355 #nanoseconds

#%%

def H1_coeffGauss(t, args):

    return np.exp(-((t-Input_Tstart)/(sqrt(2)*Input_TWidth)) ** 2)

sizeT=10000

T=2

tlist = np.linspace(0,T,sizeT)

a  = tensor(destroy(N))

H0 = -delta1*a.dag()*a

H1 = 1j*E*(a.dag()-a)

H = [H0,[H1,H1_coeffGauss]]

    

c_op_list = []

c_op_list.append(sqrt(KappaOUT) * a)

 

output = mesolve(H, psi0, tlist,c_op_list, [])

N_c = expect(a.dag() * a,output.states)

plt.rc('text', usetex=True)

plt.rc('font', family='arial',size=20)

fig, axes = plt.subplots(1, 1, figsize=(10,6))

plt.plot(tlist, N_c, label="Cavity")