Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F63549328
QRSP_2modes_SchrodingerCatStates.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, May 20, 22:32
Size
3 KB
Mime Type
text/x-python
Expires
Wed, May 22, 22:32 (2 d)
Engine
blob
Format
Raw Data
Handle
17786784
Attached To
R10106 BosonicCodes
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
(
"
\n
spectral 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
)
Event Timeline
Log In to Comment