Calculating Kraus Operators

[Erin Flk]

Dear all,

I want to calculate the Kraus Operators, from my unitary matrix which I get from sesolve. The Kraus Operators don't seem to fulfill the identity-condition. How can I correct this code? 

The Code I present here also includes the calculation of the unitary matrix for completeness:

from scipy import *

from numpy import *

import scipy.constants as c

import numpy as np

import matplotlib as mpl

import matplotlib.pyplot as plt

from scipy. linalg import expm, sinm, cosm, eigh

import math

%matplotlib inline

import qutip as qt

from numpy.random import RandomState

from scipy.interpolate import CubicSpline

from scipy import integrate

from numpy.linalg import inv

import time

from qutip import tensor, basis, qeye, destroy, sigmax, sigmam,sigmay, sigmaz, thermal_dm, coherent_dm, fock_dm, mesolve,mcsolve , ket2dm, Options

import sys

import tqdm

from scipy.signal import find_peaks

import sympy as sp

from qutip import *

delta = 3.71*1e3*2*np.pi#800 * 1e3 * 2 * np.pi #motional detuning in Hz

Omega = 1.071*1e3*2*np.pi#5.2 * 1e6 * 2 * np.pi * eta

t_max = delta * 5e-2 # = 1200 ca

num_points = 1000 # Number of points for time evolution

Omega_tilde = Omega / delta

# H = S_x * (a*e^(i delta t) + a^(dagger)*e^(-i delta t)) -> not Magnus expanded

def ExpPos(t, args):

#delta = args['delta']

return np.exp(1j * t)

def ExpNeg(t, args):

#delta = args['delta']

return np.exp(-1j * t)

def compute_time_evolution_ideal(Omega, t_max, num_points):

# Define the dimension of the motional mode Hilbert space

dim_motional = 2 # Truncation level for the motional mode

dim_qubit = 2 # Qubit has 2 states

# Define annihilation and creation operators for the motional mode

a = tensor(qeye(dim_qubit), qeye(dim_qubit), destroy(dim_motional))

adag = tensor(qeye(dim_qubit), qeye(dim_qubit), create(dim_motional))

# Define collective spin operator along x for the qubits

sx1 = tensor(sigmax(), qeye(dim_qubit)) # σ_x on qubit 1

sx2 = tensor(qeye(dim_qubit), sigmax()) # σ_x on qubit 2

# Define collective angular momentum operator along x

S_x = (sx1 + sx2) / 2

S_x_nb = tensor(S_x, qeye(dim_motional))

# Hamiltonian H = Omega * S_x * (ae^(i delta t) + a^(dagger)e^(-i delta t))

H1 = Omega * S_x_nb * a

H2 = Omega * S_x_nb * adag

H = [[H1, ExpPos], [H2, ExpNeg]]

# Time evolution using the sesolve solver for unitary evolution

tlist_ideal = np.linspace(0, t_max, num_points) # Time list from 0 to t_max with num_points

# Use the identity as the initial state to calculate the unitary operator

identity_full = qeye(dim_qubit * dim_qubit * dim_motional)

# define args

#args = {'delta': delta}

# Solve the Schrödinger equation to get the unitary evolution operator

result = sesolve(H, identity_full, tlist_ideal, options = Options(nsteps=1000000))

# Convert to reduced unitary operator acting on the qubit subspace

U_t_ideal = [Qobj(state.full(), dims=[[dim_qubit, dim_qubit, dim_motional], [dim_qubit, dim_qubit, dim_motional]]) for state in result.states]

return U_t_ideal, tlist_ideal

# Compute the time evolution once

U_t_ideal, tlist_ideal = compute_time_evolution_ideal(Omega_tilde, t_max, num_points)

phi_gate = np.pi / 2

alpha = np.exp(-1j * phi_gate) + 1

beta = np.exp(-1j * phi_gate) - 1

spin_unitary = Qobj(0.5 * np.array([[alpha, 0, 0, beta],

[0, alpha, beta, 0],

[0, beta, alpha, 0],

[beta, 0, 0, alpha]]),

dims=[[2, 2], [2, 2]])

motion_unitary = qeye(2)

target_unitary = tensor(spin_unitary, motion_unitary)

