國立臺灣大學應用物理研究所 碩士論文

## Graduate Institute of Applied Physics National Taiwan University

## Master Thesis

非馬可夫開放系統中含流失態的 CNOT 量子邏輯閘最 佳化控制

Optimal Control of CNOT Gate Operation with Leakage states and Non-Markovian Environments

楊存毅 Tsun-Yi Yang

指導教授：管希聖 博士 Advisor: Hsi-Sheng Goan, Ph.D.

中華民國 102 年 7 月

口口口試試試審審審定定定書書書

誌誌誌謝謝謝

碩士論文能夠如期完成，首先要感謝我的老師管希聖教授這麼詳細的幫我改論文還 有指導我，真的比我預想的要詳細很多，老師真的很用心，還麻煩了胡崇德教授和 蘇正耀教授當口試委員，很感謝兩位百忙之中還抽空前來，還有黃上瑜學長不厭其 煩的指導我，還常常指出我一些錯誤的觀念，沒有學長我大概連程式也很難做出 來，另外還要感謝我的良師益友戴榮身常常與我討論問題，感謝好戰友林冠廷在我 脆弱的時候還有很多地方幫助我，感謝跟周宜之前合作讓我學到很多，還有感謝簡 崇欽學長、陳柏文學長、張晏瑞學長、黃斌學長、萬哥、洪常力、李彥賢，還有許 多學長學弟們在501室的交流和合作都讓我感動和感恩。

感謝我的爸爸和媽媽和弟弟，感謝我在西方極樂世界的阿嬤和阿公，還有所有 的親戚朋友，雖然爸爸在我讀大學和碩士的這段期間受到了身體上不小的打擊，但 和媽媽總是兩人一起很堅強，日後是我面對一切困難的榜樣，這是比其他一切都更 重要的課程。 有太多人要感謝，但是一時之間我也無法一一數盡，碩士這兩年我 看了很多人、很多事，遇過很多好人，也看過我無法理解的人，在專業素養的養成 以外，我覺得做研究的心態才是重要的，如果只是為了學歷或錢去做研究，那肯定 是做不好的，我想以後我會盡力的調整自己的心態，去匹配像老師和學長們那樣的 傑出研究者，在生活上，但求無愧於心，了然於胸的是自己的心願，如果以後一生 我也能堂堂正正的面對自己的一切抉擇，我想我也無遺憾了，而碩士這兩年，確確 實實的幫助我了解了更進一步的人生觀，期望自己以後能遇到共度一生的人，分享 我此時此刻的心願，大概也就值得了吧。

另外還要感謝我的室友和好朋友們，陳凱評、馬伯倫、施景堯、沙政亞，都是 我的好朋友和人生導師，因為遇見你們讓我變得更加成熟，也受過你們的幫助，讓 我覺得認識你們真好。

再一次的，感謝所有人。

中中中文文文摘摘摘要要要

為了建造一個量子電腦，高精準度的邏輯閘運算是必要的。特別是兩個量子位 元(two qubit)所能進行的控制非門(CNOT)的運算，是一個很重要的量子邏輯閘。

然而，現實世界中有需多的問題導致實驗上的結果和理論上的預測是不盡相同的。

在這篇碩論中，我們探索了在量子最佳化控制理論下，考慮了流失態(leakage)和非 馬可夫(non-Markovian)環境下的控制非門的運算。首先，我們給了一個簡短的超 導體量子位元介紹，然後決定了單量子位元的漢米爾頓(Hamiltonian)矩陣。然後 我們描述了環境的雜訊能譜(noise power spectrum)，接著就用費米黃金定律(Fermi golden rule)去連結馳耗速率(relaxation rate)和關聯函數(correlation function)之間 的關係。然後我們描述了一個非馬可夫環境下的控制方程式，接著用我們的模型實 行了量子最佳化控制。科羅多夫最佳化控制是我們在這篇碩論中採用的方法，目的 是為了最小化邏輯閘運算的誤差，然後找到相對應的控制參數，並且考慮非同時的 記憶效應(non-local memory effect)的非馬可夫開放系統。我們也討論了不同控制參 數的不同波形，以及控制參數和關聯函數影響最佳化邏輯閘控制的行為。我們發現 在環境的參數影響下(此參數參考真實實驗系統結果)，使用超導體電量量子位元是

可以完成一個誤差約為10^{−4} ∼ 10^{−5}的控制非門的高精準度邏輯閘控制的。

中中中文文文關關關鍵鍵鍵字字字：：：最最最佳佳佳化化化、、、位位位元元元、、、邏邏邏輯輯輯閘閘閘、、、科科科羅羅羅多多多夫夫夫、、、約約約瑟瑟瑟芬芬芬

Abstract

When building a quantum computer, high precision gate operations are needed. In particular, the controlled-not (CNOT) operation regarded as a crucial universal two- qubit gate is a very important quantum gate to implement. However, real world contains a lot of problems and causes the difference between experimental results and theoretical simulations. In this thesis, we investigate CNOT gate operation using quantum optimal control theory for superconducting charge qubit system taking into account the effects of leakage states and a non-Markovian environment. First, we give a brief introduction to superconducting qubits and decide the Hamiltonian of a simple one-qubit model. Then we describe the noise power spectrum of environments, which gives the relation between the relaxation rate and the bath correlation function through the Fermi golden rule . After that, we describe the non-Markovian master equation approach and apply it to our model together with the quantum optimal control theory. The Krotov optimal control method that we used in this thesis can minimize the error of gate operations and find the corresponding optimal control pulses to realize the gate operations. Considering the non-local memory effect in non-Markovian open quantum systems. We also discuss the effect of different shapes and behaviors in the bath correlation function on the optimal control gate fidelity.

We find that it is possible to implement high-fidelity CNOT gates with error about
10^{−4} ∼ 10^{−5} in superconducting charge qubit system with environment parameters
extracted from the realistic noise power spectrum of experiments.

Keywords: optimal, qubit, gate, Krotov, Josephson

Contents

口口口試試試審審審定定定書書書 i

誌誌誌謝謝謝 ii

中中中文文文摘摘摘要要要 iii

Abstract iv

1 Introduction and Background 1

1.1 Introduction . . . 1

1.2 Review of Josephson Charge Qubit . . . 2

1.3 Hamiltonian of the One Josephson Charge Qubit Model . . . 2

1.4 Power Spectrum and Noise Fluctuation . . . 5

1.5 Transition Rate . . . 8

2 Superconducting Qubit and Noise 11 2.1 Overview . . . 11

2.2 Hamiltonian of Two Coupled Josephson Charge Qubits . . . 12

