simulating PT symmetric coupled cavities

[mehrosadat ebrahimi]

Hi,

I am studying a system which has the properties of PT symmetric. In fact, two bosonic modes are coupled  where one of them is dissipated and the other one is amplified. For more details, you can take a look at

https://journals.aps.org/pra/abstract/10.1103/PhysRevA.92.053837.

I want to obtain the second order correlation function, the same as the reference. To this end, I have used the following code:

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from qutip import *
from IPython.display import Image

g  = 0.5        
kappab = -1    
kappaa =1    
X =0.1  
w=0.01
N=7
M=7

  def calulcate_avg_photons(N, D):

    b  = tensor(destroy(N), identity(M))
    a = tensor(identity(N), destroy(M))
    H = l *(b.dag() * b + a.dag() *a) + g * (b.dag() * a +  a.dag() * b) + w * (a + a.dag()) + X*(a.dag()*a)**2
   
    c_ops = [(kappaa)*b, (kappam)*a]
    
    rho_ss = steadystate(H,c_ops)
   
    n_cavity = expect(a.dag() * a, rho_ss)
 

    g2_cavity = expect(a.dag() * a.dag() * a * a, rho_ss) / (n_cavity ** 2)
  

    return n_cavity, g2_cavity

Gamma_vec = np.linspace(-5, 5, 200)
n_avg_vec = []
g2_vec = []
for l in Gamma_vec:
    n_avg, g2 = calulcate_avg_photons(N, l)
    n_avg_vec.append(n_avg)
    g2_vec .append(g2)

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

axes.plot(Gamma_vec, g2_vec, color="blue", alpha=0.6, label="numerical")

axes.set_xlabel(r'$\Delta/\kappa_{m}$', fontsize=18)
axes.set_ylabel(r'$g^{(2)}(0)$', fontsize=18)
axes.ticklabel_format(useOffset=False)

Although I have selected "kappab=-1",  meaning that mode "b" is amplified instead of dissipating, the result is exactly the same as the case where "kappaa = kappab = 1"!

Can anybody please help me to obtain the correct result the same as the reference?

Bests

Mehri

[Diego Alejandro Rasero Causil]

Best regard

I once had a similar problem. Simon and Neil suggested that I create my own Lindbladian. That way I managed to fix it.

I send part of the code that worked for me. I hope that, with the appropriate modifications, it will serve you

hbar = 1
N = 10 
gamma = 0.1 
g = 5.0*gamma  
kappa = gamma
omega = 0.1*gamma 

# cavity operators
a  = tensor(destroy(N), identity(2))
ad = tensor(create(N), identity(2))

# atomic operators
sz = tensor(identity(N), sigmaz())   #sigma_{z}
sx = tensor(identity(N), sigmax())   #sigma_{x}
sm = tensor(identity(N), destroy(2)) #sigma_{-}
sp = tensor(identity(N), create(2))

v_expected= [ad*a, sz]

tlist = np.linspace(0.0, 100000, 2000)
taulist = tlist

def calulcate_avg_photons(N, Delta):  

# Hamiltoniano

H = hbar*Delta*ad*a  +  hbar*Delta*sp*sm  +  hbar*g*(sm*ad + sp*a)  +  hbar*omega*(a + ad)
 # Liouvilliano

L  = - 1j*(spre(H) - spost(H)) + kappa*(-sprepost(a,ad) + 0.5*spre(ad*a) + 0.5*spost(ad*a)) + gamma*(sprepost(sm,sp) - 0.5*spre(sp*sm) - 0.5*spost(sp*sm))  

# Ground state and steady state

rho_ss = steadystate(L)  

# cavity photon number
n_cavity = expect(ad * a, rho_ss) 

# cavity second order coherence function
g2_cavity = expect(ad * ad * a * a, rho_ss) / (n_cavity * n_cavity)
return n_cavity, g2_cavity

Delta_max = 30
Delta_vec = np.linspace(-Delta_max, Delta_max, 10000)

n_avg_vec = []
g2_vec = []

for Delta in Delta_vec:

n_avg, g2 = calulcate_avg_photons(N, Delta)
n_avg_vec.append(n_avg)
g2_vec.append(g2)

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

axes.plot(Delta_vec / (gamma), np.abs(n_avg_vec), color="blue", alpha=1.0, label="numerical", linewidth=4.5)
axes.set_xlabel(r'$\Delta/\gamma$', fontsize=24)
axes.set_ylabel(r'Occupation probability $\langle n \rangle$', fontsize=24)
axes.set_xlim(-Delta_max, Delta_max);

[mehrosadat ebrahimi]

Thanks for your answer.

I used the Lindblad equation according to your suggestion. First of all, I checked the result for positive dissipation rate. As is expected it works. But for negative dissipation rate, it does not work correctly.

In fact, I have selected kappab = -1   and   kappaa =1. In this case, the second order correlation function is negative!

 I have also checked the commutation relation  (the standard bosonic commutator) in this case. It is near to one, but there are so fluctuations .

You can see the result of the commutation relation in the attached file, which is obtained by choosing the dimension N=M=10.

I have to mention that by considering positive dissipation rates, kappaa=kappab=1, the commutation relation is smooth and exactly equal to one even by opting smaller dimensions N=M=5.

Does anyone have any idea about what happened?

Bests

Mehri

[Fabrizio Minganti]

Hey guys,

I gave a quick look to the code and the paper. In my opinion, the problem is that in the article you cite the authors are using a non-Hermitian Hamiltonian, while you are using a Liouvillian. They say "Two approaches [non-Hermitian Hamiltonian and Liouvillian] can offer us the same solutions and here we demonstrate only the results obtained from the Schrödinger equation, as this method is much faster and computationally less demanding." I'm quite skeptical about this claim, especially for second-order correlation functions and in the presence of Kerr nonlinearity. Some examples of this are given in:

• Phys. Rev. A 100, 062131
  
• Phys. Rev. A 104, 012205 
  
• Phys. Rev. A 102, 033715
  

Nonetheless, some fundamental remarks are in order: 

1. The dissipation rates should always be positive numbers. Gain and dissipation in a Lindbladian picture are not described by choosing negative and positive numbers, but rather by their jump operators. Equivalently: the dissipation rates are the answer to the question "how often is particle gained or lost"? This cannot be a negative number.  Gain injects particle and is the operator b.dag(), loss subtract particle, and it is a. Thus, in your case, you should fix  c_ops = [(kappa)*b.dag(), (kappa)*a], with kappa=1. 
2. A non-Hermitian Hamiltonian description works for zero nonlinearity, where the equation of motion for the bosonic field moments are described, to the lowest order, by a non-Hermitian operator. Notice also that, in order to obtain negative and positive dissipation rates, you need to explicitly gauge out some decaying process. What I would try is to build a non-Hermitian Hamiltonian, obtain its ``steady state'', and see if you manage to obtain the same results in the article.

Hope this has been useful,

Fabrizio