print(spin_unitary) # my target unitary that I am using for the calculation of the fidelity

# Define parameters

dim_qubit = 2

dim_motional = 2

n_max = dim_motional - 1 # Maximum environment state index

nbar = 0.5 # Average occupation number for thermal distribution

# Function to compute the thermal occupation probability mu_n

def mu_n(n, nbar):

return np.sqrt((nbar**n) / ((1 + nbar)**(n + 1)))

def compute_kraus_operators_for_all_times(U_t_kraus, tlist_kraus, n_max, nbar):

kraus_operators_all_times = []

for idx, t in enumerate(tlist_kraus):

U_sm = U_t_kraus[idx] # Unitary operator at this time

kraus_operators = []

for n_dash in range(n_max + 1): # Iterate over final motional states

# Define the initial and final motional states

n = 0

initial_motional_state = basis(dim_motional, n)

final_motional_state = basis(dim_motional, n_dash)

# Construct the outer product of the initial and final motional states

vector = initial_motional_state * final_motional_state.dag()

# Construct the tensor product with the identity operator for the qubit subsystem

vector2 = tensor(qeye(2), qeye(2), vector)

kraus_operator = vector2 * U_sm * mu_n(n, nbar)

# Perform the partial trace over the motional subsystem

kraus_operator = kraus_operator.ptrace([0, 1])

kraus_operators.append(Qobj(kraus_operator))

K_multiplied_summed = 0

for i in range(2):

K_multiplied = kraus_operators[i].dag()*kraus_operators[i]

K_multiplied_summed += K_multiplied

#print(K_multiplied_summed)

kraus_operators_all_times.append(kraus_operators)

return kraus_operators_all_times

def get_kraus_operators_at_time_t(kraus_operators_all_times, tlist_kraus, t):

idx = (np.abs(tlist_kraus - t)).argmin()

print("idx:", idx, "t =", t)

return kraus_operators_all_times[idx]

kraus_operators_all_times = compute_kraus_operators_for_all_times(U_t_ideal, tlist_ideal, n_max, nbar)

t = 2 * np.pi #/ delta

kraus_operators_at_t = get_kraus_operators_at_time_t(kraus_operators_all_times, tlist_ideal, t)

#print(kraus_operators_at_t)

super_kraus_at_t = kraus_to_super(kraus_operators_at_t)

avg_fidelity_at_t = average_gate_fidelity(super_kraus_at_t, spin_unitary)

print("Average Gate Fidelity at t:", avg_fidelity_at_t)

[Paul Menczel]

Hi, 

If I may give some general advice: you posted a lot of code and very little explanation. More explanation and code trimmed down to what is necessary will make it more likely that people engage with your question. 

A comment that I can give without looking at the code in detail is that you might have some kind of mathematical misunderstanding. You say you calculate your evolution with sesolve, meaning it is a unitary evolution. However, for unitary evolution, there is nothing to do to calculate Kraus operators: there is only one Kraus operator, which is the unitary evolution operator. 

I also note that qutip has a function to calculate Kraus operators automatically, call "to_super(qobj)" on a superoperator Qobj: https://qutip.readthedocs.io/en/stable/guide/guide-states.html#choi-kraus-stinespring-and-chi-representations

Cheers

[Paul Menczel]

Hi, it seems that you accidentally sent your reply to my personal email address instead of the mailing list. For the sake of other readers, I will copy the reply here:

> Hi Paul,
> thank you for your advice. I do understand that a unitary operator can be a Kraus-Operator. But as my unitary operator is of the joined system spin x motion and I want to calculate E(rho) (quantum operation) to use for my fidelity calculation, shouldn’t I need to calculate „new“ Kraus-Operators that are on the subsystem of the spin, as I take the partial trace over the motional part?
> Best regards

Yes, that is correct, you have to calculate new Kraus operators. In your original question, you did not mention that you have a joint (spin - motion) system, so I did not know that.
Beyond that, I am sorry, but I do not have any advice other than what I wrote above. Please remember that we are volunteers answering questions in our free time. If there is only little effort into explaining to us your problem and your code, it gives us little motivation for spending a long time trying to understand your code and find the mistake (which might not even be related to qutip at all).