2.3 Spectral Density and Correlation Function . . . 18

2.4 Diagram of simulation . . . 22

3 Optimal Control 25 3.1 Introduction . . . 25

3.2 Error/Fidelity Definition . . . 26

3.4 Two Qubit Optimal Control . . . 31

3.5 State-independent Optimal Control . . . 31

3.5.1 Computational States and Leakage States . . . 31

3.5.2 Superoperator and Column Representation . . . 37

4 Results 41 4.1 Overview . . . 41

4.2 Closed System with No Leakage Levels . . . 43

4.2.1 Overview . . . 43

4.2.2 Pulse Shape . . . 43

4.3 Open System with No Leakage Levels . . . 44

4.3.1 Overview . . . 44

4.4 Closed System with Leakage Levels . . . 51

4.4.1 Overview . . . 51

4.4.2 Pulse Shape . . . 54

4.5 Open System with Leakage Levels . . . 55

5 Conclusion 59 A RK4 and expm 61 B Derivatives of Matrix, Traces 65 B.1 Matrix Multiplication . . . 65

B.2 Derivatives of Matrix . . . 66

B.3 Derivatives of Trace . . . 66

B.3.1 Derivation of Trace REAL Matrix . . . 66

B.3.1.1 First Order . . . 67

B.3.1.2 Second Order . . . 67

B.3.2 Derivation of Trace COMPLEX Matrix . . . 67

B.3.2.1 Generalized Complex Derivative: . . . 68

B.3.2.2 Useful Formula . . . 68

C Leakage Optimal Calculation 69

Chapter 1

Introduction and Background

## 1.1 Introduction

The field of solid-state quantum computation grows very fast. Since some quantum algorithms (e.g. Shor’s algorithm) are shown to outperform significantly their best known classical counterparts, the idea and implementation of quantum information processing becomes important. The research of quantum bits, or qubits, is rapidly growing. There are many types of physical objects that have potential to be imple- mented as qubits. However, solid-state circuits, and superconducting circuit are of great interest as they offer great scalability. That is, the possibility of making circuits with a larger number of interacting qubits. The Josephson charge qubit circuit is a good candidate because it can implement qubit, i.e., quantum two-level systems.

Such systems possess coherence, can be initialized and read out. Also, the logic gates operations can be constructed and gate operation can be done. First, In chapter 1, we introduce the background such as Josephson charge qubit Hamiltonian, power spec- trum, and transition rate. In chapter 2, we describe our main purpose, i.e., two-qubit gate operations and present some necessary material that are useful to achieve our goal. We will give an overview of the performing two qubit Hamiltonian. We will dis- cuss the noise spectral density and the bath correlation functions. We will try to use a more realistic bath spectral density parameters extracted from relevant experiments.

In chapter 3, we introduce the Krotov optimal control method and the derivation to

obtain the iteration equations. In chapter 4, we present the optimal control results of four different kinds of situations. They are the combination of open/closed systems and with/without leakage levels. Finally, we conclude this whole thesis in Chapter 5.

Some details about how to perform the calculation are put in appendices.

## 1.2 Review of Josephson Charge Qubit

The Josephson charge qubit consists basically of a Cooper pair box as a supercon-
ducting island which is connected via a Josephson junction to a large electrode, called
a superconducting reservoir. Figure 1.1 shows a schematic illustration of a Joseph-
son charge qubit. There are some paper reviewing the property of superconducting
qubits, [1, 16]. Typically there are around 10^{7} − 10^{8} conduction electrons in the is-
land and the dimensions are about 1000 nm × 50 nm × 20 nm. The Cooper pairs can
tunnel into or out of the box with a tunneling amplitude described by the Josephson
coupling E_{J} between the box and the reservoir. The number of the Cooper pairs on
the island can be tuned by the voltage V_{g} applied to the gate coupled to the box
through the gate capacitor C_{g}.

## 1.3 Hamiltonian of the One Josephson Charge Qubit Model

The Hamiltonian of a Josephson charge qubit can be written as

H = 4Ec(ˆn − ng)^{2}− Ejcos ˆΘ, (1.1)

where n is the number operator of (excess) Cooper-pair charges on the island, and
Θ, the phase of the superconducting order parameter of the island (phase difference
across a Josephson junction), is its quantum-mechanical conjugate, n = −i~_{∂(~Θ)}^{∂} , Ej

Figure 1.1: Schematic illustration of a Josephson charge qubit; from [1]. The bold

“n” is the number of Cooper pairs that tunnel through the insulator into the island or

“box”. C_{J} is the junction capacitance, C_{g} is the gate capacitance. E_{J} is the junction
energy and Vg is the gate voltage.

is the Josephson coupling energy.

E_{C} = e^{2}

2(C_{g} + C_{J}) (1.2)

is the single electron charging energy and depends on the total capacitance of the island. The dimensionless gate charge

n_{g} = C_{g}V_{g}

2e (1.3)

accounts for the effect of the gate voltage. We will vary the gate voltage as our external control.

Since we will control the tunneling Cooper pairs, we choose the number states describing the number of the Cooper pairs on the island as our basis, the so called the charge basis. In order to express the Hamiltonian by the basis, we introduce the relation between n and Θ:

[ ˆΘ, ˆn] = i (1.4)

with ˆΘ ∈ [0, 2π]. One can use Eq.(1.4) to derive [e^{i ˆ}^{Θ}, ˆn] = −e^{i ˆ}^{Θ} and [e^{−i ˆ}^{Θ}, ˆn] =
e^{−i ˆ}^{Θ}. This is similar to the commutation relation between the creation operator, the
annihilation operator and the number operator of a harmonic oscillator. So one can
regard e^{±i ˆ}^{Θ}as creation and annihilation operators, respectively. Inserting an identity
I =P

n

|nihn| onto the both sides of the Hamiltonian, Eq.(1.1), we obtain

H = X

n

{4E_{c}(n − n_{g})^{2}|nihn| − 1

2E_{j}(|n + 1ihn| + |nihn + 1|)}. (1.5)
One may have seen different form of the Hamiltonian from the form of Eq.(1.3), but
they are actually the same just because sometimes people would like to shift identities
to make the equation look more simple than it’s original form.

Usually, one considers the charge states corresponding to n = 0, 1 as the com- putational states of the qubit to perform a gate operation. However, as in [2], we know that the dynamics of the Josephson charge qubit system is affected by the other number state, which is the leaky qubit effect, or so called leakage. That is, the con- sideration including other states into the calculations may have large difference to the case by considering only the computational states {|0i, |1i}. To study the leakage effect, in this thesis, we consider the charge states {| − 1i, |0i, |1i, |2i}. Many related papers point out that there is no need to consider even more charge states.[2, 3, 4].

One will see that the control problem is not physical enough if one does not consider the leakage effect. One may ask what n = −1 state is. What is minus number of Cooper pairs? We have to see this in a bigger picture. Since the number of Cooper pairs on the island should be a huge amount instead of just one or two, the number of Cooper pairs on the island we discussed before is the change in the number of the Cooper pairs as compared with the equilibrium number, i.e., the number of the excess Cooper pairs.

Note that some people might use energy eigenstates like Figure 1.2 to construct their two-level system. Here we just use the charge basis which are eigenstates of the number operator with eigenvalue describing the number of Cooper-pair on the island.

Figure 1.2: The charge energy of the Cooper-pair box as a function of n_{g}. From [1].

Since n_{g} is not an actual number of charge, it can be a float number. When n_{g} = ^{1}_{2},
the system with states corresponding to the energy splitting caused by E_{j} can be
approximate as a two-level system. The dashed lines means the changing energy
with definite Cooper pairs n on the island, and the solid line stands for the energy
eigenstates.

## 1.4 Power Spectrum and Noise Fluctuation

There must exist some noise in every quantum system to prevent us from perfectly controlling the quantum system. For Josephson charge qubit, there are some fluctu- ation of the voltage that will cause noise.

Let’s rewrite Eq.(1.1) as in [7]

H_{box}= 4E_{C}(n − Q_{g}

2e)^{2}− E_{j}cos Θ, (1.6)

where Q_{g} = C_{g}V_{g} is the gate charge which is controlled by gate voltage and E_{C} =

e^{2}
2(Cj+Cg).

By properly choosing the right parameters such as voltage, temperature, fre- quency, and the strength of control pulses, only two charge states |n = 0iand |n = 1i play an important role. In the spin notation, the number operator becomes n =

1

2(1 − σ_{z}), while cos Θ = ^{1}_{2}σ_{x}, and the effective spin Hamiltonian after inserting
V_{g}+ δV , where δV denotes the voltage fluctuation, into Eq. (1.6) reads

H = −1

2B_{z}(V_{g})σ_{z}−1

2B_{x}σ_{x}− 1

2Xσ_{z}+ H_{bath}, (1.7)

where the parameters are B_{z} = _{C} ^{2e}_{+C} (e − C_{g}V_{g}), B_{x} = E_{j}, and B_{z}(V_{g} + δV ) =

Figure 1.3: From [7]. (a) Circuit of a Josephson charge qubit. (b) Circuit representing noise coming from fluctuation denoted by δV .

B_{z}(V_{g}) + X. The sum of the voltage V_{g} + δV couples to the charge in the circuit.

The coupling term involving X, which comes from the fluctuation of voltage rep- resented as the interaction with a bath results in decoherence and dissipation. Here

X = −(2eC_{t}

C_{j} )δV. (1.8)

Since we assume interaction with the bath or environment is sufficient weak, [19, 20, 21], we can represent the environment as a set of harmonic oscillators with a coupling linear in the oscillator coordinates and/or momenta to the system; by appropriate canonical transformation and related tricks, it is possible to ensure that the “coupling is only to the coordinates”. That is

H_{bath} =X

i

1 2mi

p^{2}_{i} +m_{i}ω_{i}
2 x^{2}_{i}

(1.9)

and the coupling part of the environment is

δV =X

i

λ_{i}x_{i}, (1.10)

where λ_{i} is the coupling coefficient with each coordinate x_{i} of the oscillators, where

xi = 1

√2mω_{i}(ai+ a^{†}_{i}) (1.11)

We will see later in Sec. 2.2 that X represents the coupling part of δn there:

X =X

i

¯λ_{i}x_{i} =X

i

¯λ_{i}(a_{i}+ a^{†}_{i}), (1.12)

where ¯λ_{i} absorbs the rest of the parameters of Eq. (1.8).

From the Hamiltonian of eigenbasis, one can easily see the energy splitting of the system. Although one may not always choose eigenbasis, this representation can let one see the relation between the power spectrum and the relaxation rate from the Fermi golden rule. The Hamiltonian, Eq.(1.7), written in the qubit energy eigenbasis becomes

H = −1

2 4 Eσz− 1

2X [cos(η)σz+ sin(η)σx] + Hbath (1.13)
with the energy difference 4E = pB_{x}^{2}+ B_{z}^{2}, tan η = ^{B}_{B}^{x}

z, X_{tran} = −^{1}_{2}X sin η, and
V (t) = X_{tran}(t). Why do we want to connect the power spectrum with the relaxation
rate? The reason is that we normally measure the relaxation rate to define how much
the environment affects our system, and the noise power spectrum and the bath
spectral density are related to the bath correlation function. The following equations
tell us how the relaxation rate is related to the noise power spectrum.

First we know from the Fermi golen rule that the transition rate between two different states |ai and |bi is

Γ_{ba} = 1

~^{2}
ˆ ∞

−∞

dthV_{ab}(t)V_{ba}(0)i_{B}e^{−iω}^{ba}^{t}, (1.14)

where V denotes the interaction Hamiltonian which is assumed to be weak. This equation will be discuss in the next section. And we know there must be two directions

of the transition

Γ↓ = Γ_{10}= sin^{2}η

4~^{2} hX^{2}_{ω=}4E

~

i, (1.15)

Γ↑ = Γ_{01}= sin^{2}η

4~^{2} hX^{2}_{ω=−}4E

~

i, (1.16)

where hX^{2}_{ω}i = ´

dte^{iωt}hX(t)X(0)i. Combining these two rates, we obtain the relax-
ation rate

Γ_{1} = Γ↓+ Γ↑ = π sin^{2}η

2~^{2} S_{U}(ω = 4E

~ ) (1.17)

with the noise power spectrum S_{U}(ω) = _{2π}^{1} (hX^{2}_{ω}i + hX^{2}_{−ω}i). The relation between the
noise power spectrum and the bath spectral density will be discussed in Sec. 2.3 .
We can describe the evolution of the density matrix of the qubit via the equation of
motion:

T r ˆρ = ρ_{00}+ ρ_{11} = 1, (1.18)

˙

ρ_{00} = − ˙ρ_{11} = −Γ↑ρ_{00}+ Γ↓ρ_{11,} (1.19)

## 1.5 Transition Rate

Since different papers have their own definition of the noise spectrum, power spec- trum, spectral density, with or without 2π of the Fourier transformation, we can check their definition of the transition rate, and relate it to the noise terms men- tioned above. Then we can compute the environment effect using the noise power spectrum obtained from the experimental data.

