superradiant_light_emission.rst

Superradiant Light Emission

We consider a system of N two-level systems (TLSs) with identical frequency ω₀, which can emit collectively at a rate γ_CE, and suffer from dephasing and local losses at rates γ_D and γ_E, respectively. The dynamics can be written as

$$\dot{\rho} =-i\lbrack \omega_{0}J_z,\rho \rbrack +\frac{\gamma_\text {CE}}{2}\mathcal{L}_{J_{-}}[\rho] +\sum_{n=1}^{N}\frac{\gamma_\text{D}}{2}\mathcal{L}_{J_{z,n}}[\rho] +\frac{\gamma_\text{E}}{2}\mathcal{L}_{J_{-,n}}[\rho].$$

When γ_E = γ_D = 0 this dynamics is the classical superradiant master equation. In this limit, a system initially prepared in the fully-excited state undergoes superradiant light emission whose peak intensity scales proportionally to N².

from qutip import *
from piqs import *
import matplotlib.pyplot as plt

N = 20
[jx, jy, jz] = jspin(N)
jp = jspin(N, "+")
jm = jp.dag()

# spin hamiltonian
w0 = 1
H = w0 * jz

# dissipation
gCE, gD, gE = 1, 1, 0

# set initial conditions for spins
system = Dicke(N=N, hamiltonian=h, dephasing=gD,
               collective_emission=gCE)

# build the Liouvillian matrix  
liouv = system.liouvillian()

Now that the system Liouvillian is defined, we can use QuTiP to solve the dynamics. We use as integration time a multiple of the superradiant delay time, t_D = log (N)/(Nγ_CE). We specify the operators for which the expectation values should be calculated to mesolve with the keyword argument e_ops. In this case, we are interested in J_x, J₊J₋, J_z².

nt = 1001
td0 = np.log(N)/(N*gCE)
tmax = 10 * td0

t = np.linspace(0, tmax, nt)

# initial state
excited_rho = excited(N)

# alternative states 
superradiant_rho = dicke(N, N/2, 0)
subradiant_rho = dicke(N, 0, 0)
css_symmetric = css(N)

a = 1/np.sqrt(2)
css_antisymmetric = css(N, a, -a)
ghz_rho = ghz(N)
rho0 = excited_rho
result = mesolve(liouv, rho0, t, [], e_ops = [jz, jp*jm, jz**2], 
                 options = Options(store_states=True))
rhot = result.states

We can then plot the results of the time evolution of the expectation values of the collective spin operators for different initial states.

jz_t = result.expect[0]
jpjm_t = result.expect[1]
jz2_t = result.expect[2]

jmax = (0.5 * N)
fig1 = plt.figure()
plt.plot(t/td0, jz_t/jmax)
plt.show()
plt.close()

References: