# This file is part of QuTiP: Quantum Toolbox in Python.
#
# Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# 1. Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
# of its contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
###############################################################################
"""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.
"""
# Authors: Nathan Shammah, Shahnawaz Ahmed
# Contact: nathan.shammah@gmail.com, shahnawaz.ahmed95@gmail.com
from math import factorial
from decimal import Decimal
import numpy as np
from scipy.integrate import odeint
from scipy import constants
from scipy.sparse import dok_matrix, block_diag, lil_matrix
from qutip.solver import Options, Result
from qutip import (Qobj, spre, spost, tensor, identity, ket2dm,
vector_to_operator)
from qutip import sigmax, sigmay, sigmaz, sigmap, sigmam
from qutip.cy.piqs import Dicke as _Dicke
from qutip.cy.piqs import (jmm1_dictionary, _num_dicke_states,
_num_dicke_ladders, get_blocks, j_min,
m_vals, j_vals)
# Functions necessary to generate the Lindbladian/Liouvillian
[docs]def num_dicke_states(N):
"""Calculate the number of Dicke states.
Parameters
----------
N: int
The number of two-level systems.
Returns
-------
nds: int
The number of Dicke states.
"""
return _num_dicke_states(N)
[docs]def 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: int
The number of Dicke ladders.
"""
return _num_dicke_ladders(N)
[docs]def num_tls(nds):
"""Calculate the number of two-level systems.
Parameters
----------
nds: int
The number of Dicke states.
Returns
-------
N: int
The number of two-level systems.
"""
if np.sqrt(nds).is_integer():
# N is even
N = 2*(np.sqrt(nds)-1)
else:
# N is odd
N = 2*(np.sqrt(nds + 1/4)-1)
return int(N)
[docs]def isdiagonal(mat):
"""
Check if the input matrix is diagonal.
Parameters
==========
mat: ndarray/Qobj
A 2D numpy array
Returns
=======
diag: bool
True/False depending on whether the input matrix is diagonal.
"""
if isinstance(mat, Qobj):
mat = mat.full()
return np.all(mat == np.diag(np.diagonal(mat)))
[docs]class Dicke(object):
"""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: :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` 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
Attributes
----------
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).
"""
def __init__(self, N, hamiltonian=None,
emission=0., dephasing=0., pumping=0.,
collective_emission=0., collective_dephasing=0.,
collective_pumping=0.):
self.N = N
self.hamiltonian = hamiltonian
self.emission = emission
self.dephasing = dephasing
self.pumping = pumping
self.collective_emission = collective_emission
self.collective_dephasing = collective_dephasing
self.collective_pumping = collective_pumping
self.nds = num_dicke_states(self.N)
self.dshape = (num_dicke_states(self.N), num_dicke_states(self.N))
def __repr__(self):
"""Print the current parameters of the system."""
string = []
string.append("N = {}".format(self.N))
string.append("Hilbert space dim = {}".format(self.dshape))
string.append("Number of Dicke states = {}".format(self.nds))
string.append("Liouvillian space dim = {}".format((self.nds**2,
self.nds**2)))
if self.emission != 0:
string.append("emission = {}".format(self.emission))
if self.dephasing != 0:
string.append("dephasing = {}".format(self.dephasing))
if self.pumping != 0:
string.append("pumping = {}".format(self.pumping))
if self.collective_emission != 0:
string.append(
"collective_emission = {}".format(self.collective_emission))
if self.collective_dephasing != 0:
string.append(
"collective_dephasing = {}".format(self.collective_dephasing))
if self.collective_pumping != 0:
string.append(
"collective_pumping = {}".format(self.collective_pumping))
return "\n".join(string)
[docs] def lindbladian(self):
"""Build the Lindbladian superoperator of the dissipative dynamics.
Returns
-------
lindbladian: :class: qutip.Qobj
The Lindbladian matrix as a `qutip.Qobj`.
"""
cythonized_dicke = _Dicke(int(self.N),
float(self.emission),
float(self.dephasing),
float(self.pumping),
float(self.collective_emission),
float(self.collective_dephasing),
float(self.collective_pumping))
return cythonized_dicke.lindbladian()
[docs] def liouvillian(self):
"""Build the total Liouvillian using the Dicke basis.
Returns
-------
liouv: :class: qutip.Qobj
The Liouvillian matrix for the system.
"""
lindblad = self.lindbladian()
if self.hamiltonian is None:
liouv = lindblad
else:
hamiltonian = self.hamiltonian
hamiltonian_superoperator = - 1j * \
spre(hamiltonian) + 1j * spost(hamiltonian)
liouv = lindblad + hamiltonian_superoperator
return liouv
[docs] def pisolve(self, initial_state, tlist, options=None):
"""
Solve for diagonal Hamiltonians and initial states faster.
Parameters
==========
initial_state: :class: qutip.Qobj
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: :class: qutip.solver.Options
The options for the solver.
Returns
=======
result: list
A dictionary of the type `qutip.solver.Result` which holds the
results of the evolution.
"""
if isdiagonal(initial_state) == False:
msg = "`pisolve` requires a diagonal initial density matrix. "
msg += "In general construct the Liouvillian using "
msg += "`piqs.liouvillian` and use qutip.mesolve."
raise ValueError(msg)
if self.hamiltonian and isdiagonal(self.hamiltonian) == False:
msg = "`pisolve` should only be used for diagonal Hamiltonians. "
msg += "Construct the Liouvillian using `piqs.liouvillian` and"
msg += " use `qutip.mesolve`."
raise ValueError(msg)
if initial_state.full().shape != self.dshape:
msg = "Initial density matrix should be diagonal."
raise ValueError(msg)
pim = Pim(self.N, self.emission, self.dephasing, self.pumping,
self.collective_emission, self.collective_pumping,
self.collective_dephasing)
result = pim.solve(initial_state, tlist, options=None)
return result
[docs] def prune_eigenstates(self, 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: list
The list with the correct eigenvalues and eigenvectors of the
Liouvillian.
"""
liouvillian_eigenstates = liouvillian.eigenstates()
N = self.N
block_mat = block_matrix(N)
nnz_tuple_bm = [(i, j) for i, j in zip(*block_mat.nonzero())]
# 0. Create a copy of the eigenvalues to approximate values
eig_val, eig_vec = liouvillian_eigenstates
tol = 10
eig_val_round = np.round(eig_val, tol)
# 2. Use 'block_matrix(N)' to remove eigenvectors with matrix
# elements
# outside of the block matrix.
forbidden_eig_index = []
for k in range(0, len(eig_vec)):
dm = vector_to_operator(eig_vec[k])
nnz_tuple = [(i, j) for i, j in zip(*dm.data.nonzero())]
for i in nnz_tuple:
if i not in nnz_tuple_bm:
if np.round(dm[i], tol) != 0:
forbidden_eig_index.append(k)
forbidden_eig_index = np.array(list(set(forbidden_eig_index)))
# 3. Remove the forbidden eigenvalues and eigenvectors.
correct_eig_val = np.delete(eig_val, forbidden_eig_index)
correct_eig_vec = np.delete(eig_vec, forbidden_eig_index)
correct_eigenstates = correct_eig_val, correct_eig_vec
return correct_eigenstates
[docs] def c_ops(self):
"""Build collapse operators in the full Hilbert space 2^N.
Returns
-------
c_ops_list: list
The list with the collapse operators in the 2^N Hilbert space.
"""
ce = self.collective_emission
cd = self.collective_dephasing
cp = self.collective_pumping
c_ops_list = collapse_uncoupled(N=self.N,
emission=self.emission,
dephasing=self.dephasing,
pumping=self.pumping,
collective_emission=ce,
collective_dephasing=cd,
collective_pumping=cp)
return c_ops_list
[docs] def coefficient_matrix(self):
"""Build coefficient matrix for ODE for a diagonal problem.
Returns
-------
M: ndarray
The matrix M of the coefficients for the ODE dp/dt = Mp.
p is the vector of the diagonal matrix elements
of the density matrix rho in the Dicke basis.
"""
diagonal_system = Pim(N=self.N,
emission=self.emission,
dephasing=self.dephasing,
pumping=self.pumping,
collective_emission=self.collective_emission,
collective_dephasing=self.collective_dephasing,
collective_pumping=self.collective_pumping)
coef_matrix = diagonal_system.coefficient_matrix()
return coef_matrix
# Utility functions for properties of the Dicke space
def 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: int
The energy degeneracy
"""
numerator = Decimal(factorial(N))
d1 = Decimal(factorial(N/2 + m))
d2 = Decimal(factorial(N/2 - m))
degeneracy = numerator/(d1 * d2)
return int(degeneracy)
[docs]def state_degeneracy(N, j):
"""Calculate the degeneracy of the Dicke state.
Each state :math:`|j, m\\rangle` includes D(N,j) irreducible
representations :math:`|j, m, \\alpha\\rangle`.
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: int
The state degeneracy.
"""
if j < 0:
raise ValueError("j value should be >= 0")
numerator = Decimal(factorial(N)) * Decimal(2*j + 1)
denominator_1 = Decimal(factorial(N/2 + j + 1))
denominator_2 = Decimal(factorial(N/2 - j))
degeneracy = numerator/(denominator_1 * denominator_2)
degeneracy = int(np.round(float(degeneracy)))
return degeneracy
[docs]def m_degeneracy(N, m):
"""Calculate the number of Dicke states :math:`|j, m\\rangle` 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: int
The m-degeneracy.
"""
jvals = j_vals(N)
maxj = np.max(jvals)
if m < -maxj:
e = "m value is incorrect for this N."
e += " Minimum m value can be {}".format(-maxj)
raise ValueError(e)
degeneracy = N/2 + 1 - abs(m)
return int(degeneracy)
[docs]def ap(j, m):
"""Calculate the coefficient `ap` by applying J_+ |j, m>.
The action of ap is given by:
:math:`J_{+}|j, m\\rangle = A_{+}(j, m)|j, m+1\\rangle`
Parameters
----------
j, m: float
The value for j and m in the dicke basis |j,m>.
Returns
-------
a_plus: float
The value of :math:`a_{+}`.
"""
a_plus = np.sqrt((j-m) * (j+m+1))
return a_plus
[docs]def 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
----------
j: float
The value for j.
m: float
The value for m.
Returns
-------
a_minus: float
The value of :math:`a_{-}`.
"""
a_minus = np.sqrt((j+m) * (j-m+1))
return a_minus
[docs]def 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 :math:`2^N`-dimensional Hilbert space.
Notes
-----
sx[i] is :math:`\\frac{\\sigma_x}{2}` in the composite Hilbert space.
Parameters
----------
N: int
The number of two-level systems.
Returns
-------
spin_operators: list or :class: qutip.Qobj
A list of `qutip.Qobj` operators - [sx, sy, sz] or the
requested operator.
"""
# 1. Define N TLS spin-1/2 matrices in the uncoupled basis
N = int(N)
sx = [0 for i in range(N)]
sy = [0 for i in range(N)]
sz = [0 for i in range(N)]
sp = [0 for i in range(N)]
sm = [0 for i in range(N)]
sx[0] = 0.5 * sigmax()
sy[0] = 0.5 * sigmay()
sz[0] = 0.5 * sigmaz()
sp[0] = sigmap()
sm[0] = sigmam()
# 2. Place operators in total Hilbert space
for k in range(N - 1):
sx[0] = tensor(sx[0], identity(2))
sy[0] = tensor(sy[0], identity(2))
sz[0] = tensor(sz[0], identity(2))
sp[0] = tensor(sp[0], identity(2))
sm[0] = tensor(sm[0], identity(2))
# 3. Cyclic sequence to create all N operators
a = [i for i in range(N)]
b = [[a[i - i2] for i in range(N)] for i2 in range(N)]
# 4. Create N operators
for i in range(1, N):
sx[i] = sx[0].permute(b[i])
sy[i] = sy[0].permute(b[i])
sz[i] = sz[0].permute(b[i])
sp[i] = sp[0].permute(b[i])
sm[i] = sm[0].permute(b[i])
spin_operators = [sx, sy, sz]
if not op:
return spin_operators
elif op == 'x':
return sx
elif op == 'y':
return sy
elif op == 'z':
return sz
elif op == '+':
return sp
elif op == '-':
return sm
else:
raise TypeError('Invalid type')
def _jspin_uncoupled(N, op=None):
"""Construct the the collective spin algebra in the uncoupled basis.
jx, jy, jz, jp, jm are constructed in the uncoupled basis of the
two-level system (TLS). Each collective operator is placed in a
Hilbert space of dimension 2^N.
Parameters
----------
N: int
The 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',].
Returns
-------
collective_operators: list or :class: qutip.Qobj
A list of `qutip.Qobj` representing all the operators in
uncoupled" basis or a single operator requested.
"""
# 1. Define N TLS spin-1/2 matrices in the uncoupled basis
N = int(N)
sx, sy, sz = spin_algebra(N)
sp, sm = spin_algebra(N, "+"), spin_algebra(N, "-")
jx = sum(sx)
jy = sum(sy)
jz = sum(sz)
jp = sum(sp)
jm = sum(sm)
collective_operators = [jx, jy, jz]
if not op:
return collective_operators
elif op == 'x':
return jx
elif op == 'y':
return jy
elif op == 'z':
return jz
elif op == '+':
return jp
elif op == '-':
return jm
else:
raise TypeError('Invalid type')
[docs]def jspin(N, op=None, basis="dicke"):
"""
Calculate the list of collective operators of the total algebra.
The Dicke basis :math:`|j,m\\rangle\\langle j,m'|` is used by
default. Otherwise with "uncoupled" the operators are in a
:math:`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: list or :class: qutip.Qobj
A list of `qutip.Qobj` representing all the operators in
the "dicke" or "uncoupled" basis or a single operator requested.
"""
if basis == "uncoupled":
return _jspin_uncoupled(N, op)
nds = num_dicke_states(N)
num_ladders = num_dicke_ladders(N)
jz_operator = dok_matrix((nds, nds))
jp_operator = dok_matrix((nds, nds))
jm_operator = dok_matrix((nds, nds))
s = 0
for k in range(0, num_ladders):
j = 0.5 * N - k
mmax = int(2*j + 1)
for i in range(0, mmax):
m = j-i
jz_operator[s, s] = m
if (s+1) in range(0, nds):
jp_operator[s, s+1] = ap(j, m-1)
if (s-1) in range(0, nds):
jm_operator[s, s-1] = am(j, m+1)
s = s+1
jx_operator = 1/2 * (jp_operator+jm_operator)
jy_operator = 1j/2 * (jm_operator-jp_operator)
jx = Qobj(jx_operator)
jy = Qobj(jy_operator)
jz = Qobj(jz_operator)
jp = Qobj(jp_operator)
jm = Qobj(jm_operator)
if not op:
return [jx, jy, jz]
if op == '+':
return jp
elif op == '-':
return jm
elif op == 'x':
return jx
elif op == 'y':
return jy
elif op == 'z':
return jz
else:
raise TypeError('Invalid type')
[docs]def collapse_uncoupled(N, emission=0., dephasing=0., pumping=0.,
collective_emission=0., collective_dephasing=0.,
collective_pumping=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 :math:`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: list
The list of collapse operators as `qutip.Qobj` for the system.
"""
N = int(N)
if N > 10:
msg = "N > 10. dim(H) = 2^N. "
msg += "Better use `piqs.lindbladian` to reduce Hilbert space "
msg += "dimension and exploit permutational symmetry."
raise Warning(msg)
[sx, sy, sz] = spin_algebra(N)
sp, sm = spin_algebra(N, "+"), spin_algebra(N, "-")
[jx, jy, jz] = jspin(N, basis="uncoupled")
jp, jm = (jspin(N, "+", basis = "uncoupled"),
jspin(N, "-", basis="uncoupled"))
c_ops = []
if emission != 0:
for i in range(0, N):
c_ops.append(np.sqrt(emission) * sm[i])
if dephasing != 0:
for i in range(0, N):
c_ops.append(np.sqrt(dephasing) * sz[i])
if pumping != 0:
for i in range(0, N):
c_ops.append(np.sqrt(pumping) * sp[i])
if collective_emission != 0:
c_ops.append(np.sqrt(collective_emission) * jm)
if collective_dephasing != 0:
c_ops.append(np.sqrt(collective_dephasing) * jz)
if collective_pumping != 0:
c_ops.append(np.sqrt(collective_pumping) * jp)
return c_ops
# State definitions in the Dicke basis with an option for basis transformation
[docs]def 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
:math:`|j, m\\rangle \\langle j, m^{\\prime}|`. We create coefficients for each
(j, m, m1) value in the dictionary jmm1. The mapping for the (i, k)
index of the density matrix to the |j, m> values is given by the
cythonized function `jmm1_dictionary`. A density matrix is created from
the given dictionary of coefficients for each (j, m, m1).
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: :class: qutip.Qobj
The density matrix in the Dicke basis.
"""
if jmm1 is None:
msg = "Please specify the jmm1 values as a dictionary"
msg += "or use the `excited(N)` function to create an"
msg += "excited state where jmm1 = {(N/2, N/2, N/2): 1}"
raise AttributeError(msg)
nds = _num_dicke_states(N)
rho = np.zeros((nds, nds))
jmm1_dict = jmm1_dictionary(N)[1]
for key in jmm1:
i, k = jmm1_dict[key]
rho[i, k] = jmm1[key]
return Qobj(rho)
[docs]def dicke(N, j, m):
"""
Generate a Dicke state as a pure density matrix in the Dicke basis.
For instance, the superradiant state given by
:math:`|j, m\\rangle = |1, 0\\rangle` for N = 2,
and the state is represented as a density matrix of size (nds, nds) or
(4, 4), with the (1, 1) element set to 1.
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: :class: qutip.Qobj
The density matrix.
"""
nds = num_dicke_states(N)
rho = np.zeros((nds, nds))
jmm1_dict = jmm1_dictionary(N)[1]
i, k = jmm1_dict[(j, m, m)]
rho[i, k] = 1.
return Qobj(rho)
# Uncoupled states in the full Hilbert space. These are returned with the
# choice of the keyword argument `basis="uncoupled"` in the state functions.
def _uncoupled_excited(N):
"""
Generate the density matrix of the excited Dicke state in the full
:math:`2^N` dimensional Hilbert space.
Parameters
----------
N: int
The number of two-level systems.
Returns
-------
psi0: :class: qutip.Qobj
The density matrix for the excited state in the uncoupled basis.
"""
N = int(N)
jz = jspin(N, "z", basis="uncoupled")
en, vn = jz.eigenstates()
psi0 = vn[2**N - 1]
return ket2dm(psi0)
def _uncoupled_superradiant(N):
"""
Generate the density matrix of a superradiant state in the full
:math:`2^N`-dimensional Hilbert space.
Parameters
----------
N: int
The number of two-level systems.
Returns
-------
psi0: :class: qutip.Qobj
The density matrix for the superradiant state in the full Hilbert
space.
"""
N = int(N)
jz = jspin(N, "z", basis="uncoupled")
en, vn = jz.eigenstates()
psi0 = vn[2**N - (N+1)]
return ket2dm(psi0)
def _uncoupled_ground(N):
"""
Generate the density matrix of the ground state in the full\
:math:`2^N`-dimensional Hilbert space.
Parameters
----------
N: int
The number of two-level systems.
Returns
-------
psi0: :class: qutip.Qobj
The density matrix for the ground state in the full Hilbert space.
"""
N = int(N)
jz = jspin(N, "z", basis="uncoupled")
en, vn = jz.eigenstates()
psi0 = vn[0]
return ket2dm(psi0)
def _uncoupled_ghz(N):
"""
Generate the density matrix of the GHZ state in the full 2^N
dimensional Hilbert space.
Parameters
----------
N: int
The number of two-level systems.
Returns
-------
ghz: :class: qutip.Qobj
The density matrix for the GHZ state in the full Hilbert space.
"""
N = int(N)
rho = np.zeros((2**N, 2**N))
rho[0, 0] = 1/2
rho[2**N - 1, 0] = 1/2
rho[0, 2**N - 1] = 1/2
rho[2**N - 1, 2**N - 1] = 1/2
spin_dim = [2 for i in range(0, N)]
spins_dims = list((spin_dim, spin_dim))
rho = Qobj(rho, dims=spins_dims)
return rho
def _uncoupled_css(N, a, b):
"""
Generate the density matrix of the CSS state in the full 2^N
dimensional Hilbert space.
The CSS states are non-entangled states given by
:math:`|a, b\\rangle = \\prod_i (a|1\\rangle_i + b|0\\rangle_i)`.
Parameters
----------
N: int
The number of two-level systems.
a: complex
The coefficient of the :math:`|1_i\rangle` state.
b: complex
The coefficient of the :math:`|0_i\rangle` state.
Returns
-------
css: :class: qutip.Qobj
The density matrix for the CSS state in the full Hilbert space.
"""
N = int(N)
# 1. Define i_th factorized density matrix in the uncoupled basis
rho_i = np.zeros((2, 2), dtype=complex)
rho_i[0, 0] = a * np.conj(a)
rho_i[1, 1] = b * np.conj(b)
rho_i[0, 1] = a * np.conj(a)
rho_i[1, 0] = b * np.conj(b)
rho_i = Qobj(rho_i)
rho = [0 for i in range(N)]
rho[0] = rho_i
# 2. Place single-two-level-system density matrices in total Hilbert space
for k in range(N - 1):
rho[0] = tensor(rho[0], identity(2))
# 3. Cyclic sequence to create all N factorized density matrices
# |CSS>_i<CSS|_i
a = [i for i in range(N)]
b = [[a[i - i2] for i in range(N)] for i2 in range(N)]
# 4. Create all other N-1 factorized density matrices
# |+><+| = Prod_(i=1)^N |CSS>_i<CSS|_i
for i in range(1, N):
rho[i] = rho[0].permute(b[i])
identity_i = Qobj(np.eye(2**N), dims=rho[0].dims, shape=rho[0].shape)
rho_tot = identity_i
for i in range(0, N):
rho_tot = rho_tot * rho[i]
return rho_tot
[docs]def 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: :class: qutip.Qobj
The excited state density matrix in the requested basis.
"""
if basis == "uncoupled":
state = _uncoupled_excited(N)
return state
jmm1 = {(N/2, N/2, N/2): 1}
return dicke_basis(N, jmm1)
[docs]def 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: :class: qutip.Qobj
The superradiant state density matrix in the requested basis.
"""
if basis == "uncoupled":
state = _uncoupled_superradiant(N)
return state
if N % 2 == 0:
jmm1 = {(N/2, 0, 0): 1.}
return dicke_basis(N, jmm1)
else:
jmm1 = {(N/2, 0.5, 0.5): 1.}
return dicke_basis(N, jmm1)
[docs]def css(N, x=1/np.sqrt(2), y=1/np.sqrt(2),
basis="dicke", coordinates="cartesian"):
"""
Generate the density matrix of the Coherent Spin State (CSS).
It can be defined as,
:math:`|CSS \\rangle = \\prod_i^N(a|1\\rangle_i + b|0\\rangle_i)`
with :math:`a = sin(\\frac{\\theta}{2})`,
:math:`b = e^{i \\phi}\\cos(\\frac{\\theta}{2})`.
The default basis is that of Dicke space
:math:`|j, m\\rangle \\langle j, m'|`.
The default state is the symmetric CSS,
:math:`|CSS\\rangle = |+\\rangle`.
Parameters
----------
N: int
The number of two-level systems.
x, y: float
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: :class: qutip.Qobj
The CSS state density matrix.
"""
if coordinates == "polar":
a = np.cos(0.5 * x) * np.exp(1j * y)
b = np.sin(0.5 * x)
else:
a = x
b = y
if basis == "uncoupled":
return _uncoupled_css(N, a, b)
nds = num_dicke_states(N)
num_ladders = num_dicke_ladders(N)
rho = dok_matrix((nds, nds))
# loop in the allowed matrix elements
jmm1_dict = jmm1_dictionary(N)[1]
j = 0.5*N
mmax = int(2*j + 1)
for i in range(0, mmax):
m = j-i
psi_m = np.sqrt(float(energy_degeneracy(N, m))) * \
a**(N*0.5 + m) * b**(N*0.5 - m)
for i1 in range(0, mmax):
m1 = j - i1
row_column = jmm1_dict[(j, m, m1)]
psi_m1 = np.sqrt(float(energy_degeneracy(N, m1))) * \
np.conj(a)**(N*0.5 + m1) * np.conj(b)**(N*0.5 - m1)
rho[row_column] = psi_m*psi_m1
return Qobj(rho)
[docs]def 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 :math:`2^N`-dimensional Hilbert space.
Parameters
----------
N: int
The number of two-level systems.
basis: str
The basis to use. Either "dicke" or "uncoupled".
Returns
-------
state: :class: qutip.Qobj
The GHZ state density matrix in the requested basis.
"""
if basis == "uncoupled":
return _uncoupled_ghz(N)
nds = _num_dicke_states(N)
rho = dok_matrix((nds, nds))
rho[0, 0] = 1/2
rho[N, N] = 1/2
rho[N, 0] = 1/2
rho[0, N] = 1/2
return Qobj(rho)
[docs]def 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
:math:`2^N`-dimensional Hilbert space.
Parameters
----------
N: int
The number of two-level systems.
basis: str
The basis to use. Either "dicke" or "uncoupled"
Returns
-------
state: :class: qutip.Qobj
The ground state density matrix in the requested basis.
"""
if basis == "uncoupled":
state = _uncoupled_ground(N)
return state
nds = _num_dicke_states(N)
rho = dok_matrix((nds, nds))
rho[N, N] = 1
return Qobj(rho)
def identity_uncoupled(N):
"""
Generate the identity in a :math:`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: :class: qutip.Qobj
The identity matrix.
"""
N = int(N)
rho = np.zeros((2**N, 2**N))
for i in range(0, 2**N):
rho[i, i] = 1
spin_dim = [2 for i in range(0, N)]
spins_dims = list((spin_dim, spin_dim))
identity = Qobj(rho, dims=spins_dims)
return identity
def block_matrix(N):
"""Construct the block-diagonal matrix for the Dicke basis.
Parameters
----------
N: int
Number of two-level systems.
Returns
-------
block_matr: ndarray
A 2D block-diagonal matrix of ones with dimension (nds,nds),
where nds is the number of Dicke states for N two-level
systems.
"""
nds = _num_dicke_states(N)
# create a list with the sizes of the blocks, in order
blocks_dimensions = int(N/2 + 1 - 0.5*(N % 2))
blocks_list = [(2 * (i+1 * (N % 2)) + 1*((N+1) % 2))
for i in range(blocks_dimensions)]
blocks_list = np.flip(blocks_list, 0)
# create a list with each block matrix as element
square_blocks = []
k = 0
for i in blocks_list:
square_blocks.append(np.ones((i, i)))
k = k + 1
return block_diag(square_blocks)
# ============================================================================
# Adding a faster version to make a Permutational Invariant matrix
# ============================================================================
def 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.
"""
# In the notes, we indexed from k = 1, here we do it from k = 0
k = k + 1
mapping = {"tau3": k-(2 * j + 3),
"tau2": k-1,
"tau4": k+(2 * j - 1),
"tau5": k-(2 * j + 2),
"tau1": k,
"tau6": k+(2 * j),
"tau7": k-(2 * j + 1),
"tau8": k+1,
"tau9": k+(2 * j + 1)}
# we need to decrement k again as indexing is from 0
return int(mapping[tau] - 1)
[docs]class Pim(object):
"""
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
Attributes
----------
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
"""
def __init__(self, N, emission=0., dephasing=0, pumping=0,
collective_emission=0, collective_pumping=0,
collective_dephasing=0):
self.N = N
self.emission = emission
self.dephasing = dephasing
self.pumping = pumping
self.collective_pumping = collective_pumping
self.collective_dephasing = collective_dephasing
self.collective_emission = collective_emission
self.M = {}
[docs] def isdicke(self, 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_row, dicke_col : int
Index of the element in Dicke space which needs to be checked
"""
rows = self.N + 1
cols = 0
if (self.N % 2) == 0:
cols = int(self.N/2 + 1)
else:
cols = int(self.N/2 + 1/2)
if (dicke_row > rows) or (dicke_row < 0):
return (False)
if (dicke_col > cols) or (dicke_col < 0):
return (False)
if (dicke_row < int(rows/2)) and (dicke_col > dicke_row):
return False
if (dicke_row >= int(rows/2)) and (rows - dicke_row <= dicke_col):
return False
else:
return True
[docs] def tau_valid(self, 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_row, dicke_col : int
Index of the element in Dicke space which needs to be checked.
Returns
-------
taus: dict
A dictionary of key, val as {tau: value} consisting of the valid
taus for this row and column of the Dicke space element.
"""
tau_functions = [self.tau3, self.tau2, self.tau4,
self.tau5, self.tau1, self.tau6,
self.tau7, self.tau8, self.tau9]
N = self.N
if self.isdicke(dicke_row, dicke_col) is False:
return False
# The 3x3 sub matrix surrounding the Dicke space element to
# run the tau functions
indices = [(dicke_row + x, dicke_col + y) for x in range(-1, 2)
for y in range(-1, 2)]
taus = {}
for idx, tau in zip(indices, tau_functions):
if self.isdicke(idx[0], idx[1]):
j, m = self.calculate_j_m(idx[0], idx[1])
taus[tau.__name__] = tau(j, m)
return taus
[docs] def calculate_j_m(self, dicke_row, dicke_col):
"""
Get the value of j and m for the particular Dicke space element.
Parameters
----------
dicke_row, dicke_col: int
The row and column from the Dicke space matrix
Returns
-------
j, m: float
The j and m values.
"""
N = self.N
j = N/2 - dicke_col
m = N/2 - dicke_row
return(j, m)
[docs] def calculate_k(self, dicke_row, dicke_col):
"""
Get k value from the current row and column element in the Dicke space.
Parameters
----------
dicke_row, dicke_col: int
The row and column from the Dicke space matrix.
Returns
-------
k: int
The row index for the matrix M for given Dicke space
element.
"""
N = self.N
if dicke_row == 0:
k = dicke_col
else:
k = int(((dicke_col)/2) * (2 * (N + 1) - 2 * (dicke_col - 1)) +
(dicke_row - (dicke_col)))
return k
[docs] def coefficient_matrix(self):
"""
Generate the matrix M governing the dynamics for diagonal cases.
If the initial density matrix and the Hamiltonian is diagonal, the
evolution of the system is given by the simple ODE: dp/dt = Mp.
"""
N = self.N
nds = num_dicke_states(N)
rows = self.N + 1
cols = 0
sparse_M = lil_matrix((nds, nds), dtype=float)
if (self.N % 2) == 0:
cols = int(self.N/2 + 1)
else:
cols = int(self.N/2 + 1/2)
for (dicke_row, dicke_col) in np.ndindex(rows, cols):
if self.isdicke(dicke_row, dicke_col):
k = int(self.calculate_k(dicke_row, dicke_col))
row = {}
taus = self.tau_valid(dicke_row, dicke_col)
for tau in taus:
j, m = self.calculate_j_m(dicke_row, dicke_col)
current_col = tau_column(tau, k, j)
sparse_M[k, int(current_col)] = taus[tau]
return sparse_M.tocsr()
[docs] def solve(self, rho0, tlist, options=None):
"""
Solve the ODE for the evolution of diagonal states and Hamiltonians.
"""
if options is None:
options = Options()
output = Result()
output.solver = "pisolve"
output.times = tlist
output.states = []
output.states.append(Qobj(rho0))
rhs_generate = lambda y, tt, M: M.dot(y)
rho0_flat = np.diag(np.real(rho0.full()))
L = self.coefficient_matrix()
rho_t = odeint(rhs_generate, rho0_flat, tlist, args=(L,))
for r in rho_t[1:]:
diag = np.diag(r)
output.states.append(Qobj(diag))
return output
[docs] def tau1(self, j, m):
"""
Calculate coefficient matrix element relative to (j, m, m).
"""
yS = self.collective_emission
yL = self.emission
yD = self.dephasing
yP = self.pumping
yCP = self.collective_pumping
N = float(self.N)
spontaneous = yS * (1 + j - m) * (j + m)
losses = yL * (N/2 + m)
pump = yP * (N/2 - m)
collective_pump = yCP * (1 + j + m) * (j - m)
if j == 0:
dephase = yD * N/4
else:
dephase = yD * (N/4 - m**2 * ((1 + N/2)/(2 * j * (j+1))))
t1 = spontaneous + losses + pump + dephase + collective_pump
return -t1
[docs] def tau2(self, j, m):
"""
Calculate coefficient matrix element relative to (j, m+1, m+1).
"""
yS = self.collective_emission
yL = self.emission
N = float(self.N)
spontaneous = yS * (1 + j - m) * (j + m)
losses = yL * (((N/2 + 1) * (j - m + 1) * (j + m))/(2 * j * (j+1)))
t2 = spontaneous + losses
return t2
[docs] def tau3(self, j, m):
"""
Calculate coefficient matrix element relative to (j+1, m+1, m+1).
"""
yL = self.emission
N = float(self.N)
num = (j + m - 1) * (j + m) * (j + 1 + N/2)
den = 2 * j * (2 * j + 1)
t3 = yL * (num/den)
return t3
[docs] def tau4(self, j, m):
"""
Calculate coefficient matrix element relative to (j-1, m+1, m+1).
"""
yL = self.emission
N = float(self.N)
num = (j - m + 1) * (j - m + 2) * (N/2 - j)
den = 2 * (j + 1) * (2 * j + 1)
t4 = yL * (num/den)
return t4
[docs] def tau5(self, j, m):
"""
Calculate coefficient matrix element relative to (j+1, m, m).
"""
yD = self.dephasing
N = float(self.N)
num = (j - m) * (j + m) * (j + 1 + N/2)
den = 2 * j * (2 * j + 1)
t5 = yD * (num/den)
return t5
[docs] def tau6(self, j, m):
"""
Calculate coefficient matrix element relative to (j-1, m, m).
"""
yD = self.dephasing
N = float(self.N)
num = (j - m + 1) * (j + m + 1) * (N/2 - j)
den = 2 * (j + 1) * (2 * j + 1)
t6 = yD * (num/den)
return t6
[docs] def tau7(self, j, m):
"""
Calculate coefficient matrix element relative to (j+1, m-1, m-1).
"""
yP = self.pumping
N = float(self.N)
num = (j - m - 1) * (j - m) * (j + 1 + N/2)
den = 2 * j * (2 * j + 1)
t7 = yP * (float(num)/den)
return t7
[docs] def tau8(self, j, m):
"""
Calculate coefficient matrix element relative to (j, m-1, m-1).
"""
yP = self.pumping
yCP = self.collective_pumping
N = float(self.N)
num = (1 + N/2) * (j-m) * (j + m + 1)
den = 2 * j * (j+1)
pump = yP * (float(num)/den)
collective_pump = yCP * (j-m) * (j+m+1)
t8 = pump + collective_pump
return t8
[docs] def tau9(self, j, m):
"""
Calculate coefficient matrix element relative to (j-1, m-1, m-1).
"""
yP = self.pumping
N = float(self.N)
num = (j+m+1) * (j+m+2) * (N/2 - j)
den = 2 * (j+1) * (2*j + 1)
t9 = yP * (float(num)/den)
return t9