Return to how we get the transition rate from the Fermi golden rule [18]. Consider a simple model

H = (H_{S}+ H_{B}) + V = H_{0}+ V, (1.20)

where V is the interaction Hamiltonian describing weak coupling the between system and the bath.

Assume that the system has two eigenstates |ai and |bi and their corresponding
eigenenergies are Ea and Eb, respectively. Similarly, we assume that the bath mode
γ has eigenstates |γi and eigenenergy E_{γ}. The Hamiltonian H_{S} and H_{B} written in
their eigenbasis become

HS = Ea|aiha| + Eb|bihb|, (1.21)

H_{B} =X

γ

E_{γ}|γihγ| (1.22)

with γ = · · · α, β · · · . Thus, we have

H_{0}|aαi = (E_{a}+ E_{α})|aαi, (1.23)

where |aαi = |ai ⊗ |αi is the product state of the system and the bath. Set |aαi the initial state and |bβi the final state. Consider the case where the system energy decreases and the bath energy rises.

Γ_{f i} = 2π

~ X

i,f

p_{f}|hi|V|f i|^{2}δ(E_{f} − E_{i})

= 2π

~ X

a,b,α,β

p_{bβ}|haα|V|bβi|^{2}δ((E_{b}+ E_{β}) − (E_{a}+ E_{α}))

= 1

~^{2}
ˆ ∞

−∞

dt X

a,b,α,β

p_{bβ}hbβ|V|aαihaα|V|bβie^{−}^{~}^{i}^{[(E}^{b}^{−E}^{a}^{)+(E}^{β}^{−E}^{α}^{)]t}, (1.24)

where p_{f} is the probability of the final state |f i. Therefore, the transition rate system
from the system eigenstate |ai to |bi is

Γ_{ba} = 1

~^{2}
ˆ ∞

−∞

dtX

α,β

p_{β}hβ|e^{~}^{i}^{E}^{α}^{t}V_{ba}e^{−}^{~}^{i}^{E}^{β}^{t}|αihα|V_{ba}|βie^{−iω}^{ba}^{t} (1.25)

Define

V_{ab} = ha|V|bi. (1.26)

In interaction picture, one can get

V_{ab}(t) = e^{~}^{i}^{H}^{B}^{t}V_{ab}e^{−}^{~}^{i}^{H}^{B}^{t}. (1.27)

The bath correlation function is defined as

C_{ba}(t) = hV_{ab}(t)V_{ba}(0)i_{B}, (1.28)

where

h· · · iB =X

β

hβ| · · · |βi (1.29)

means an equilibrium thermal average over the bath states. Then, we have

Γ_{ba} = 1

~^{2}
ˆ ∞

−∞

dthV_{ab}(t)V_{ba}(0)i_{B}e^{−iω}^{ba}^{t}. (1.30)

So we can use this form to the relation between the relaxation rate and the bath correlation function.

Chapter 2

Superconducting Qubit and Noise

## 2.1 Overview

The basic element of quantum information processing is qubit. There are many physical systems that can implement qubits. Here we choose the Josephson super- conducting charge qubit system to study. Like other qubits, the Josephson charge qubit has its own weakness; especially it is susceptible to charge or voltage noise and has the problem of qubit state leakage. However it is still a good physical system for quantum computation, if we can perform quantum gate operation with high fidelity taking into account the influence of noise and environment. So the goal of this thesis is to find control field pulse to achieve high-fidelity controlled-not (CNOT) gate oper- ation by the tool called the Krotov optimal control method, which will be introduced in Chapter 3. In this Chapter, we first present the Hamiltonian of two Josephson charge qubits. Then we describe the decoherence effect from the environment and the leakage states.

## 2.2 Hamiltonian of Two Coupled Josephson Charge Qubits

The total Hamiltonian contain three parts: the Hamiltonian of the system, environ- ment (bath), and their interaction

H_{T} = H_{S}+ H_{bath}+ H_{int}. (2.1)

Since our theoretical model is based on the experimental setup in [1, 5], we follow the form of their Hamiltonian. So the system Hamiltonian of two capacitively coupled charge qubits is

H_{S} =

N

X

n1,n2=NL

{E_{c1}(n_{1}− n_{g1}(t))^{2}+ E_{c2}(n_{2}− n_{g2}(t))^{2}
+E_{m}(n_{1}− n_{g1}(t))(n_{2}− n_{g2}(t))}|n_{1}, n_{2}ihn_{1}, n_{2}|

−1
2E_{j1}

N −1

X

n1=NL

N

X

n2=NL

{(|n_{1}+ 1ihn_{1}| + |n_{1}ihn_{1}+ 1|) ⊗ |n_{2}ihn_{2}|}

−1
2E_{j2}

N

X

n1=NL

N −1

X

n2=NL

{|n_{1}ihn_{1}| ⊗ (|n_{2}+ 1ihn_{2}| + |n_{2}ihn_{2}+ 1|)}, (2.2)

where n_{1}, n_{2} denote the charge states of the two qubits, and n_{g1}, n_{g2} are the gate
charges, NLmeans the lowest charge state, and the N is the highest charge state. We
choose the values of the parameters according to the experimental setup in Refs.[1, 5]

as

Ec1

h = 140.2 GHz,

Ec2

h = 162.2 GHz,

Ej1

h = 10.9 GHz,

Ej2

h = 9.9 GHz,

Em

h = 23.0 GHz.

(2.3)

The control parameters are

ngν(t) = n^{0}_{gν}+ ngν(t) (2.4)

with ν ∈ {1, 2}. The values of n^{0}_{g1} = 0.24, n^{0}_{g2} = 0.26 are chosen according to
Refs.[1, 5].

We consider two cases for the Hamiltonian Eq. (2.1). The first case is to treat
each of the qubits as a two-level system with N = 1, and N_{L}= 0. The second case is
to include the leakage state for the qubits with N = 2, and N_{L}= −1. Also, one can
split Eq. (2.2) into

H_{S} = H_{d}+ n_{g1}(t) × H_{1}+ n_{g2}(t) × H_{2}, (2.5)

where H_{d}, H_{1}, H_{2} are time independent. In the case of N = 1, N_{L} = 0, the Hamilto-
nians may be written in terms of Pauli matrices σ_{i} as

H_{d} = −(Em

4 +Ec1

2 )(σ_{z}^{(1)}⊗ I) − Ej1

2 (σ_{x}^{(1)}⊗ I)

−(E_{m}

4 +E_{c2}

2 )(I ⊗ σ^{(2)}_{z} ) − E_{j2}

