Homec4science
Diffusion BosonicCodes (master) Browse

R10106/g2_testcode.py
ff0cd64ecf07master

Last Change
Raw File

g2_testcode.py
View Options

import numpy as np
import scipy.sparse as sp
from scipy.integrate import complex_ode
# Number of Fock states and modes
from scipy.sparse import coo_matrix
N_max, N_modes = 8, 2 # choose even number of modes!
N_states = (N_max + 1) ** N_modes
"""initial parameters"""
U = 0.5 *1e-2
F0, theta_in, phi_in = 1e-1, np.pi/2, 0
pump = F0 * np.array([np.cos(theta_in/2), np.sin(theta_in/2) * np.exp(1j*phi_in)])
"""------------------"""
# 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, aa_dag = [], [] # initialize two empty lists for creation and annihilation operators for the corr. modes
for i in 1 + np.arange(N_modes):
a_tensor = sp.kron(sp.eye((N_max + 1) ** (N_modes - i)),
sp.kron(a, sp.eye((N_max + 1) ** (i - 1))))
aa.append(a_tensor)
aa_dag.append(a_tensor.transpose())
del a_tensor # free memory
nn = [] # initialize list for number of photon number occupation in each mode
for i in np.arange(N_modes):
nn.append(aa_dag[i].dot(aa[i]))
def hamiltonian(U, pump):
ham: coo_matrix = sp.coo_matrix((N_states, N_states), dtype=complex)
for j in range(N_modes):
ham += U/2 * aa_dag[j].dot(aa_dag[j]).dot(aa[j]).dot(aa[j]) \
+ pump[j] * aa_dag[j] + np.conj(pump[j]) * aa[j]
return ham
def liouvillian(rho, h, kappa):
l = -1j * (h.dot(rho) - rho.dot(h))
for i in range(N_modes):
l += kappa * (aa[i].dot(rho).dot(aa_dag[i]) - nn[i].dot(rho) - rho.dot(nn[i]))
return l
def super_liouvillian(h):
H = -1j * (sp.kron(h, sp.eye(N_states), ) - sp.kron(sp.eye(N_states), h.transpose()))
for i in range(N_modes):
H += sp.kron(aa[i], aa[i]) \
- 0.5 * sp.kron(sp.eye(N_states), nn[i]) - 0.5 * sp.kron(nn[i].transpose(), sp.eye(N_states))
return H
"""----------------"""
""" initialization """
"""----------------"""
# initial density matrix
rho = sp.lil_matrix((N_states, N_states), dtype=complex)
rho[0, 0] = 1
U = 1e-2
F0, theta_in, phi_in = 1e-1, np.pi/2, 0
pump = F0 * np.array([np.cos(theta_in/2), np.sin(theta_in/2) * np.exp(1j*phi_in)]) # np.random.uniform(-1e4, 1e4, N_modes)
h = hamiltonian(U, pump)
c4science · Help