Homec4science
Diffusion BosonicCodes (master) Browse

R10106/QRSP_2modes_SchrodingerCatStates.py
22796502c783master

Last Change
Raw File

QRSP_2modes_SchrodingerCatStates.py
View Options

import numpy as np
import scipy.sparse as sp
from math import floor
from scipy.integrate import complex_ode
# Number of Fock states and modes
from scipy.sparse import coo_matrix
N_max, N_modes = 9, 2 # choose even number of modes!
N_states = (N_max + 1) ** N_modes
# Define annihilation and creation operator for a single mode
a = sp.diags(np.sqrt(1 + np.arange(N_max)), offsets=1)
a_dag = a.transpose()
# Define annihilation and creation operator for entire Hilbert space
aa_list, aa_dag_list = [], [] # initialize two empty lists for creation and annihilation operators for the corr. modes
for i in 1 + np.arange(N_modes):
aa = sp.kron(sp.eye((N_max + 1) ** (N_modes - i)),
sp.kron(a, sp.eye((N_max + 1) ** (i - 1))))
aa_list.append(aa)
aa_dag_list.append(aa.transpose())
del aa # free memory
nn_list = [] # initialize list for number of photon number occupation in each mode
for i in np.arange(N_modes):
nn_list.append(aa_dag_list[i].dot(aa_list[i]))
# Dimensions of the reservoir lattice
Lx = floor(np.sqrt(N_modes))
Ly = floor(N_modes / Lx)
def spectral_rad(lx, ly, onsite, hop):
mat = np.diag(onsite)
ctr = 0
for ix in np.arange(lx):
for iy in np.arange(ly):
# print(ix, iy)
i = ly * ix + iy # site-index in the operator list
# define neighbors:
for jx, jy in [(0, 1), (1, 0)]:
ii = ly * (ix + jx) + iy + jy
if 0 <= ix + jx < Lx and 0 <= iy + jy < Ly: # make sure not to leave the boundaries
mat[i, ii] = hop[ctr]
mat[ii, i] = hop[ctr] # make sure the connection is bidirectional (symmetric)
ctr += 1
print("coupling matrix:")
print(mat)
rad = max(abs(np.linalg.eigvals(mat) ** 2))
print("\nspectral radius:", rad)
return rad
def hamiltonian(lx, ly, onsite, hop, alpha, pump):
ham: coo_matrix = sp.coo_matrix((N_states, N_states), dtype=complex)
for j in range(N_modes):
# 1. onsite-energies 2. nonlinearities 3. pumping
ham += onsite[j] * nn_list[j] # onsite-energies
# Loop over site indices
ctr = 0 # set counter for hopping array
for ix in np.arange(lx):
for iy in np.arange(ly):
i = ly * ix + iy # site-index in the operator list
# define neighbors:
for jx, jy in [(0, 1), (1, 0)]:
ii = ly * (ix + jx) + iy + jy
if 0 <= ix + jx < Lx and 0 <= iy + jy < Ly:
print(i, ii, ctr)
ham += hop[ctr] * (aa_dag_list[i].dot(aa_list[ii]) + aa_dag_list[ii].dot(aa_list[i])) # Hopping
print("done")
ctr += 1
# compute spectral radius of the onsite + hopping part of the Hamiltonian
ham.multiply(1 / spectral_rad(lx, ly, onsite, hop)) # rescale onsite + hopping term by spectral radius
for j in range(N_modes):
ham += alpha[j] * aa_dag_list[j].dot(aa_dag_list[j]).dot(aa_list[j]).dot(aa_list[j]) \
+ pump[j] * aa_dag_list[j] + np.conj(pump[j]) * aa_list[j]
return ham
def liouvillian(rho, h, gam):
l = -1j * (h.dot(rho) - rho.dot(h))
for i in range(N_modes):
l += gam * (aa_list[i].dot(rho).dot(aa_dag_list[i]) - nn_list[i].dot(rho) - rho.dot(nn_list[i]))
return l.reshape((N_states**2))
"""----------------"""
""" initialization """
"""----------------"""
# initial density matrix
rho = sp.lil_matrix((1, N_states**2), dtype=complex)
rho[0, 0] = 1
# random initialization parameters
onsite = np.random.random(N_modes)
hop = np.random.random((Lx - 1) * Ly + (Ly - 1) * Lx)
alpha = np.random.random(N_modes)
pump = np.random.random(N_modes)
gamma = 1
h = hamiltonian(Lx, Ly, onsite, hop, alpha, pump)
ode = complex_ode(liouvillian)
ode.set_f_params(h, gamma)
ode.set_integrator("dopri5")
ode.set_initial_value(rho)
c4science · Help