2 (I ⊗ σ^{(2)}_{x} ) (2.6)

+E_{m}

4 (σ_{z}^{(1)}⊗ σ_{z}^{(2)}),

H_{1} = E_{c1}(σ_{z}^{(1)}⊗ I) +E_{m}

2 (I ⊗ σ_{z}^{(2)}) (2.7)

H_{2} = Em

2 (σ_{z}^{(1)}⊗ I) + E_{c2}(I ⊗ σ_{z}^{(2)}) (2.8)
where σ^{(ν)}_{i} are Pauli matrices acting on the Hilbert space of the νth qubit. In the
case of N = 2, and N_{L} = −1, we choose to use the original form of Hamiltonian in
number basis states. After appropriate identity shift

HS0

= HS− {...}I (2.9)

and H_{1} = ^{∂H}_{∂n}^{S}^{0}

g1, H_{2} = ^{∂H}_{∂n}^{S}^{0}

g2, we write H_{d} = H^{0}_{S} − n_{g1}(t)H_{1} − n_{g2}(t)H_{2}, and it’s

possible to get the form as

H_{1} = E_{c1}(S_{1}) + E_{m}

2 (S_{2}), (2.10)

H_{2} = Em

2 (S_{1}) + E_{c2}(S_{2}). (2.11)

where S1 and S2 are system operators similar to those in Eqs. (2.7) and (2.8).

In order to describe the decoherence effect, we focus on the charge noise and the interaction between the system and the bath is

Hint = δng1S1+ δng2S2 (2.12)

The symbol δn_{g} denotes the fluctuation of the control ,like Eq.(1.10), or other forms of
charge noise. Each of the qubit couples to its own independent charge noise environ-
ment or independent voltage fluctuation source. We model each of the independent
environments as a large set of harmonic oscillators, each of which interacts weakly
with its corresponding qubit. So the bath Hamiltonian is

H_{bath} =

2

X

ν=1

X

i

[ 1

2m^{(ν)}_{i} (p^{(ν)}_{i} )^{2}+m^{(ν)}_{i} ω_{i}^{(ν)}

2 (x^{(ν)}_{i} )^{2}] (2.13)

with

x^{(1)}_{i} = 1
q

2m^{(1)}_{i} ω^{(1)}_{i}

(a_{i}+ a^{†}_{i}), x^{(2)}_{i} = 1
q

2m^{(2)}_{i} ω_{i}^{(2)}

(b_{i}+ b^{†}_{i}). (2.14)

The bath operators δn_{g1} and δn_{g2} are

δn_{g1} ≡X

i

λ_{i}(a_{i}+ a^{†}_{i}) ≡ Y_{1}⊗ I = Γ_{1}, (2.15)

δng2 ≡X

j

λj(bj+ b^{†}_{j}) ≡ I ⊗ Y2 = Γ2, (2.16)

where Y_{1} and Y_{2} are bath operators acting on the Hilbert space of the first and

second qubits, respectively. The system-bath coupling Hamiltonian is then

H_{int} = S_{1}⊗ Γ_{1}+ S_{2}⊗ Γ_{2}. (2.17)

The master equation for the reduced density matrix ˜ρ(t) of the two qubit system in the interaction picture can be written as [23]

˙˜

ρ(t) = −X

i,j

ˆt

0

dt^{0}{[ ˜S_{i}(t) ˜S_{j}(t^{0}) ˜ρ(t^{0}) − ˜S_{j}(t^{0}) ˜ρ(t^{0}) ˜S_{j}(t)] < ˜Γ_{i}(t)˜Γ_{j}(t^{0}) >_{R}

+[ ˜ρ(t^{0}) ˜S_{j}(t^{0}) ˜S_{i}(t) − ˜S_{i}(t) ˜ρ(t^{0}) ˜S_{j}(t^{0})] < ˜Γ_{j}(t^{0})˜Γ_{i}(t) >_{R}}, (2.18)

where

˜

ρ(t) = e^{~}^{i}^{´}^{0}^{t}^{H}^{S}^{(τ )dτ}ρ(t)e^{−}^{i}^{~}^{´}^{0}^{t}^{H}^{S}^{(τ )dτ} (2.19)

Transforming back to the Schrodinger picture with

ρ(t) = e^{−}^{~}^{i}^{´}^{0}^{t}^{H}^{S}^{(τ )dτ}ρ(t)e˜ ^{i}^{~}^{´}^{0}^{t}^{H}^{S}^{(τ )dτ} (2.20)

we obtain

˙

ρ(t) = −i

~H_{S}(t) + e^{−}^{~}^{i}^{´}^{0}^{t}^{H}^{S}^{(τ )dτ}ρe˙˜ ^{~}^{i}^{´}^{0}^{t}^{H}^{S}^{(τ )dτ} + e^{−}^{i}^{~}^{´}^{0}^{t}^{H}^{S}^{(τ )dτ}ρ˜i

~H_{S}(t)e^{~}^{i}^{´}^{0}^{t}^{H}^{S}^{(τ )dτ}

= −i

~

[H_{S}(t), ρ(t)] + U^{†}(t) ˙˜ρU (t)

= −i

~

[H_{S}(t), ρ(t)]

− 1

~^{2}U^{†}(t)

2

X

i=1

ˆt

0

dt^{0}{[ ˜S_{i}(t) ˜S_{i}(t^{0}) ˜ρ(t^{0}) − ˜S_{i}(t^{0}) ˜ρ(t^{0}) ˜S_{j}(t)] < ˜Γ_{i}(t)˜Γ_{j}(t^{0}) >_{R}
+[ ˜ρ(t^{0}) ˜S_{i}(t^{0}) ˜S_{i}(t) − ˜S_{i}(t) ˜ρ(t^{0}) ˜S_{i}(t^{0})] < ˜Γ_{i}(t^{0})˜Γ_{i}(t) >_{R}}U (t)

= −i

~

[H_{S}(t), ρ(t)]

− 1

~^{2}

2

X

i=1

{SiU^{†}(t)
ˆt

0

dt^{0}( ˜Si(t^{0}) ˜ρ(t^{0}) < ˜Γi(t)˜Γi(t^{0}) >R−˜ρ(t^{0}) ˜Si(t) < ˜Γi(t^{0})˜Γi(t) >R)

−U^{†}(t)
ˆt

0

dt^{0}( ˜S_{i}(t^{0}) ˜ρ(t^{0}) < ˜Γ_{i}(t)˜Γ_{i}(t^{0}) >_{R} −˜ρ(t^{0}) ˜S_{i}(t^{0}) < ˜Γ_{i}(t^{0})˜Γ_{i}(t) >_{R})U (t)S_{i}}

