# setup the matplotlib graphics library and configure it to show # figures inline in the notebook %matplotlib inline import matplotlib.pyplot as plt import numpy as np # make qutip available in the rest of the notebook from qutip import * wc = 1.0 * 2 * pi # cavity frequency wa = 1.0 * 2 * pi # atom frequency g = 0.05 * 2 * pi # coupling strength kappa = 0.005 # cavity dissipation rate gamma = 0.05 # atom dissipation rate N = 15 # number of cavity fock states n_th_a = 0.0 # avg number of thermal bath excitation use_rwa = True tlist = np.linspace(0,25,101) # intial state psi0 = tensor(basis(N,0), basis(2,1)) # start with an excited atom # operators a = tensor(destroy(N), qeye(2)) sm = tensor(qeye(N), destroy(2)) # Hamiltonian if use_rwa: H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag()) else: H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() + a) * (sm + sm.dag()) c_ops = [] # cavity relaxation rate = kappa * (1 + n_th_a) if rate > 0.0: c_ops.append(sqrt(rate) * a) # cavity excitation, if temperature > 0 rate = kappa * n_th_a if rate > 0.0: c_ops.append(sqrt(rate) * a.dag()) # qubit relaxation rate = gamma if rate > 0.0: c_ops.append(sqrt(rate) * sm) output = mesolve(H, psi0, tlist, c_ops, [a.dag() * a, sm.dag() * sm]) n_c = output.expect[0] n_a = output.expect[1] fig, axes = plt.subplots(1, 1, figsize=(10,6)) axes.plot(tlist, n_c, label="Cavity") axes.plot(tlist, n_a, label="Atom excited state") axes.legend(loc=0) axes.set_xlabel('Time') axes.set_ylabel('Occupation probability') axes.set_title('Vacuum Rabi oscillations') output = mesolve(H, psi0, tlist, c_ops, []) output type(output.states) len(output.states) output.states[-1] # indexing the list with -1 results in the last element in the list # find the indices of the density matrices for the times we are interested in t_idx = where([tlist == t for t in [0.0, 5.0, 15.0, 25.0]])[1] tlist[t_idx] # get a list density matrices rho_list = array(output.states)[t_idx] # loop over the list of density matrices xvec = np.linspace(-3,3,200) fig, axes = plt.subplots(1,len(rho_list), sharex=True, figsize=(3*len(rho_list),3)) for idx, rho in enumerate(rho_list): # trace out the atom from the density matrix, to obtain # the reduced density matrix for the cavity rho_cavity = ptrace(rho, 0) # calculate its wigner function W = wigner(rho_cavity, xvec, xvec) # plot its wigner function axes[idx].contourf(xvec, xvec, W, 100, norm=mpl.colors.Normalize(-.25,.25), cmap=plt.get_cmap('RdBu')) axes[idx].set_title(r"$t = %.1f$" % tlist[t_idx][idx], fontsize=16) t_idx = where([tlist == t for t in [0.0, 5.0, 10, 15, 20, 25]])[1] rho_list = array(output.states)[t_idx] fig_grid = (2, len(rho_list)*2) fig = plt.figure(figsize=(2.5*len(rho_list),5)) for idx, rho in enumerate(rho_list): rho_cavity = ptrace(rho, 0) W = wigner(rho_cavity, xvec, xvec) ax = plt.subplot2grid(fig_grid, (0, 2*idx), colspan=2) ax.contourf(xvec, xvec, W, 100, norm=mpl.colors.Normalize(-.25,.25), cmap=plt.get_cmap('RdBu')) ax.set_title(r"$t = %.1f$" % tlist[t_idx][idx], fontsize=16) # plot the cavity occupation probability in the ground state ax = plt.subplot2grid(fig_grid, (1, 1), colspan=(fig_grid[1]-2)) ax.plot(tlist, n_c, label="Cavity") ax.plot(tlist, n_a, label="Atom excited state") ax.legend() ax.set_xlabel('Time') ax.set_ylabel('Occupation probability'); from qutip.ipynbtools import version_table version_table()