all

Structural codes/Helper functions/PT_helper.py

+import numpy as np
+from qutip import *
+from tabulate import tabulate
+from scipy.sparse.linalg import splu
+
+
+
+
+def print_params(params):
+    Us  = params['Us']
+    Ur  = params['Ur']
+    G   = params['G']
+    η   = params['η']
+    g2  = params['g2']
+    χrs = params['χrs']
+    ϵd  = params['εd']
+    κr  = params['κr']
+    ξp  = params['ξp']
+    
+    floatfmt = ".2e"
+    
+    title   = dict(Dimensional_Vars='Value [MHz]')
+    #toprint = dict(Us=Us/(2*np.pi), Ur=Ur/(2*np.pi), G=G/(2*np.pi), χrs=χrs/(2*np.pi),\
+    #               η=η/(2*np.pi), κr=κr/(2*np.pi))
+    toprint = dict(Us=Us, Ur=Ur, G=G, χrs=χrs, η=η, κr=κr)
+    title.update(toprint)
+    print(tabulate(title.items(), headers='firstrow', floatfmt=floatfmt, tablefmt="fancy_grid", numalign="decimals"))
+    print("\n")
+    
+    
+    title   = dict(Dimensionless_Vars='Value [η]')
+    toprint = dict(Us=Us/η, G=G/η, χrs=χrs/η, κr = κr/η)
+    title.update(toprint)
+    print(tabulate(title.items(), headers='firstrow', floatfmt=floatfmt, tablefmt="fancy_grid", numalign="decimals"))
+    print("\n")
+    
+
+    title   = dict(Adiabaticity_conditions='Value [κr]')
+    toprint = dict(Us=Us/κr, g2=g2/κr, ϵd=ϵd/κr, χrs=χrs/κr)
+    title.update(toprint)
+    print(tabulate(title.items(), headers='firstrow', floatfmt=floatfmt, tablefmt="fancy_grid", numalign="decimals"))
+    print("\n")
+    
+    A = κr/(4*Ur*ξp**2)
+    B = κr/Us
+    ν = ϵd/κr
+    n = 2*ν*np.sqrt(B/(A+1/A))
+    
+    title   = dict(Independent_pars='Value')
+    toprint = dict(A=A, n=n, ν=ν, ξp=ξp, ξp2=ξp**2)
+    title.update(toprint)
+    print(tabulate(title.items(), headers='firstrow', floatfmt=".4f", tablefmt="fancy_grid", numalign="decimals"))
+    
+    return 
+
+
+
+
+def get_ss_fromblock(k, parity):
+    A = k.Lsymm[parity].data.tocsc()
+    m = A.shape[0]
+    msqrt = int(np.sqrt(m))
+    
+    b = np.ones(m)*1e-12
+    B = splu(A)
+    x = B.solve(b)
+    
+    ss = Qobj(x.reshape((msqrt,msqrt), order='F'))
+    ss = k.full_completion(operator_to_vector(ss),parity)
+    return ss
+
+
+def get_fidelity(ss, get_entropy=False):
+    N = ss.shape[0]
+    a = destroy(N)
+    parity = (1j*np.pi*a.dag()*a).expm()
+    
+    p = np.round(expect(parity, ss),2)
+    α = np.sqrt(expect(a**2,ss))
+    
+    if p==1.:
+        ss_est = ket2dm( (coherent(N, α) + coherent(N, -α)).unit() ) 
+    elif p==-1.:
+        ss_est = ket2dm( (coherent(N, α) - coherent(N, -α)).unit() )
+    else:
+        ss_est = 0
+        print("Error: ss is not an eigenstate of the parity operator.")
+        return 
+    f = ( (ss.sqrtm() * ss_est * ss.sqrtm()).sqrtm() ).tr()
+    if get_entropy : 
+        s = entropy_vn(ss)
+        return f.real, s
+    else:
+        return f.real
\ No newline at end of file

Structural codes/Kerr resonator/kerr_cavity.py

@@ -20,12 +20,18 @@ import time
 from matplotlib.animation import FuncAnimation
 from scipy.interpolate import interp1d
 from scipy.integrate import cumtrapz
-from DHS1 import dynamic_mesolve
+from itertools import chain
 
 from scipy.optimize import minimize_scalar
 from scipy.optimize import minimize
 
 
+
+def Fock_to_SymmFock(n):
+    I, J = range(n), list(chain(range(0, n, 2), range(1, n, 2)))
+    return spr.coo_matrix(([1]*n, (I, J)))
+
+
 class Parameters():
     def __init__(self, delta, U=0, G=0, F=0, gamma=0, eta=0, kappa_phi=0, get_scGc=False):
         self.η  = eta
@@ -138,10 +144,52 @@ class Kerr():
         rho = Qobj(vec, dims=[[[self.msqrt],[self.msqrt]],[1]])
         rho = self.complete_vec(rho,symm,invert=True)
         rho = vector_to_operator(rho)
-        rho+=rho.dag()
-        rho/=rho.tr()
+        #rho+=rho.dag()
+        if rho.tr()!=0.:
+            rho/=rho.tr()
         return rho
-
+    
+    
+    def get_ss(self, symm_FS=False, get_left=False):
+        if symm_FS : P_SF = Fock_to_SymmFock(self.N)
+        if get_left: Jlist = []
+        
+        sslist = []
+        for i in range(int(len(self.blocks_list))):
+            block = self.blocks_list[i]
+            _, evecs = Qobj(block).eigenstates()    
+            
+            ss_block_form = evecs[-1]
+            
+            eigvs = [ss_block_form]
+            if get_left:
+                R = spr.hstack(evecs).tocsc()
+                L = spr.linalg.inv(R)
+                eigvs.append(Qobj(L[-1].H))
+            
+            for j in range(len(eigvs)):
+                vv = spr.dok_matrix((self.LL.shape[0], 1), dtype='complex')
+                vv[self.bl_indices[i]:self.bl_indices[i]+self.block_sizes[i]] = eigvs[j].data
+                vv = vv.tocsr()
+                vv = self.P.data.transpose().dot( vv )
+                vv = vv.reshape((self.N,self.N), order='F')
+                eigvs[j] = vv
+            
+            #symmetric FS
+            if symm_FS: 
+                for v in eigvs:
+                    v = P_SF.dot(v.dot(P_SF.H))
+                
+            ss = Qobj(eigvs[0]) 
+            ss/= (ss.dag()*ss).tr()
+            sslist.append(ss)
+            
+            if get_left: Jlist.append(Qobj(eigvs[1]))
+                
+        if get_left:return sslist, Jlist
+        else: return sslist
+        
+    
     def get_spectrum(self,ret_ss=False):
         done=False
         for i in range(int(len(self.blocks_list))):