= −i

~[HS(t), ρ(t)] − 1

~^{2}

2

X

i=1

{[Si, D^{(i)}(t)] − [Si, (D^{(i)}(t))^{†}]}. (2.21)

Since we are dealing with two independent bathes, so there is no cross term. If there
is no environment effect, there will only be the first term left in Eq.(2.21). The 2nd
and 3rd terms describe the environment effect and the dissipators D^{(i)}(t) are defined
as [24]

D^{(i)}(t) = − 1

~^{2}U^{†}(t)
ˆ t

0

dt^{0}U (t^{0})S_{i}ρ(t^{0})U^{†}(t^{0}) < ˜Γ_{i}(t)˜Γ_{i}(t^{0}) >_{R} U (t)

= − 1

~^{2}
ˆ _{t}

0

U_{S}(t − t^{0})S_{i}ρ(t^{0})C(t − t^{0})dt^{0} (2.22)

The definition of the bath correlation function C(t − t^{0}) =< ˜Γ_{i}(t)˜Γ_{i}(t^{0}) >_{R} will be
discussed in the next section. To overcome the difficulty of solving the time-nonlocal
integral-differential equations (2.21) and (2.22), we follow the approach in [24] to
introduce auxiliary density matrices in the extended Liouville space. The procedure
is to use many exponential terms to fit the bath correlation function, and the total

number of the exponential terms is determined by the accuracy of the fitting error.

Usually, we choose error to be less than 10^{−5}. So we can write

C(t − t^{0}) = C(τ ) =X

j

Cj(t − t^{0}) =X

j

Cj(0)e^{γ}^{j}^{τ}, (2.23)

where C_{j}(0) and γ_{j} are complex. In this case, the dissipators,

D^{(i)}(t) =X

j

K_{j}^{(i)}(t), (2.24)

where

K_{j}^{(i)}(t) = − 1

~^{2}
ˆ t

0

U_{S}(t − t^{0})S_{i}ρ(t^{0})C_{j}(t − t^{0})dt^{0}

= − 1

~^{2}
ˆ t

0

A × C_{j}(t − t^{0})dt^{0} (2.25)

with A = U_{S}(t − t^{0})S_{i}ρ(t^{0}). Differentiating Eq.(2.25) with respect to time, we obtain
d

dtK_{j}^{(i)}(t) = − 1

~^{2}U_{S}(0)S_{i}ρ(t)C_{j}(0) − 1

~^{2}
ˆ t

0

d

dt{U_{S}(t − t^{0})S_{i}ρ(t^{0})C_{j}(t − t^{0})}dt^{0}

= − 1

~^{2}C_{j}(0)S_{i}ρ(t) − 1

~^{2}
ˆ _{t}

0

dt^{0}{[−i

~H_{s}(t)A + i

~AH_{s}(t)]C_{j}(t − t^{0}) + AdC_{j}(t − t^{0})

dt }

= − 1

~^{2}

C_{j}(0)S_{i}ρ(t) − 1

~^{2}
ˆ _{t}

0

dt^{0}{−i

~

[H_{s}(t), A]C_{j}(t − t^{0}) + AdC_{j}(t − t^{0})

dt }

= − 1

~^{2}

C_{j}(0)S_{i}ρ(t) − 1

~^{2}

{L_{s}(t)K_{j}^{(i)}(t) + γ_{j}K_{j}^{(i)}(t)}. (2.26)

Eq.(2.21) with Eq.(2.24) and Eq.(2.26) form a set of coupled time-local differential
equations. Now with the above coupled time-local equations, one may solve the
evolution of the model. Remember that one can always returns to the Schrodinger
(von-Neumann) equation of the closed system by setting D = 0 in Eq.(2.21). Since
ρ(t) and K_{j}(t) are coupled and mutually needed for solving the differential equations.

One could combine them into a new vector as

−

→ρ (t) =

ρ(t)
K_{j}^{†}(t)
K_{j}(t)

...

. (2.27)

In this way, one may rewrite the coupled equations of motion in a very concise form as

d dt

−

→ρ (t) = Λ(t)−→ρ (t), (2.28)

## 2.3 Spectral Density and Correlation Function

Generally speaking, there are two processes relevant for the observed decoherence:

energy relaxation (the decay of probability of some states) and dephasing effect (the decay of the off-diagonal terms of the density matrix). We will discuss the bath correlation function and the bath spectral density in the following

In the interaction picture, the bath operator takes the form

Γ˜_{1}(t) = e^{~}^{i}^{H}^{B}^{t}Γ_{1}e^{−}^{~}^{i}^{H}^{B}^{t}. (2.29)

So the bath correlation function becomes

C(t − t^{0}) = C(τ ) = T r_{B}[˜Γ_{1}(t)˜Γ_{1}(t^{0})R_{0}]

= T r_{B}[(X

i

λ_{i}(a_{i}e^{−iω}^{i}^{t}+ a^{†}_{i}e^{iω}^{i}^{t}))(X

j

λ_{j}(a_{j}e^{−iω}^{j}^{t}^{0} + a^{†}_{j}e^{iω}^{j}^{t}^{0}))R_{0}]

= X

i

