content: import numpy as np
  import matplotlib.pyplot as plt
  import qutip
  
  N = 15
  taus = np.linspace(0,10.0,200)
  a = qutip.destroy(N)
  H = 2 * np.pi * a.dag() * a
  
  # collapse operator
  G1 = 0.75
  n_th = 2.00  # bath temperature in terms of excitation number
  c_ops = [np.sqrt(G1 * (1 + n_th)) * a, np.sqrt(G1 * n_th) * a.dag()]
  
  # start with a coherent state
  rho0 = qutip.coherent_dm(N, 2.0)
  
  # first calculate the occupation number as a function of time
  n = qutip.mesolve(H, rho0, taus, c_ops, e_ops=[a.dag() * a]).expect[0]
  
  # calculate the correlation function G1 and normalize with n to obtain g1
  G1 = qutip.correlation_2op_1t(H, rho0, taus, c_ops, a.dag(), a)
  g1 = np.array(G1) / np.sqrt(n[0] * np.array(n))
  
  plt.plot(taus, np.real(g1), 'b', lw=2)
  plt.plot(taus, n, 'r', lw=2)
  plt.title('Decay of a coherent state to an incoherent (thermal) state')
  plt.xlabel(r'$\tau$')
  plt.legend([
      r'First-order coherence function $g^{(1)}(\tau)$',
      r'Occupation number $n(\tau)$',
  ])
  plt.show()
  