Welcome to piqs’s documentation!

Introduction

Permutational Invariant Quantum Solver (PIQS)

PIQS is an open-source Python solver to study the exact Lindbladian dynamics of open quantum systems consisting of identical qubits.

In the case where local processes are included in the model of a system’s dynamics, numerical simulation requires dealing with density matrices of size \(2^N\). This becomes infeasible for a large number of qubits. We can simplify the calculations by exploiting the permutational invariance of indistinguishable quantum particles which allows the user to study hundreds of qubits.

Integrated with QuTiP

A major feature of PIQS is that it allows to build the Liouvillian of the system in an optimal way. It uses Cython to optimize performance and by taking full advangtage of the sparsity of the matrix it can deal with large systems. Since it is compatible with the quantum object class of [QuTiP] one can take full advantage of existing features of this excellent open-source library.

A wide range of applications

• The time evolution of the total density matrix of quantum optics and cavity QED systems for permutationally symmetric initial states (such as the GHZ state, Dicke states, coherent spin states).
• Quantum phase transitions (QPT) of driven-dissipative out-of-equilibrium quantum systems.
• Correlation functions of collective systems in quantum optics experiments, such as the spectral density and second-order correlation functions.
• Various quantum optics phenomena such as steady-state superradiance, superradiant light emission, superradiant phase transition, spin squeezing, boundary time crystals, resonance fluorescence.

Installation

In the terminal enter the following commands (you just need git and python installed). If you do not have git installed, just download the folder from Github and run the setup.py file with python. Please install cython, numpy, scipy and qutip as piqs depends on these packages.

We will soon publish the code in the Python Packaging Index (pip) and also make a conda package for easy installation on Windows. If you have any problems installing the tool, please open an issue or write to us.

git clone https://github.com/nathanshammah/piqs.git
cd piqs
python setup.py install

User Guide

The Permutational Invariant Quantum Solver (PIQS) is an open-source Python solver to study the exact Lindbladian dynamics of open quantum systems consisting of identical qubits. It is integrated in QuTiP and can be imported as as a model.

Using this library, the Liouvillian of an ensemble of \(N\) qubits, or two-level systems (TLSs), \(\mathcal{D}_{TLS}(\rho)\), can be built using only polynomial – instead of exponential – resources. This has many applications for the study of realistic quantum optics models of many TLSs and in general as a tool in cavity QED [1].

Consider a system evolving according to the equation

\[ \begin{align}\begin{aligned}\dot{\rho} = \mathcal{D}_\text{TLS}(\rho)=-\frac{i}{\hbar}\lbrack H,\rho \rbrack +\frac{\gamma_\text{CE}}{2}\mathcal{L}_{J_{-}}[\rho] +\frac{\gamma_\text{CD}}{2}\mathcal{L}_{J_{z}}[\rho] +\frac{\gamma_\text{CP}}{2}\mathcal{L}_{J_{+}}[\rho]\\+\sum_{n=1}^{N}\left( \frac{\gamma_\text{E}}{2}\mathcal{L}_{J_{-,n}}[\rho] +\frac{\gamma_\text{D}}{2}\mathcal{L}_{J_{z,n}}[\rho] +\frac{\gamma_\text{P}}{2}\mathcal{L}_{J_{+,n}}[\rho]\right)\end{aligned}\end{align} \]

where \(J_{\alpha,n}=\frac{1}{2}\sigma_{\alpha,n}\) are SU(2) Pauli spin operators, with \({\alpha=x,y,z}\) and \(J_{\pm,n}=\sigma_{\pm,n}\). The collective spin operators are \(J_{\alpha} = \sum_{n}J_{\alpha,n}\) . The Lindblad super-operators are \(\mathcal{L}_{A} = 2A\rho A^\dagger - A^\dagger A \rho - \rho A^\dagger A\).

The inclusion of local processes in the dynamics lead to using a Liouvillian space of dimension \(4^N\). By exploiting the permutational invariance of identical particles [2-8], the Liouvillian \(\mathcal{D}_\text{TLS}(\rho)\) can be built as a block-diagonal matrix in the basis of Dicke states \(|j, m \rangle\).

The system under study is defined by creating an object of the Dicke class, e.g. simply named system, whose first attribute is

• system.N, the number of TLSs of the system \(N\).

The rates for collective and local processes are simply defined as

• collective_emission defines \(\gamma_\text{CE}\), collective (superradiant) emission
• collective_dephasing defines \(\gamma_\text{CD}\), collective dephasing
• collective_pumping defines \(\gamma_\text{CP}\), collective pumping.
• emission defines \(\gamma_\text{E}\), incoherent emission (losses)
• dephasing defines \(\gamma_\text{D}\), local dephasing
• pumping defines \(\gamma_\text{P}\), incoherent pumping.