|λ_{i}|^{2}{(T r_{B}[a_{i}a^{†}_{i}R_{0}]e^{−iω}^{i}^{(t−t}^{0}^{)}+ T r_{B}[a^{†}_{i}a_{i}R_{0}]e^{iω}^{i}^{(t−t}^{0}^{)})

= ˆ ∞

0

dω^{0}J (ω^{0})[(n(ω^{0}) + 1)e^{−iω}^{i}^{(t−t}^{0}^{)}+ n(ω)e^{iω}^{i}^{(t−t}^{0}^{)}]

= ˆ ∞

0

dω^{0}J (ω^{0})[(n(ω^{0}) + 1){cos[ω^{0}(t − t^{0})] − i sin[ω^{0}(t − t^{0})]}

+n(ω){cos[ω^{0}(t − t^{0})] + i sin[ω^{0}(t − t^{0})]}]

= ˆ ∞

0

dω^{0}J (ω^{0}){(2n(ω^{0}) + 1) cos[ω^{0}(t − t^{0})] − i sin[ω^{0}(t − t^{0})]}

= ˆ ∞

0

dω^{0}J (ω^{0}){ 1 + e^{−}^{kB T}^{~ω}
1 − e^{−}^{kB T}^{~ω}

!

cos[ω^{0}(t − t^{0})] − i sin[ω^{0}(t − t^{0})]}

= ˆ ∞

0

dω^{0}J (ω^{0}){cos[ω^{0}(t − t^{0})] coth( ~ω^{0}

2k_{B}T) − i sin[ω^{0}(t − t^{0})]},
where J (ω^{0}) is the bath spectral density. Depending on the choice of the bath spectral
density, instead of the frequency integration limits of 0 and ∞, we may introduce the
lower bound and upper bound frequencies, also called infrared and ultraviolet cutoff
frequencies, ω_{r} and ω_{c}, respectively. Rewrite the bath correlation function concisely:

< ˜Γ_{1}(τ )˜Γ_{1}(0) > = < (X

i

λ_{i}(a_{i}e^{−iω}^{i}^{τ}+ a^{†}_{i}e^{iω}^{i}^{τ}))(X

j

λ_{j}(a_{j}+ a^{†}_{j})) >_{B}

=
ˆ _{ω}_{c}

ωr

dω^{0}J (ω^{0})(cos(ω^{0}τ ) coth( ~ω^{0}

2k_{B}T) − i sin(ω^{0}τ )) (2.30)
The noise power spectrum is defined as the Fourier transform of the symmetric bath

correlation function:

S_{U}(ω) = 1
2π

ˆ ∞

−∞

< Γ_{1}(τ )Γ_{1}(0) + Γ_{1}(0)Γ_{1}(τ ) >_{B} e^{−iωτ}dτ

= 1

2π ˆ ∞

−∞

[ ˆ ∞

0

dω^{0}J (ω^{0}) cos(ω^{0}τ ) coth( ~ω^{0}

2k_{B}T) × 2]e^{−iωτ}dτ

= 1

2π ˆ ∞

0

dω^{0}J (ω^{0}) coth( ~ω^{0}

2k_{B}T) × 2 × 1

2{πδ(ω − ω^{0}) + πδ(ω^{0} − ω)}

= 1

2π · 2πJ(ω) coth( ~ω
2k_{B}T)

= J (ω) coth( ~ω

2k_{B}T). (2.31)

The result of the following calculation has been used to derive Eq.(2.31) ˆ ∞

−∞

cos(ω^{0}τ )e^{−iωτ}dτ =
ˆ ∞

−∞

(e^{iω}^{0}^{τ} + e^{−iω}^{0}^{τ}

2 )e^{−iωτ}dτ

= 1 2[

ˆ ∞

−∞

e^{i(ω}^{0}^{−ω)τ} + e^{−i(ω}^{0}^{+ω)τ}dτ ]

= 1 2[

ˆ ∞ 0

e^{i(ω}^{0}^{−ω)τ}dτ −
ˆ −∞

0

e^{i(ω}^{0}^{−ω)τ}dτ
+

ˆ ∞ 0

e^{−i(ω}^{0}^{+ω)τ}dτ −
ˆ −∞

0

e^{−i(ω}^{0}^{+ω)τ}dτ ]

= 1 2[

ˆ ∞ 0

e^{i(ω}^{0}^{−ω)τ}dτ +
ˆ ∞

0

e^{i(ω}^{0}^{−ω)τ}dτ
+

ˆ ∞ 0

e^{−i(ω}^{0}^{+ω)τ}dτ +
ˆ ∞

0

e^{−i(ω}^{0}^{+ω)τ}dτ ]

= 1

2{[πδ(ω − ω^{0}) + i P

ω − ω^{0}] + [πδ(ω^{0}− ω) + i P
ω^{0} − ω]
+[πδ(ω^{0} + ω) + i P

ω^{0}+ ω] + [πδ(−ω^{0}− ω) + i P

−ω^{0}− ω]}

= 1

2{πδ(ω − ω^{0}) + πδ(ω^{0}− ω)}

The decoherence effect is discuss in a number of theoretical papers and the re- laxation of the excited states of the charge qubits off the degeneracy point had been discussed. The decoherence effects still lack of systematic studies. Especially the classical 1/f noise, we can only use it as a experimental data.

Figure 2.1 is the experimentally measured noise spectrum reported in [8]. The

Figure 2.1: The noise spectrum of the experimental data from [8]. (a) The two dashed
lines are Ohmic noise S_{U}^{Q} and 1/f noise S_{U}^{C} . (b) Simple diagram and the cross point
ωC. ωC is renamed as ωcross in the texts since it has be used to denote the cutoff
frequency.

measured samples. The pattern of the S4E can be seen as a combination of two
linear lines representing as Ohmic noise and 1/f noise, respectively. The cross point
of the two lines indicate that the noise contributions from 1/f and f dependencies
are the same at ω_{cross} = 2π × 2.6GHz. This also defines an effective temperature
Tcross = ^{~ω}_{k}^{cross}

B = 120mK, close to the electron temperature. Theoretical studies to construct a model to explain the measured noise spectra have been reported.

Instead of integrating the correlation function from zero to infinite, we choose
a ultraviolet cutoff frequency and an infrared cutoff frequency as the upper bound
and lower bound for the integration limits. The Ohmic noise spectral density is
proportional to the frequency, so it dominates over the 1/f noise at high frequencies
or low temperatures as shown below. For ~ω 2k^{B}T

SU(ω) = Jf(ω) coth( ~ω

2k_{B}T) ≈ Jf(ω)

= S_{U}^{Q}(ω) = 4e^{2}

π R~ω (2.32)

The symbol Q of S_{U}^{Q}(ω) means quantum. The other noise is the 1/f noise. Just like
its name, the 1/f noise is proportional to the inverse of the frequency of the bath

oscillators. For ~ω 2kBT

S_{U}(ω) = J_{1/f}(ω) coth( ~ω

2k_{B}T) ≈ J_{1/f}(ω)2k_{B}T

~ω

= S_{U}^{C}(ω) = (4E_{C}
e )^{2} α

2ω (2.33)

with α ≈ (1.3 × 10^{−3}e)^{2}. The symbol C of S_{U}^{C}(ω) means classical. So the 1/f noise
dominates at low frequencies or high temperatures as shown below. One may use the
following equations to get the above equations.

coth(x) = 1

tanh(x) = cosh(x)

sinh(x) = e^{x}+ e^{−x}

e^{x}− e^{−x} = 1 + e^{−2x}
1 − e^{−2x}

x 1 ⇒ coth(x) ≈ (1 + e^{−2x})(1 + e^{−2x}) ≈ 1 + 2e^{−2x}+ e^{−4x} ≈ 1

x 1 ⇒ e^{−2x} = 1 + (−2x) + (−2x)^{2}
2! + . . .

∴ coth(x → ∞) ≈ 2 2x = 1

x

## 2.4 Diagram of simulation

The cross point of the ohmic and 1/f power spectrum is around ω_{cross} = 2π ×2.6GHz.

So we define that the “high” and “low” frequency cutoff according to ω_{cross}, i.e.

high (ωc> ωcross) low (ωc< ωcross)

We plot the simulated bath correlation functions using the parameters of T =
10 mK, ω_{c}= 2π × 1 GHz, t_{f} = 55 ps, and dt = 0.1 ps.

As one can see in Figure 2.3, the value of the 1/f noise correlation function is bigger than the Ohmic correlation function when the cutoff frequency is low. Also, the shape of the Ohmic correlation function is steeper and the 1/f noise remains flat.

The blue lines in Figure 2.3 and Figure 2.4 indicate the real parts of the correlation

**10**^{−1}**10**^{0}**10**^{1}**10**^{2}**10**^{6}

**10**^{7}**10**^{8}**10**^{9}

ω

2π(GHz)

S(ω) ¯h2(GHz)

**ohmic (exp)**
**1/f (exp)**

Figure 2.2: The diagram of power spectrum from experimental data [8]. The cross
point is ω_{cross} = 2π × 2.6 GHz.

functions and the green lines indicate the imaginary parts of the correlation functions.

On the contrary, the bath spectral density with a high cutoff frequency gives
different responses. We plot the bath correlations in Figure 2.4 for the parameters of
T = 10 mK, ω_{c}= 2π × 100 GHz, t_{f} = 55 ps, and dt = 0.1 ps.

Both the shapes of the Ohmic and the 1/f noise correlation functions are steeper than their corresponding cases at low cutoff frequency. However, the values of the Ohmic correlation functions dominates over the 1/f noise correlation function at short time (within the time scale of the gate operation time). Therefore, we consider the Ohmic correlation function in the high cutoff frequency region.

0 20 40 60

−0.01

−0.005 0 0.005 0.01 0.015

C (Ohmic)

t (ps)

C/(1 GHz)2

0 20 40 60

−0.5 0 0.5 1 1.5

C (1/f)

t (ps)

C/(1 GHz)2

Figure 2.3: Correlation functions (blue line: real part; green line: imaginary part) for low cutoff frequency. The parameters used are T = 10 mK, ωc = 2π × 1 GHz, tf = 55 ps, and dt = 0.1 ps.

0 20 40 60

−100

−50 0 50 100 150

C (Ohmic)

t (ps)

C/(1 GHz)2

0 20 40 60

−0.5 0 0.5 1 1.5 2

C (1/f)

t (ps)

C/(1 GHz)2

Figure 2.4: Correlation functions (blue line: real part; green line: imaginary part) for
high cutoff frequency. The parameters used are T = 10 mK, ω_{c} = 2π ×100 GHz, t_{f} =
55 ps, and dt = 0.1 ps..

Chapter 3

Optimal Control

## 3.1 Introduction

The optimal control method can be applied to many different fields. But they all have
the same purpose that is to get a extreme accuracy result. There are many way to
do the optimal control such as GRAPE, Newton method, quasi-Newton method, and
Krotov method. The optimal control has an advantage to design a unitary quantum
gate in a closed system. In this thesis, we also use optimal control method to deal
with open systems, especially, non-Markovian open quantum systems. Some details
regarding the implementation of the optimal control are put in Appendices. Different
fidelity/error measures that are adopted as the objective function for optimization
and error values. We will discuss two commonly used error measures in Sec. 3.2
and will conclude in Chapter 4 why we use the error J_{2} defined in Eq. (3.2) for
the cases of open quantum systems. Since we are dealing with the gate operation,
we need to perform the operation independent of initial states. So I introduce the
state-independent optimal control and the superoperator concept latter.

## 3.2 Error/Fidelity Definition

According to [22], we can simply definition two kinds definition of error/fidelity.

J_{1} = 1 − 1

NT r(Q^{†}G(T )) (3.1)

J2 = 1

2NT r[(Q − G(T ))^{†}(Q − G(T ))] (3.2)

with Q is the target superoperator (CNOT gate), and G(T ) is the final optimized propagator in a superoperator form, and N is the dimension of the superoperator.

The concept of superoperator is presented in section 3.5.2 since it also relates to the concept of computational states.

In a closed system, the decoherence effect is not considered yet, the propagator should be unitary since the system is “perfect”. Thus in a closed system, we have

G(t)^{†}G(t) = I, (3.3)

and t ∈ [0, T ]. Expanding J_{2} , one can obtain the result that J_{2} = J_{1} in a “closed”

system.

However, we cannot promise Eq.(3.3) in a (non-Markovian) open system. And since G(T ) is no longer a unitary matrix, then the definition of J1may not be a better measure of gate error. If we consider leakage levels, we will not have a unitary result for the gate operation since we project the operation matrix onto the computational basis states and let the other elements be zeros.

So when we consider an open system or leakage levels, we will use the error
definition of J_{2}. However, to prove the optimal control based on the Krotov’ s method
in the next section. We will still use J1 as a part of the cost function.

Note that, when we replace the J_{1} into J_{2} in the cost function, the optimization
procedure can still work and the update rule just like Eq. (3.23) is obtained. However,
this change of target might cause the original derivation of Krotov unable to be done
and make the differential derivation more like the gradient ascent method. But the

step by step update is still more like a Krotov type. So even the original proof of Krotov might be a problem when we choose J2 as our error function, the result is still monotonically convergent if we choose the right penalty function λ.

## 3.3 One Qubit Optimal Control

There are many paper [9, 10, 11, 12, 13] describe the standard way of the Krotov optimal control method. However, some other papers use a different approach and obtain the same result,[14, 15]. Here we choose the most comprehensive parts in these approaches and combine them together to show the optimal iteration update equation.

Consider the equation of motion d

dt

−

→ρ (t) = Λ(t)−→ρ (t), (3.4)

where Λ(t) is a matrix operator (superoperator) acting on a column vector −→ρ (t). Let the propagator (superoperator) of −→ρ (t) be G(t). That is

−

→ρ (t) = G(t, 0)−→ρ (0) = G(t)−→ρ (0). (3.5)

So we can rewrite the equation of motion d

dtG(t)−→ρ (0) = Λ(t)G(t)−→ρ (0), (3.6)

and get d

dtG(t) = Λ(t)G(t), (3.7)

which will be Schrodinger equation in our control problem if the environment is
excluded. Suppose that the superoperator of the Hamiltonian Λ^{k}(t) has the form

Λ^{k}(t) = Λ + ^{k}(t)Λ , (3.8)