Negative second order coherence function !!!

[mehrosadate...@gmail.com]

Hello,

I want to plot the figure of second order coherence function for a special system. when I run my code, I see that some part of this figure is negative and this result is not physical !

My code is in this way :


from qutip import *

In [2]: import numpy as np

In [3]: import matplotlib.pyplot as plt

In [4]: N=5

In [5]: a=destroy(N)

In [6]: H=20.6*(a.dag()*a)+20.6*(a.dag()*a.dag()*a*a)+0.01*(a+a.dag())

In [7]: c_ops=[np.sqrt(1)*a]

In [8]: taus=np.linspace(0,100,1000)

In [9]: corr=coherence_function_g2(H,taus,c_ops,a)

In [10]: plt.plot(np.real(corr))
Out[10]: [<matplotlib.lines.Line2D at 0xb271d94c>]

In [11]: plt.legend(['two-photon blockade'])
Out[11]: <matplotlib.legend.Legend at 0xb27216cc>

In [12]: plt.xlabel(r'$\tau(ns)$')
Out[12]: <matplotlib.text.Text at 0xb2c48c0c>

In [13]: plt.ylabel(r'$g^{(2)}(\tau)$')
Out[13]: <matplotlib.text.Text at 0xb2c4a4ac>

In [14]: plt.show()



Is there any problem in my code, or the problem is from another things ?

Can you please help me ?


Regards
Mehrosadat

[Kevin Fischer]

Hi Mehrosadate,

This is the result of a numerical integration error. QuTiP doesn't automatically check for convergence and doesn't always select a sufficiently small integration step size. It uses a few heuristics but they don't necessarily work for some complicated Hamiltonians.

Importantly, the number of steps in taulist doesn't usually affect the integration step size. Instead, you should pass the Options() class to the coherence_function_g2() function. Decrease the max_step and increase nsteps and you will eventually get a smooth positive function for this simulation.

The other option (which I highly recommend here) is to change the solver from 'me' to 'es'. The eseries solver decomposes the time evolution of the master equation into exponentials, which works very well for steady state problems like this (and basically removes all numerical artifacts).

Kevin

[mehrosadate...@gmail.com]

Thanks you very much.

Now, I have a question about the second option that you suggested.
I did not exactly understand. what should I change in my code?

Bests
Mehrosadat

[Kevin Fischer]

Hi Mehrosadat,

You need to change:

coherence_function_g2(H,taus,c_ops,a)

->

coherence_function_g2(H,taus,c_ops,a, solver='es')

Kevin

[mehrosadate...@gmail.com]

Thank a lot.
I changed my code according to what you suggested, but I think that the first option was better.
Thanks because of this good group.

Regards
Mehrosadat