Then the system.lindbladian() creates the total TLS Linbladian superoperator matrix. Similarly, system.hamiltonian defines the TLS hamiltonian of the system \(H_\text{TLS}\).

The system’s Liouvillian can be built using system.liouvillian(). The properties of a Piqs object can be visualized by simply calling system. We give two basic examples on the use of PIQS. In the first example the incoherent emission of N driven TLSs is considered.

from piqs import Dicke
from qutip import steadystate
N = 10
system = Dicke(N, emission = 1, pumping = 2)
L = system.liouvillian()
steady = steadystate(L)

Superradiant Light Emission

We consider a system of \(N\) two-level systems (TLSs) with identical frequency \(\omega_{0}\), which can emit collectively at a rate \(\gamma_\text{CE}\), and suffer from dephasing and local losses at rates \(\gamma_\text{D}\) and \(\gamma_\text{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 \(\gamma_\text{E}=\gamma_\text{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^2\).

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_\text{D}=\log(N)/(N \gamma_\text{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^2\).

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:

[1]	Dicke, R. H. “Coherence in Spontaneous Radiation Processes”. Phys. Rev. 93, 91 (1954) doi/10.1103/PhysRev.93.99.

[2]	Bonifacio, R. and Schwendimann, P. and Haake, Fritz. “Quantum Statistical Theory of Superradiance. I.” Phys. Rev. A 4, 302 (1971) doi:10.1103/PhysRevA.4.302

Superradiance: Qubits in a cavity

We consider a system of \(N\) two-level systems (TLSs) coupled to a cavity mode. This is known as the Dicke model

\[H = \omega_{0}J_z + \omega_{c}a^\dagger a + g\left(a^\dagger + a\right)\left(J_{+} + J_{-}\right)\]

where each TLS has identical frequency \(\omega_0\). The light matter coupling can be in the ultrastrong coupling (USC) regime, \(g/ \omega_0 >0.1\).

If we study this model as an open quantum system, the cavity can leak photons and the TLSs are subject to local processes. For example the system can be incoherently pumped at a rate \(\gamma_\text{P}\), the TLSs are subject to dephaisng at a rate \(\gamma_\text{D}\), and local incoherent emission occurs at a rate \(\gamma_\text{E}\). The dynamics of the coupled light-matter system is governed by

\[ \begin{align}\begin{aligned}\dot{\rho} = -i\lbrack \omega_{0}J_z + \omega_{a}a^\dagger a + g\left(a^\dagger + a\right)\left(J_{+} + J_{-}\right),\rho \rbrack\\+\frac{\kappa}{2}\mathcal{L}_{a}[\rho] +\sum_{n=1}^{N}\left(\frac{\gamma_\text{P}}{2}\mathcal{L}_{J_{+,n}}[\rho] +\frac{\gamma_\text{E}}{2}\mathcal{L}_{J_{+,n}}[\rho] +\frac{\gamma_\text{D}}{2}\mathcal{L}_{J_{+,n}}[\rho]\right) \ \ \ \ \ \ (1)\end{aligned}\end{align} \]

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

#TLS parameters
N = 6
nds = num_dicke_states(N)
[jx, jy, jz] = jspin(N)
jp, jm = jspin(N, "+"), jspin(N, "-")
w0 = 1
gE, gD = 0.1, 0.01

# Hamiltonian
h = w0 * jz

#photonic parameters
nphot = 20
wc = 1
kappa = 1
ratio_g = 2
g = ratio_g/np.sqrt(N)
a = destroy(nphot)

After defining all the parameters, we can build a Liouvillian for the TLS ensemble and the photonic cavity. In order to study this system using QuTiP and \(PIQS\), we will first build the TLS Liouvillian, then we will build the photonic Liouvillian and finally we will build the light-matter interaction. The total dynamics of the system is thus defined in a Liouvillian space that has both TLS and photonic degrees of freedom.

#TLS liouvillian
ensemble = dicke(N = N, hamiltonian=h, emission=gE, dephasing=gD)
liouv = ensemble.liouvillian()

#photonic liouvilian
h_phot = wc * a.dag() * a
c_ops_phot = [np.sqrt(kappa) * a]
liouv_phot = liouvillian(h_phot, c_ops_phot)

We can then make a light-matter superoperator to address the total system of N spins and the photonic cavity by the super_tensor function in QuTiP. Similarly, the Liouvillian for the interaction Hamiltonian can be constructed with the spre and spost functions representing pre and post multiplication super-operators to finally construct the total Liouvillian of the combined light-matter system.

A similar treatment is possible for any operator and we construct the total \(J_z, J_+ J_-\) superoperators.

When only the dissipation of the cavity is present, beyond a critical value of the coupling \(g\), the steady state of the system becomes superradiant. This is visible by looking at the Wigner function of the photonic part of the density matrix, which displays two displaced lobes in the \(x\) and \(p\) plane.

rho_steady_state = steadystate(liouv_tot)
jz_steady_state = expect(jz_tot, rho_steady_state)
jpjm_steady_state = expect(jpjm_tot, rho_steady_state)

nphot_steady_state = expect(nphot_tot, rho_steady_state)
psi = rho_steady_state.ptrace(0)
xvec = np.linspace(-6, 6, 100)
W = wigner(psi, xvec, xvec)

wmap = wigner_cmap(W)  # Generate Wigner colormap
nrm = mpl.colors.Normalize(0, W.max())
plt.contourf(xvec, xvec, W, 100, cmap=wmap, norm=nrm)
plt.show()

As it has been shown in Ref. [1], the presence of dephasing suppresses the superradiant phase transition, while the presence of local emission restores it [2].

References:

[1]	Kirton, Peter, and Jonathan Keeling. “Suppressing and restoring the dicke superradiance transition by dephasing and decay.” Physical review letters 118.12 (2017): 123602.

Spin squeezing

PIQS can be used to study spin squeezing and the effect of collective and local processes on a spin squeezing Hamiltonian such as:

\[H = -i\Lambda\left(J_{+}^2-J_{-}^2\right)\]

which evolves under the dynamics given by:

\[\dot{\rho} = -\frac{i}{\hbar} \lbrack H,\rho \rbrack +\frac{\gamma_\text{CE}}{2}\mathcal{L}_{J_{-}} + \frac{\gamma_\text{E}}{2}\sum_{n=1}^{N}\mathcal{L}_{J_{-,n}}[\rho].\]

In [1] it has been shown that the collective emmission (\(\gamma_\text{CE}\)) affects the spin squeezing in a system in a different way than the homogeneous local emission (\(\gamma_\text{E}\)). In PIQS, we can study these effects easily by adding these rates to an ensemble constructed as a Dicke object.

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

# general parameters
N = 20
nds = num_dicke_states(N)
[jx, jy, jz] = jspin(N)
jp, jm = jspin(N, "+"), jspin(N, "-")
jpjm = jp*jm

lam = 1
# spin hamiltonian
h = -1j*lam*(jp**2-jm**2)

gamma = 0.2

# Ensemble with collective emission only
ensemble_ce = Dicke(N=N, hamiltonian=h, collective_emission=gamma)

# Ensemble with local emission only
ensemble_le = Dicke(N=N, hamiltonian=h, emission=gamma)

# Build the Liouvillians for both ensembles
liouv_collective = ensemble_ce.liouvillian()
liouv_local = ensemble_le.liouvillian()

Once we have defined our ensembles and constructed their Liouvillians, we can plot the time evolution of the spin squeezing parameter given by \(\xi^2= \frac{N \langle\Delta J_y^2\rangle}{\langle J_z\rangle^2}\) starting from any initial state.

# set initial state for spins (Dicke basis)
rho0 = dicke(N, 10, 10)
t = np.linspace(0, 2.5, 1000)

result_collective = mesolve(liouv_collective, excited, t, [],
                 e_ops = [jz, jy, jy**2,jz**2, jx])
result_local = mesolve(liouv_local, excited, t, [],
                 e_ops = [jz, jy, jy**2,jz**2, jx])

# Get the expectation values
jzt_c, jyt_c, jy2t_c, jz2t_c, jxt_c = result_collective.expect
jzt_l, jyt_l, jy2t_l, jz2t_l, jxt_l = result_local.expect

del_jy_c = jy2t_c - jyt_c**2
del_jy_l = jy2t_l - jyt_l**2

xi2_c = N * del_jy_c/(jzt_c**2 + jxt_c**2)
xi2_l = N * del_jy_l/(jzt_l**2 + jxt_l**2)

# Generate the plots
plt.plot(t*N*lam, xi2_c, 'k-', label="Collective emission")
plt.plot(t*N*lam, xi2_l, 'r--', label="Local_emission")
plt.plot(t*N*lam, 1+0*t, '--k')
plt.ylabel(r'$\xi^2$')
plt.xlabel(r'$ N \Lambda t$')
plt.legend()

plt.xlim([0, 2])
plt.ylim([0, 2])
plt.show()

References:

[1]	Chase and J. Geremia, Collective processes of an ensemble of spin-1 particles, Phys. Rev. A 78, 052101 2 (2008).

PIQS functions

Operators	Command	Description
Collective spin Jx	jspin(N, "x")	The collective spin operator Jx. Requires N number of TLS
Collective spin J+	jspin(N, "+")	The collective spin operator J+.
Collective spin J-	jspin(N, "-")	The collective spin operator Jz.
Collective spin Jx in uncoupled basis	jspin(N, "z", basis='uncoupled')	The collective spin operator Jz in the uncoupled basis
Dicke state |j, m>	dicke(N, j, m)	A Dicke state given by |j, m>
Excited state in uncoupled basis	excited(N, basis="uncoupled")	The excited state in the uncoupled basis
GHZ state in the Dicke basis	ghz(N)	The GHZ state in the Dicke (default) basis for N number of TLS
Collapse operators of the ensemble	Dicke.c_ops()	The collapse operators for the ensemble can be called by the c_ops method of the dicke class.

Documentation

Dicke module

Permutational Invariant Quantum Solver (PIQS)

This module calculates the Liouvillian for the dynamics of ensembles of identical two-level systems (TLS) in the presence of local and collective processes by exploiting permutational symmetry and using the Dicke basis.

class piqs.dicke.Dicke(N, hamiltonian=None, emission=0.0, dephasing=0.0, pumping=0.0, collective_emission=0.0, collective_dephasing=0.0, collective_pumping=0.0)

The Dicke class which builds the Lindbladian and Liouvillian matrix.

Example

>>> from piqs import Dicke, jspin
>>> N = 2
>>> jx, jy, jz = jspin(N)
>>> jp = jspin(N, "+")
>>> jm = jspin(N, "-")
>>> ensemble = Dicke(N, emission=1.)
>>> L = ensemble.liouvillian()

Parameters:

• N (int) -- The number of two-level systems.
• hamiltonian --

  A Hamiltonian in the Dicke basis.

  The matrix dimensions are (nds, nds), with nds being the number of Dicke states. The Hamiltonian can be built with the operators given by the jspin functions.

• emission (float) -- Incoherent emission coefficient (also nonradiative emission). default: 0.0
• dephasing (float) -- Local dephasing coefficient. default: 0.0
• pumping (float) -- Incoherent pumping coefficient. default: 0.0
• collective_emission (float) -- Collective (superradiant) emmission coefficient. default: 0.0
• collective_pumping (float) -- Collective pumping coefficient. default: 0.0
• collective_dephasing (float) -- Collective dephasing coefficient. default: 0.0

N

int -- The number of two-level systems.

hamiltonian

:class: qutip.Qobj -- A Hamiltonian in the Dicke basis.

The matrix dimensions are (nds, nds), with nds being the number of Dicke states. The Hamiltonian can be built with the operators given by the jspin function in the "dicke" basis.

emission

float -- Incoherent emission coefficient (also nonradiative emission). default: 0.0

dephasing

float -- Local dephasing coefficient. default: 0.0

pumping

float -- Incoherent pumping coefficient. default: 0.0

collective_emission

float -- Collective (superradiant) emmission coefficient. default: 0.0

collective_dephasing

float -- Collective dephasing coefficient. default: 0.0

collective_pumping

float -- Collective pumping coefficient. default: 0.0

nds

int -- The number of Dicke states.

dshape

tuple -- The shape of the Hilbert space in the Dicke or uncoupled basis. default: (nds, nds).

__repr__()

Print the current parameters of the system.

c_ops()

Build collapse operators in the full Hilbert space 2^N.

Returns:c_ops_list -- The list with the collapse operators in the 2^N Hilbert space. Return type:list

coefficient_matrix()

Build coefficient matrix for ODE for a diagonal problem.

Returns:M -- The matrix M of the coefficients for the ODE dp/dt = M p. p is the vector of the diagonal matrix elements of the density matrix rho in the Dicke basis. Return type:ndarray

lindbladian()

Build the Lindbladian superoperator of the dissipative dynamics.

Returns:lindbladian -- The Lindbladian matrix as a qutip.Qobj. Return type:

class:qutip.Qobj

liouvillian()

Build the total Liouvillian using the Dicke basis.

Returns:liouv -- The Liouvillian matrix for the system. Return type:

class:qutip.Qobj

pisolve(initial_state, tlist, options=None)

Solve for diagonal Hamiltonians and initial states faster.

Parameters:

• initial_state -- An initial state specified as a density matrix of qutip.Qbj type
• tlist (ndarray) -- A 1D numpy array of list of timesteps to integrate
• options -- The options for the solver.

Returns:

result -- A dictionary of the type qutip.solver.Result which holds the results of the evolution.

Return type:

list

prune_eigenstates(liouvillian)

Remove spurious eigenvalues and eigenvectors of the Liouvillian.

Spurious means that the given eigenvector has elements outside of the block-diagonal matrix.

Parameters:liouvillian_eigenstates (list) -- A list with the eigenvalues and eigenvectors of the Liouvillian including spurious ones. Returns:correct_eigenstates -- The list with the correct eigenvalues and eigenvectors of the Liouvillian. Return type:list

class piqs.dicke.Pim(N, emission=0.0, dephasing=0, pumping=0, collective_emission=0, collective_pumping=0, collective_dephasing=0)

The Permutation Invariant Matrix class.

Initialize the class with the parameters for generating a Permutation Invariant matrix which evolves a given diagonal initial state p as:

dp/dt = Mp

Parameters:

• N (int) -- The number of two-level systems.
• emission (float) -- Incoherent emission coefficient (also nonradiative emission). default: 0.0
• dephasing (float) -- Local dephasing coefficient. default: 0.0
• pumping (float) -- Incoherent pumping coefficient. default: 0.0
• collective_emission (float) -- Collective (superradiant) emmission coefficient. default: 0.0
• collective_pumping (float) -- Collective pumping coefficient. default: 0.0
• collective_dephasing (float) -- Collective dephasing coefficient. default: 0.0

N

int -- The number of two-level systems.

emission

float -- Incoherent emission coefficient (also nonradiative emission). default: 0.0

dephasing

float -- Local dephasing coefficient. default: 0.0

pumping

float -- Incoherent pumping coefficient. default: 0.0

collective_emission

float -- Collective (superradiant) emmission coefficient. default: 0.0

collective_dephasing

float -- Collective dephasing coefficient. default: 0.0

collective_pumping

float -- Collective pumping coefficient. default: 0.0

M

dict -- A nested dictionary of the structure {row: {col: val}} which holds non zero elements of the matrix M

calculate_j_m(dicke_row, dicke_col)

Get the value of j and m for the particular Dicke space element.

Parameters:dicke_col (dicke_row,) -- The row and column from the Dicke space matrix Returns:j, m -- The j and m values. Return type:float

calculate_k(dicke_row, dicke_col)

Get k value from the current row and column element in the Dicke space.

Parameters:dicke_col (dicke_row,) -- The row and column from the Dicke space matrix Returns:k -- The row index for the matrix M for given Dicke space element Return type:int

coefficient_matrix()

Generate the matrix M governing the dynamics.

If the initial density matrix and the Hamiltonian is diagonal, the evolution of the system is given by the simple ODE: dp/dt = Mp.

isdicke(dicke_row, dicke_col)

Check if an element in a matrix is a valid element in the Dicke space. Dicke row: j value index. Dicke column: m value index. The function returns True if the element exists in the Dicke space and False otherwise.

Parameters:dicke_col (dicke_row,) -- Index of the element in Dicke space which needs to be checked

solve(rho0, tlist, options=None)

Solve the ODE for the evolution of diagonal states and Hamiltonians.

tau1(j, m)

Calculate tau1 for value of j and m.

tau2(j, m)

Calculate tau2 for given j and m

tau3(j, m)

Calculate tau3.

tau4(j, m)

Calculate tau4.

tau5(j, m)

Calculate tau5.

tau6(j, m)

Calculate tau6.

tau7(j, m)

Calculate tau7.

tau8(j, m)

Calculate tau8.

tau9(j, m)

Calculate tau9.

tau_valid(dicke_row, dicke_col)

Find the Tau functions which are valid for this value of (dicke_row, dicke_col) given the number of TLS. This calculates the valid tau values and reurns a dictionary specifying the tau function name and the value.

Parameters:dicke_col (dicke_row,) -- Index of the element in Dicke space which needs to be checked. Returns:taus -- A dictionary of key, val as {tau: value} consisting of the valid taus for this row and column of the Dicke space element. Return type:dict

piqs.dicke.am(j, m)

Calculate the operator am used later.

The action of ap is given by: J_{-}|j, m> = A_{-}(jm)|j, m-1>

Parameters:m (j,) -- The value for j and m in the dicke basis |j, m>. Returns:a_minus -- The value of a_minus. Return type:float

piqs.dicke.ap(j, m)

Calculate the operator ap used later.

The action of ap is given by: J_{+}|j, m> = A_{+}(jm)|j, m+1>

Parameters:m (j,) -- The value for j and m in the dicke basis |j,m>. Returns:a_plus -- The value of a_plus. Return type:float

piqs.dicke.block_matrix(N)

Construct the block-diagonal matrix for the Dicke basis.

Parameters:N (int) -- Number of two-level systems. Returns:block_matr -- A 2D block-diagonal matrix of ones with dimension (nds,nds), where nds is the number of Dicke states for N two-level systems. Return type:ndarray

piqs.dicke.collapse_uncoupled(N, emission=0.0, dephasing=0.0, pumping=0.0, collective_emission=0.0, collective_dephasing=0.0, collective_pumping=0.0)

Create the collapse operators (c_ops) of the Lindbladian in the uncoupled basis.

These operators are in the uncoupled basis of the two-level system (TLS) SU(2) Pauli matrices.

Notes

The collapse operator list can be given to qutip.mesolve. Notice that the operators are placed in a Hilbert space of dimension 2^N. Thus the method is suitable only for small N (of the order of 10).

Parameters:

• N (int) -- The number of two-level systems.
• emission (float) -- Incoherent emission coefficient (also nonradiative emission). default: 0.0
• dephasing (float) -- Local dephasing coefficient. default: 0.0
• pumping (float) -- Incoherent pumping coefficient. default: 0.0
• collective_emission (float) -- Collective (superradiant) emmission coefficient. default: 0.0
• collective_pumping (float) -- Collective pumping coefficient. default: 0.0
• collective_dephasing (float) -- Collective dephasing coefficient. default: 0.0

Returns:

c_ops -- The list of collapse operators as qutip.Qobj for the system.

Return type:

list

piqs.dicke.css(N, x=0.70710678118654746, y=0.70710678118654746, basis='dicke', coordinates='cartesian')

Generate the density matrix of the Coherent Spin State (CSS).

It can be defined as |CSS>= Prod_i^N(a|1>_i + b|0>_i) with a = sin(theta/2), b = exp(1j*phi) * cos(theta/2). The default basis is that of Dicke space |j, m> < j, m'|. The default state is the symmetric CSS, |CSS> = |+>.

Parameters:

• N (int) -- The number of two-level systems.
• y (x,) -- The coefficients of the CSS state.
• basis (str) -- The basis to use. Either "dicke" or "uncoupled".
• coordinates (str) -- Either "cartesian" or "polar". If polar then the coefficients are constructed as sin(x/2), cos(x/2)e^(iy).

Returns:

rho -- The CSS state density matrix.

Return type:

class:qutip.Qobj

piqs.dicke.dicke(N, j, m)

Generate a Dicke state as a pure density matrix in the Dicke basis.

For instance, if the superradiant state is given |j, m> = |1, 0> for N = 2, the state is represented as a density matrix of size (nds, nds) or (4, 4),

0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0

Parameters:

• N (int) -- The number of two-level systems.
• j (float) -- The eigenvalue j of the Dicke state |j, m>.
• m (float) -- The eigenvalue m of the Dicke state |j, m>.

Returns:

rho -- The density matrix.

Return type:

class:qutip.Qobj

piqs.dicke.dicke_basis(N, jmm1=None)

Initialize the density matrix of a Dicke state for several (j, m, m1).

This function can be used to build arbitrary states in the Dicke basis |j, m><j, m1|. We create coefficients for each (j, m, m1) value in the dictionary jmm1. For instance, if we start from the most excited state for N = 2, we have the following state represented as a density matrix of size (nds, nds) or (4, 4).

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

The mapping for the (i, k) index of the density matrix to the |j, m> values is given by the cythonized function jmm1_dictionary.

Parameters:

• N (int) -- The number of two-level systems.
• jmm1 (dict) -- A dictionary of {(j, m, m1): p} that gives a density p for the (j, m, m1) matrix element.

Returns:

rho -- The density matrix in the Dicke basis.

Return type:

class:qutip.Qobj

piqs.dicke.energy_degeneracy(N, m)

Calculate the number of Dicke states with same energy.

The use of the Decimals class allows to explore N > 1000, unlike the built-in function scipy.special.binom

Parameters:

• N (int) -- The number of two-level systems.
• m (float) -- Total spin z-axis projection eigenvalue. This is proportional to the total energy.

Returns:

degeneracy -- The energy degeneracy

Return type:

int

piqs.dicke.excited(N, basis='dicke')

Generate the density matrix for the excited state.

This state is given by |N/2, N/2> in the default Dicke basis. If the argument basis is "uncoupled" then it generates the state in a 2**N dim Hilbert space.

Parameters:

• N (int) -- The number of two-level systems.
• basis (str) -- The basis to use. Either "dicke" or "uncoupled".

Returns:

state -- The excited state density matrix in the requested basis.

Return type:

class:qutip.Qobj

piqs.dicke.ghz(N, basis='dicke')

Generate the density matrix of the GHZ state.

If the argument basis is "uncoupled" then it generates the state in a 2**N dim Hilbert space.

Parameters:

• N (int) -- The number of two-level systems.
• basis (str) -- The basis to use. Either "dicke" or "uncoupled".

Returns:

state -- The GHZ state density matrix in the requested basis.

Return type:

class:qutip.Qobj

piqs.dicke.ground(N, basis='dicke')

Generate the density matrix of the ground state.

This state is given by |N/2, -N/2> in the Dicke basis. If the argument basis is "uncoupled" then it generates the state in a 2**N dim Hilbert space.

Parameters:

• N (int) -- The number of two-level systems.
• basis (str) -- The basis to use. Either "dicke" or "uncoupled"

Returns:

state -- The ground state density matrix in the requested basis.

Return type:

class:qutip.Qobj

piqs.dicke.identity_uncoupled(N)

Generate the identity in a 2**N dimensional Hilbert space.

The identity matrix is formed from the tensor product of N TLSs.

Parameters:N (int) -- The number of two-level systems. Returns:identity -- The identity matrix. Return type:

class:qutip.Qobj

piqs.dicke.isdiagonal(mat)

Check if the input matrix is diagonal

Parameters:mat (ndarray/Qobj) -- A 2D numpy array Returns:diag -- True/False depending on whether the input matrix is diagonal Return type:bool

piqs.dicke.jspin(N, op=None, basis='dicke')

Calculate the list of collective operators of the total algebra.

The Dicke basis |j,m><j,m'| is used by default. Otherwise with "uncoupled" the operators are in a 2^N space.

Parameters:

• N (int) -- Number of two-level systems.
• op (str) -- The operator to return 'x','y','z','+','-'. If no operator given, then output is the list of operators for ['x','y','z'].
• basis (str) -- The basis of the operators - "dicke" or "uncoupled" default: "dicke".

Returns:

j_alg -- A list of qutip.Qobj representing all the operators in the "dicke" or "uncoupled" basis or a single operator requested.

Return type:

list or :class: qutip.Qobj

piqs.dicke.m_degeneracy(N, m)

Calculate the number of Dicke states |j, m> with same energy.

Parameters:

• N (int) -- The number of two-level systems.
• m (float) -- Total spin z-axis projection eigenvalue (proportional to the total energy).

Returns:

degeneracy -- The m-degeneracy.

Return type:

int

piqs.dicke.num_dicke_ladders(N)

Calculate the total number of ladders in the Dicke space.

For a collection of N two-level systems it counts how many different "j" exist or the number of blocks in the block-diagonal matrix.

Parameters:N (int) -- The number of two-level systems. Returns:Nj -- The number of Dicke ladders. Return type:int

piqs.dicke.num_dicke_states(N)

Calculate the number of Dicke states.

Parameters:N (int) -- The number of two-level systems. Returns:nds -- The number of Dicke states. Return type:int

piqs.dicke.num_tls(nds)

Calculate the number of two-level systems.

Parameters:nds (int) -- The number of Dicke states. Returns:N -- The number of two-level systems. Return type:int

piqs.dicke.spin_algebra(N, op=None)

Create the list [sx, sy, sz] with the spin operators.

The operators are constructed for a collection of N two-level systems (TLSs). Each element of the list, i.e., sx, is a vector of qutip.Qobj objects (spin matrices), as it cointains the list of the SU(2) Pauli matrices for the N TLSs. Each TLS operator sx[i], with i = 0, ..., (N-1), is placed in a 2^N-dimensional Hilbert space.

Notes

sx[i] is sigmax()/2 in the composite Hilbert space.

Parameters:N (int) -- The number of two-level systems. Returns:spin_operators -- A list of qutip.Qobj operators - [sx, sy, sz] or the requested operator. Return type:list or :class: qutip.Qobj

piqs.dicke.state_degeneracy(N, j)

Calculate the degeneracy of the Dicke state.

Each state |j, m> includes D(N,j) irreducible representations |j, m,alpha> Uses Decimals to calculate higher numerator and denominators numbers.

Parameters:

• N (int) -- The number of two-level systems.
• j (float) -- Total spin eigenvalue (cooperativity).

Returns:

degeneracy -- The state degeneracy.

Return type:

int

piqs.dicke.superradiant(N, basis='dicke')

Generate the density matrix of the superradiant state.

This state is given by |N/2, 0> or |N/2, 0.5> in the Dicke basis. If the argument basis is "uncoupled" then it generates the state in a 2**N dim Hilbert space.

Parameters:

• N (int) -- The number of two-level systems.
• basis (str) -- The basis to use. Either "dicke" or "uncoupled".

Returns:

state -- The superradiant state density matrix in the requested basis.

Return type:

class:qutip.Qobj

piqs.dicke.tau_column(tau, k, j)

Determine the column index for the non-zero elements of the matrix for a particular row k and the value of j from the Dicke space.

Parameters:

• tau (str) -- The tau function to check for this k and j.
• k (int) -- The row of the matrix M for which the non zero elements have to be calculated.
• j (float) -- The value of j for this row.

Cythonized Dicke module

Cythonized code for permutationally invariant Lindbladian generation

class piqs.cy.dicke.Dicke

A faster Cythonized Dicke state class to build the Lindbladian.

Parameters:

• N (int) -- The number of two-level systems.
• emission (float) -- Incoherent emission coefficient (also nonradiative emission). default: 0.0
• dephasing (float) -- Local dephasing coefficient. default: 0.0
• pumping (float) -- Incoherent pumping coefficient. default: 0.0
• collective_emission (float) -- Collective (superradiant) emmission coefficient. default: 0.0
• collective_pumping (float) -- Collective pumping coefficient. default: 0.0
• collective_dephasing (float) -- Collective dephasing coefficient. default: 0.0

N

int -- The number of two-level systems.

emission

float -- Incoherent emission coefficient (also nonradiative emission). default: 0.0

dephasing

float -- Local dephasing coefficient. default: 0.0

pumping

float -- Incoherent pumping coefficient. default: 0.0

collective_emission

float -- Collective (superradiant) emmission coefficient. default: 0.0

collective_pumping

float -- Collective pumping coefficient. default: 0.0

collective_dephasing

float -- Collective dephasing coefficient. default: 0.0

gamma1()

Calculate gamma1 for value of j, m, m'.

gamma2()

Calculate gamma2 for given j, m, m'.

gamma3()

Calculate gamma3 for given j, m, m'.

gamma4()

Calculate gamma4 for given j, m, m'.

gamma5()

Calculate gamma5 for given j, m, m'.

gamma6()

Calculate gamma6 for given j, m, m'.

gamma7()

Calculate gamma7 for given j, m, m'.

gamma8()

Calculate gamma8 for given j, m, m'.

gamma9()

Calculate gamma9 for given j, m, m'.

lindbladian()

Build the Lindbladian superoperator of the dissipative dynamics as a sparse matrix.

Returns:lindblad_qobj -- The matrix size is (nds**2, nds**2) where nds is the number of Dicke states. Return type:

class:qutip.Qobj

piqs.cy.dicke.get_blocks()

Calculate the number of cumulative elements at each block boundary.

Parameters:N (int) -- The number of two-level systems. Returns:blocks -- An array with the number of cumulative elements at the boundary of each block. Return type:np.ndarray

piqs.cy.dicke.get_index()

Get the index in the density matrix for this j, m, m1 value.

Parameters:

• N (int) -- The number of two-level systems.
• m, m1 (j,) -- The j, m, m1 values.
• blocks (np.ndarray) -- An 1D array with the number of cumulative elements at the boundary of each block.

Returns:

mvals -- The m values for given j.

Return type:

array

piqs.cy.dicke.j_min()

Calculate the minimum value of j for given N.

Parameters:N (int) -- Number of two-level systems. Returns:jmin -- The minimum value of j for odd or even number of two level systems. Return type:float

piqs.cy.dicke.j_vals()

Get the valid values of j for given N.

Parameters:N (int) -- The number of two-level systems. Returns:jvals -- The j values for given N as a 1D array. Return type:np.ndarray

piqs.cy.dicke.jmm1_dictionary()

Get the index in the density matrix for this j, m, m1 value.

The (j, m, m1) values are mapped to the (i, k) index of a block diagonal matrix which has the structure to capture the permutationally symmetric part of the density matrix. For each (j, m, m1) value, first we get the block by using the "j" value and then the addition in the row/column due to the m and m1 is determined. Four dictionaries are returned giving a map from the (j, m, m1) values to (i, k), the inverse map, a flattened map and the inverse of the flattened map.

piqs.cy.dicke.m_vals()

Get all the possible values of m or m1 for given j.

Parameters:N (int) -- The number of two-level systems. Returns:mvals -- The m values for given j as a 1D array. Return type:np.ndarray