# Mitigating controller noise in quantum gates using optimal control theory

Aviv Aroch<EMAIL_ADDRESS>The Institute of Chemistry, The Hebrew
University of Jerusalem, Jerusalem 9190401, Israel  Shimshon Kallush  Sciences
Department, Holon Academic Institute of Technology, 52 Golomb Street, Holon
58102, Israel  Ronnie Kosloff  The Institute of Chemistry, The Hebrew
University of Jerusalem, Jerusalem 9190401, Israel

###### Abstract

All quantum systems are subject to noise from the environment or external
controls. This noise is a major obstacle to the realization of quantum
technology. For example, noise limits the fidelity of quantum gates. Employing
optimal control theory, we study the generation of quantum single and two-
qubit gates. Specifically, we explore a Markovian model of phase and amplitude
noise, leading to the degradation of the gate fidelity. We show that optimal
control with such noise models generates control solutions to mitigate the
loss of gate fidelity. The problem is formulated in Liouville space employing
an extremely accurate numerical solver and the Krotov algorithm for solving
the optimal control equations.

##  I Introduction

In reality, any quantum system is open. External intervention transforms the
unitary dynamics of a closed system into a non-unitary evolution of an open
one. We can classify three significant sources of such intervention: I. A
thermal environment. II. Back-action due to quantum measurement. III. Noise
originating from the external controller. This associated loss of coherence is
the main obstacle in achieving quantum technology [1, 2]. The aim of this
study is to explore the use of optimal control theory to mitigate the fidelity
destruction of a quantum gate due to controller noise.

We formulate the dynamics of open quantum systems by considering completely
positive trace-preserving maps (CPTP) [3]. In the Markovian limit, CPTP maps
are generated by the equation:

$\frac{\partial}{\partial t}\hat{\bm{\rho}}={\cal L}\hat{\bm{\rho}}$

Where the generator ${\cal L}$ has the Gorini, Kossakowski, Lindblad,
Sudarshan (GKLS) structure [4, 5]. We adopt this framework for the present
study.

Optimal control theory (OCT) has addressed open-system control problems [6, 7,
8, 9, 10, 11, 12, 13, 14]. In this case the generator of the dynamics ${\cal
L}$ becomes explicitly time-dependent a function of the control field ${\cal
L}(\epsilon(t))$. This observation raises the issue: Are the dissipative and
unitary dynamics linked? Does a change in the control field modify the
dissipative dynamics? Consistency with thermodynamics imposes additional
restrictions on the form of the master equation of a thermal environment [15].
In this case, the dissipative and the unitary parts are linked. Employing a
thermodynamically consistent master equation, optimal control of quantum
state-to-state and quantum gates was studied [16]. It was demonstrated that
both entropy-changing and preserving transformations can be achieved. In
particular, the negative aspect of dissipation on the fidelity of quantum
gates could be mitigated. The mechanism identified can be described as active
refrigeration accompanying the control process leading to a unitary quantum
gate. A common approach to mitigating environmental noise is dynamical
decoupling [17, 11]. Nevertheless, this method is not immune to noise on the
controller.

In the present study, we deal with a different dissipative mechanism noise
originating from the controller. This immediately imposes a link between the
dissipative and unitary dynamics. Previous studies have shown that such noise
leads to loss of purity and is therefore harmful to quantum control [18, 19].
Specifically, we address quantum systems with amplitude and phase noise, [20].

Our primary concern is exploring quantum control strategies to mitigate noise
and maintain high-fidelity gates. We address this issue by employing optimal
control theory first to generate unitary quantum gates in a pure Hamiltonian
setup. We study both single and two-qubit gates forming a universal
computation set. High-fidelity unitary gates are obtained. We now add noise to
the controller and observe the degradation of fidelity. On this reference
gate, we apply open system control, including noise, to mitigate the influence
of this noise as much as possible. The paper is arranged as follows: Section
II describes the noise model and the optimal control equations. Section III
describes the results, which are discussed and summarized in the concluding
Sec. IV.

##  II Quantum noisy dynamics

The noisy dynamics are modeled within open quantum systems [21]. Expressly, we
assume a quantum dynamical semi-group describes Markovian dynamics. The
generator of such dynamics has the form of Gorini, Kossakowski, Lindblad, and
Sudarshan master equation (GKLS) [5, 22]:

$\frac{d}{dt}\hat{\bm{\rho}}(t)=\mathcal{L}\hat{\bm{\rho}}=-\frac{i}{\hbar}[\hat{\mathbf{H}},\hat{\bm{\rho}}]+{\mathcal{D}}\hat{\bm{\rho}}\\\
$ (1)

$\hat{\mathbf{H}}$ is the systems Hamiltonian, ${\cal D}$ is the dissipative
generator. $\mathcal{L}$, and $\mathcal{D}$ are termed superoperators as they
map operators in Hilbert space to other operators.

###  II.1 Noise model

We consider a quantum system controlled by an external field. The Hamiltonian
of the system:

$\hat{\mathbf{H}}~{}=~{}\hat{\mathbf{H}}_{0}+\varepsilon(t)\hat{\mathbf{H}}_{c}$
(2)

where $\hat{\mathbf{H}}_{0}$ is the drift Hamiltonian and
$\hat{\mathbf{H}}_{c}$ is the control Hamiltonian and $\varepsilon(t)$
describes the time dependent control field. We aim to study the influence of
noise originating from the controller (The pulse generator). We differentiate
between two models of control noise: amplitude and phase noise. For amplitude
noise, we consider random fluctuations in the amplitude of the control field
$\varepsilon(t)=\varepsilon_{c}(t)+\xi(t)$. Assuming the controller is fast
relative to the system dynamics We chose to model the noise as a Gaussian
random process
$\langle\xi(t)\xi(t^{\prime})\rangle=\sqrt{2\gamma}\delta(t-t^{\prime})$.
Under these conditions, the dissipative part of Eq. (1) is generated from
$\hat{\mathbf{H}}_{c}$ [4, 23, 24]:

$\begin{array}[]{l}{\cal
D}_{A}=-{\gamma}_{A}\varepsilon_{c}^{2}(t)[\hat{\mathbf{H}}_{c},[\hat{\mathbf{H}}_{c},\bullet]]\\\
\end{array}\ $ (3)

where $\gamma$ is the noise strength. Phase noise originates from timing
errors in the controller [25, 26] or fluctuations in frequency [27]. Assuming
a Gaussian white noise model, the dephasing generator becomes [25]:

$\begin{array}[]{l}{\cal
D}_{P}=-{\gamma}_{P}[\hat{\mathbf{H}},[\hat{\mathbf{H}},\bullet]]\\\
\end{array}\ $ (4)

In both cases, the dissipator depends on the control field.

###  II.2 Vectoraizing Louivile space

We describe the dynamics employing the Hilbert space of system operators
defined by the scalar product.
$\left(\hat{\bm{A}}\cdot\hat{\bm{B}}\right)=\text{tr}\\{\hat{\bm{A}}^{\dagger}\hat{\bm{B}}\\}$.
Superoperators map operators into the same Hilbert space. We, therefore,
choose to describe operators as vectors and superoperators by matrix-vector
multiplications. We can now decompose a general operator $\hat{\bm{X}}$
employing a full orthogonal operator basis,

$\\{\hat{\bm{A}}\\}=[{\hat{\bm{A}}_{1}},\hat{{\bm{A}}_{2}},....,\hat{{\bm{A}}_{N}}]~{}~{}~{}~{}~{}~{}~{},~{}~{}~{}~{}\text{tr}\\{{\hat{\bm{A}}_{i}}^{\dagger}\cdot\hat{\bm{A}}_{j}\\}=\delta_{ij}$
(5)

Using this basis, we can represent $\hat{\bm{X}}$ as;

$\displaystyle\hat{\bm{X}}=\sum_{i}\chi_{i}\hat{\bm{A}}_{i}$ (6)

where $\chi_{i}$ is an element of a vector $\vec{x}$ of expansion coefficients
of size $N^{2}$ that characterizes $\hat{\bm{X}}$. Using this notation, we can
define any operator of the dimension $N\times N$ with a vector of size
$N^{2}$.

In this formalism the mapping of a super-operator
$\Lambda\hat{\bm{X}}=\hat{\bm{Y}}$ translates to a matrix-vector
multiplication: $\tilde{\bm{\mathcal{G}}}\vec{\chi}=\vec{\upsilon}$ where
$\hat{\bm{Y}}=\sum_{i}\upsilon_{i}\hat{\bm{A}}_{i}$. The dynamical equation
Eq. (1) can be reduced to a series of matrix-vector multiplications where
$\overrightarrow{\rho}$ will denote to $\hat{\bm{\rho}}$ in vector space, and
$\tilde{\bm{\mathcal{L}}}$ will denote the superoperator $\mathcal{L}$. For
Additional details, Cf. ref. [28]. ( Throughout this paper, we will use the
notation ” $\tilde{}$ ” to assign a super operator in Louiville space ). The
action of the superoperator on the density operator will read as follows:

$\displaystyle\mathcal{L}\hat{\bm{\rho}}~{}\Rightarrow~{}\tilde{\bm{\mathcal{L}}}\overrightarrow{\rho}$
(7)

The dephasing super-operators,$\cal{D}_{A}$ and $\cal{D}_{P}$, are defined by
the following equation:

$\displaystyle\tilde{\bm{\mathcal{D}}}_{A}=-{\gamma}_{A}\varepsilon_{c}^{2}(t)({\bm{\mathcal{H}}}_{c}{\bm{\mathcal{H}}}_{c}\bullet)$
(8)
$\displaystyle\tilde{\bm{\mathcal{D}}}_{P}=-{\gamma}_{A}({\bm{\mathcal{H}}}{\bm{\mathcal{H}}}\bullet)$

Where ${\tilde{\bm{\mathcal{H}}}}^{\prime}$ corresponds to the super-operator
$-\frac{i}{\hbar}[\hat{\mathbf{H}}^{\prime},\bullet]$ in Louiville vector
space.

More details on vectorization are described in Appendix A.

###  II.3 Optimal Control Theory of Open systems (OCT)

Our general task is to mitigate the influence of noise on a quantum gate. We
employ quantum optimal control theory (OCT) to provide methods to compute such
controls. The external fields $\left\\{\mathcal{E}_{k}\right\\}$ interacting
with the quantum system have the task of steering the system’s dynamics from
an initial state to a final state. A more ambitious task is to find the
driving field that generates a quantum map $\Lambda(T)$ [29]. The dynamical
equation of motion for the map becomes:

$\frac{d\Lambda(t)}{dt}=\mathcal{L}(t){\Lambda}(t)~{}~{}\Rightarrow~{}~{}\frac{d\tilde{\bm{\mathcal{G}}}}{dt}=\tilde{\bm{\mathcal{L}}}(t)\tilde{\bm{\mathcal{G}}},$
(9)

where $\tilde{\bm{\mathcal{G}}}$ is the evolution operator in the vector
space. The initial condition $\Lambda(t=0)=\mathbb{I}$, where $\mathbb{I}$ is
the identity super-operator. The objective is to obtain the optimal driving
field $\epsilon(t)$ that induces a given transformation $\hat{\mathrm{O}}$ at
$t=T$. This can be interpreted as the mapping of a complete set of operators
$\\{\hat{\bm{A}}\\}$ as close as possible to the desired map
$\hat{\mathrm{O}}$.

In optimal control theory (OCT), this task is translated to optimizing the
objective functional [14]:

$\displaystyle\mathcal{J}_{max}=\text{Tr}\\{\hat{\mathrm{O}}^{\dagger}\Lambda(T)\\}=\sum_{j}\text{tr}\\{(\hat{\mathrm{O}}\hat{\bm{A}}_{j})^{\dagger}(\Lambda(T)\hat{\bm{A}}_{j})\\}\Rightarrow\sum_{j}\text{tr}\\{(\tilde{\bm{\mathrm{O}}}\hat{\bm{A}}_{j})^{\dagger}(\tilde{\bm{\mathcal{G}}}(T)\hat{\bm{A}}_{j})\\}$
(10)

Where $\tilde{\bm{\mathcal{G}}}$ corresponds to $\Lambda$ in Hilbert space. We
use the symbol Tr for a trace of super-operators, while tr to trace operators
in Hilbert space. Two constraints are added to the objective. The first
restricts the dynamics to comply with the Liouville-von Noiman equation Eq.
(9),

${\mathcal{J}_{con}}=\int_{0}^{T}\operatorname{Tr}\left\\{\left(\frac{\partial\Lambda(t)}{\partial
t}-{\mathcal{L}}(t)\Lambda(t)\right)\Upsilon(t)\right\\}dt$ (11)

where $\Upsilon(t)$ is a super-operator Lagrange multiplier. The second
constraint restricts the total field energy

$\mathcal{J}_{penal}=\lambda\int_{0}^{T}\frac{1}{s(t)}|\epsilon(t)|^{2}dt$
(12)

$\lambda$ is a scalar Lagrange multiplier, and $s(t)$ is a shape function that
smoothly turns the pulse on and off. For this task, we choose a Gaussian
profile. Adding these constraints, we obtain the overall functional,

$\mathcal{J}_{Tot}=\left[\text{Tr}\\{\hat{{O}}^{\dagger}\Lambda_{T}\\}-\int_{0}^{T}\text{Tr}\left\\{\left(\frac{\partial\Lambda(t)}{\partial
t}-\hat{\mathcal{L}}(t)\Lambda(t)\right)^{\dagger}\Upsilon(t)\right\\}dt\right]\\\
-\lambda\int_{0}^{T}\frac{1}{s(t)}|\epsilon(t)|^{2}dt,$ (13)

The control task is translated to the maximization of the generalized
objective ${\mathcal{J}_{Tot}}$, $\delta\bar{J}=0$. Functional derivatives
with respect to the various functional elements are then taken,
$\Upsilon,\Lambda$, and $\epsilon$, resulting in the following set of control
equations:

  1. 1.

The Louivile equation Eq. (9) with the initial condition
$\Lambda(t=0)=\mathbb{I}$ for $\Lambda$.

  2. 2.

The inverse Louivile equation

$\frac{\partial\Upsilon(t)}{\partial
t}=\mathcal{L}^{*}(t)\Upsilon(t)~{}~{}\Rightarrow~{}~{}\frac{\partial\tilde{\bm{Y}}(t)}{\partial
t}=\tilde{\bm{\mathcal{L}}}^{\dagger}(t)\tilde{Y}(t)$ (14)

Where $\tilde{\bm{Y}}$ corresponds to $\Upsilon$ in the vector space, with the
condition $\Upsilon(t=T)=\hat{\mathrm{O}}^{\dagger}$. Note that the loss of
information in the reverse equation is also inverted so that dephasing occurs
from both time ends.

  3. 3.

The control field update equation for amplitude noise:

$\Delta\bm{\epsilon}(t)=-\operatorname{Im}\frac{s(t)[\operatorname{Tr}\\{\Upsilon(t)\mathcal{\hat{H}}_{c}\Lambda(t)\\}]}{2(\lambda+{\gamma}_{A}[\operatorname{Tr}\\{\Upsilon(t)\mathcal{\hat{H}}_{c}^{2}\Lambda(t)\\}])}~{}~{}\Rightarrow~{}~{}-{s(t)}\operatorname{Im}\frac{\text{T}r\\{\tilde{\bm{Y}}^{*}\tilde{\bm{\mathcal{H}}}_{c}\tilde{\bm{\mathcal{G}}}(t)\\}}{2(\lambda+{\gamma}_{A}Tr\\{\tilde{\bm{Y}}^{*}\tilde{\bm{\mathcal{H}}}_{c}^{2}\tilde{\bm{\mathcal{G}}}(t)\\})}$
(15)

$\frac{\partial}{\partial\epsilon}\mathcal{L}={\cal H}_{c}$ where
$\tilde{\bm{\mathcal{H}}}_{c}$ is a super operator and corresponds to
$\frac{i}{\hbar}[\hat{{H}}_{c},\bullet]$.

The dynamical equations are solved iteratively (counted by k) employing the
Krotov method [30, 29], which leads to an update to $\bm{\epsilon}^{(k)}$:

$\Delta\bm{\epsilon}^{(k)}(t)=-\frac{s(t)}{2\lambda}\operatorname{Im}\left[\operatorname{Tr}\left\\{\Upsilon^{(k-1)}(t)\mathcal{L}_{c}\Lambda^{(k)}(t)\right\\}\right]\Rightarrow-\frac{s(t)}{2\lambda}\operatorname{Im}\left[\text{T}r\left\\{\tilde{\bm{Y}}^{(k-1)}(t)\tilde{\bm{\mathcal{H}}}_{c}\tilde{\bm{\mathcal{G}}}^{(k)}(t)\right\\}\right].$
(16)

This evaluation continues until certain fidelity $(1-\mathcal{J}_{max})$ is
reached. Note that for most cases $\lambda\gg\gamma_{A}$, thus one can neglect
its term at the denominator.

We solve numerically the dynamics in Eq. (9) and (14) employing a semi-global
propagation method [31] described in Appendix B. With these solutions, we
update Eq. (16) to update the field and repeat until the desired fidelity is
reached. The Krotov method is known to stagger after many iterations.
Therefore, we halt this iterative process when stagnation is reached.

##  III Results

Optimal control theory is employed first to obtain the desired unitary gate
without dissipation. This solution serves as a reference to study the
influence of noise on fidelity. At this point, optimal control is used again,
including the dissipation to search for control fields that mitigate the
impact of noise. The generator of the dynamical map has the form:

$\tilde{\bm{\mathcal{L}}}=\tilde{\bm{\mathcal{H}}}_{0}+\tilde{\bm{\mathcal{H}}}_{c}+\tilde{\bm{\mathcal{D}}}_{A/P}.$
(17)

Where $\tilde{\bm{\mathcal{D}}}$ generates noise where the indices ’A’ and ’P’
represent the two types of quantum noise defined at Eq. (3) and Eq. (4),
respectively.

We choose a driven system that is completely unitary and controllable. This
controllability results from commutators of the drift Hamiltonian,
$\hat{\mathbf{H}}_{0}$, and the control Hamiltonian, $\hat{\mathbf{H}}_{c}$,
which generate the full rank algebra [32].

Three representative target gates are studied:

  * •

A Hadamard gate on a single-qubit SU(2)

  * •

A CNOT gate in SU(2) space embedded in SU(4),

  * •

An entangling gate on a 2-qubit system SU(4).

The generator of the dynamical map is Eq. (17), and the cases differ in their
Hamiltonians and, as a consequence, their dephasing elements: Eq.’s (3) and
(4).

The first step in OCT is to obtain the reference unitary gate, which also
establishes a good guess of the control field $\varepsilon(t)$. A good guess
can become essential for convergence to high fidelity. In this case, the OC
dynamics is generated by:

$\tilde{\bm{\mathcal{L}}}=\tilde{\bm{\mathcal{H}}}_{0}+\tilde{\bm{\mathcal{H}}}_{c}.$
(18)

We define the fidelity as
$F=\text{Tr}\\{\hat{\mathrm{O}}^{\dagger}\Lambda(T)\\}$, and infidelity,
$IF=1-\text{Tr}\\{\hat{\mathrm{O}}^{\dagger}\Lambda(T)\\}$. The unitary
infidelity is assigned by $IF_{U}$.

To test the effect of noise on the infidelity $IF_{U}$, we used the solution
obtained from $IF_{U}$. We propagated a system with noise, Eq. (17) as its
generator of the dynamical map.

In figure 1, the degradation of the fidelity by the noise from the infidelity
of $IF_{U}$ to $IF_{n}$. We plot the log of ratio $R=\frac{IF_{n}}{IF_{U}}$ to
the noise rate $\gamma$.

Figure 1: Destruction of fidelity by noise. The change in infidelity of a
quantum gate $IF=1-\text{Tr}\\{\hat{\mathrm{O}}^{\dagger}\Lambda(T)\\}$. The
reference is a converged unitary control without noise with infidelity
$IF_{U}$. The same field is used to obtain the gate with noise $IF_{n}$. The
log of the ratio $R=\frac{IF_{n}}{IF_{U}}$,
$\log(R)=\log(IF_{n})-\log(IF_{U})$ is displayed as a function of the noise
strength $\gamma$. Both phase and amplitude noise are displayed. The insert
displays a zoom on the low noise level. Left panel: The single qubit Hadamared
gate. Right panel: The two-qubit entangling gate.

Figure 1 shows that amplitude noise $\tilde{D}_{A}$ is more resilient than the
phase noise $\tilde{D}_{P}$, for the Hadamard gate. The degradation is
significant when the noise strength increases, saturating with three orders of
magnitude in infidelity. The two-qubit entangling gate is two orders of
magnitude more sensitive than the single-qubit gate with respect to the noise
strength $\gamma$.

###  III.1 Hadamard gate-SU(2) space

The Hadamard gate performs a rotation of $\pi$ about the axis
$(\hat{x}+\hat{z})/\sqrt{2}$ of the Bloch sphere. Expressed as a superoperator
in the operator basis $\hat{I},\hat{S}_{X},\hat{S}_{Y},\hat{S}_{Z}$, the
Hadamard gate is given by:

$\Lambda_{U}=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&0&0&-1\\\ 0&0&-1&0\\\
0&-1&0&0\end{array}\right)~{}~{}.$ (19)

We chose a Hamiltonian, allowing high fidelity for this type of target gate.
The drift Hamiltonian:

$\hat{\mathbf{H}}_{0}=u{\hat{S}}_{Z}+a_{X}{\hat{S}}_{X}$ (20)

and the control Hamiltonian:

$\hat{\mathbf{H}}_{c}=\varepsilon(t){\hat{S}}_{XY}$ (21)

where ${\hat{S}}_{XY}$ is a rotation with an axis in the $XY$ plane (it was
found that any direction in this plane leads to control). With this generator,
we first obtained the reference unitary gate by OCT, with unitary infidelity
$IF_{U}=1\times 10^{-5}$. Using the optimal unitary control field as a guess,
solutions under noise are obtained as a function of $\gamma$, the noise
amplitude. For each value of $\gamma$, we apply OCT to mitigate the influence
of noise. Figure 2, shows the infidelity improvement as a function of the
noise rate, $\gamma$ for a short control period (2-Rabi cycles
$\tau=4\pi/\Omega$ , $\Omega=\sqrt{u^{2}+a_{X}^{2}}$). The improvement ratio
is defined as $NC=\frac{IF_{n}}{IF_{F}}$, where $IF_{n}$ is the infidelity
obtained using the unitary guess field subject to noise and $IF_{F}$ is the
infidelity brought by optimal control mitigating the noise.

Figures 3 shows the infidelity improvement for a longer control period
(6-Rabi-cycles). When increasing the time allocated for control as shown in
Fig. 3, the reference infidelity was maintained to $IF_{U}\sim 1\times
10^{-5}$. However, to keep comparable infidelity for both cases, the dephasing
rate $\gamma$ had to be reduced by two orders of magnitude. In both cases, we
find that OCT mitigates the influence of noise at low values of $\gamma$. When
the dissipation increases, the ability to mitigate its influence vanishes.

Figure 2: Noise reduction via Optimal Control for the Hadamard gate for short
times (2-Rabi cycles): The noise reduction is defined as
$NC=\frac{IF_{n}}{IF_{F}}$ Vs the dephasing rate, $\gamma$ in logarithmic
scale. We show here the two dephasing mechanisms, Amplitude noise
,$\tilde{D}_{A}$, (red,A) and Phase noise, $\tilde{D}_{P}$, (blue, P). (inset)
The accumulated energy of the control field vs. the dephasing rate in the two
dephasing mechanisms.

Infidelity can be mitigated with OC by a factor of two in low noise rates
($<$0.2$\times 10^{-3}$) for phase noise and on amplitude noise systems in
high noise rate ($\gamma>$0.2$\times 10^{-3}$) systems.

The optimal control formulation restricts the accumulated power Cf. Eq. (12).
Nevertheless, we observe that the invested pulse energy in the optimal
solution varies with the noise strength $\gamma$. The inset of figure 2 shows
this dependence. Optimal control must reduce the field strength and,
consequently, the controller’s noise to battle phase noise at slow noise
rates. However, more substantial noise causing significant dephasing can be
controlled more effectively only by boosting the field’s strength. This is not
the case for amplitude noise ($\tilde{D}_{A}$); to overcome the noise, OCT
reduces the strength of the field and, consequently, the effect of the noise.
(Cf. Eq. (29) ).

Figure 3: Noise reduction using Optimal Control for the Hadamard gate at
longer times (6-Rabi cycles): The two dephasing mechanisms, amplitude noise
,$\tilde{D}_{A}$, (red,A) and Phase noise, $\tilde{D}_{P}$ (blue), (blue, P)
are shown. (inset) The accumulated energy of the control field vs. the noise
rate in the two noise mechanisms.

Optimal control mitigation can reduce $IF_{F}$ for the phase noise more
effectively than $IF_{F}$ of the amplitude noise at low noise rates. This
trend is reversed at high rates. The inset of figure 3 shows that the optimal
control solution reduces the field’s energy for higher dephasing rates,
resulting in higher fidelity.

Figure 4: The trajectory of an initial operator on the $Z$ direction during a
Hadamard transformation under phase noise. The reference unitary evolution in
blue has an infidelity of $IF_{U}=1\times 10^{-5}$. The trajectory with noise
in red does not reside on the surface of the Bloch sphere. The optimal control
solution in yellow mitigates the infidelity to a value of $IF_{n}=15\times
10^{-3}$ The right panel zooms on the target state. The outer sphere
represents the purity of 1. The inner-sphere purity of $90\%$.

Figure 4 follows a trajectory of an element of the operator basis during the
Hadamard transformation. When the generator is noiseless, the dynamics stay on
the surface of the Bloch sphere, and noise pushes the dynamics into the
interior of the Bloch sphere. The optimal control mitigation pushes back the
trajectory to the surface.

###  III.2 CNOT gate - SU(2) embedded in SU(4).

The two-qubit entangling and single-qubit gates comprise a universal
computational set that can simulate any quantum circuit. In the subsection, we
optimize the control field generating the $CNOT$ gate:

$\displaystyle\hat{U}=\left(\begin{array}[]{llll}1&0&0&0\\\ 0&1&0&0\\\
0&0&0&1\\\ 0&0&1&0\end{array}\right)$ (26)

Here, the optimization is carried only for the second qubit, embedded in the
two-qubit Hilbert space.

The gate can be realized employing a subspace SU(2) algebra within the SU(4)
algebra of the full two-qubit gates. This can be achieved by employing the
trivial drift generator:

$\tilde{\bm{\mathcal{H}}}_{0}=\tilde{\cal I}^{1}\otimes\tilde{\cal I}^{2}$
(27)

and the control generator:

$\tilde{\bm{\mathcal{H}}}_{c}=\varepsilon(t)\sum_{i=X,Y}a_{i}(\tilde{\cal
I}^{1}-{\tilde{S}^{1}}_{Z})\otimes{\tilde{S}^{2}}_{i}$ (28)

The two noise models $\tilde{D}_{A}$ and $\tilde{D}_{P}$ coincide in this
dynamics.

First, we used OCT on the pure unitary dynamics to get $IF_{U}$ converging to
infidelities of $\sim 1\times 10^{-5}$. Figure 5 shows the fidelity
degradation when the noise strength increases.

We now employ this solution as a guess for OCT protocol on noisy systems,
partially mitigating noise degradation.

Figure 5:  Infidelity as a function of the noise strength $\gamma$ for the
C-NOT gate. The dashed line (brown, $IF_{n}$): infidelity employing the
unitary control field. The solid line (purple, $IF_{F}$): infidelity of the
OCT solution mitigating noise. The star represents the infidelity $IF_{U}$ of
the unitary control. (inset) The noise reduction, NC (log scale), Vs.
,$\gamma$.

Figure 5 shows the degradation due to noise $IF_{n}$ relative to $IF_{U}$
(star) and the mitigation due to OCT ($IF_{F}$). The inset quantifies the
noise mitigation due to OCT.

The loss of fidelity of the gate can be caused by a misguided trajectory on
the Bloch surface or an unrecoverable loss of purity toward the interior. The
instantaneous purity loss

$\frac{d}{dt}\textrm{tr}\\{\hat{\bm{\rho}}^{2}\\}~{}=~{}2\textrm{tr}\\{\hat{\bm{\rho}}\bm{\mathcal{L}}\hat{\bm{\rho}}\\}~{}=~{}2\textrm{tr}\\{\hat{\bm{\rho}}\bm{\mathcal{D}}\hat{\bm{\rho}}\\}$
(29)

The second equality is the result that the unitary generator does not change
purity. For state-to-state control for pure dephasing, the optimal solution
prefers adiabatic following [33] where
$\bm{\mathcal{D}}_{P}\hat{\bm{\rho}}=0$. We must consider the purity loss for
all possible inputs for a gate. Thus, the optimal solution is more involved.

Phase noise is equivalent to the noise generated by imperfect timekeeping
[26], where the fidelity was found to scales as:
$F=({2+e^{-\frac{\theta^{2}}{2N}}})/{3}$, where $\theta$ is the pulse area and
$N\propto 1/\gamma$, is the ticking accuracy. When calculating the fidelity as
a function of $\gamma$, our results completely agree with this scaling.

For amplitude noise, the instantaneous purity loss becomes in vector notation:

$\frac{d}{dt}\textrm{tr}\\{\overrightarrow{\rho}^{2}\\}=2\left(\overrightarrow{\rho}\cdot\tilde{\bm{\mathcal{D}}}_{A}\overrightarrow{\rho}\right)=-{\gamma}_{A}\varepsilon_{c}^{2}(t)\left(\overrightarrow{\rho}\cdot{\bm{\mathcal{H}}}_{c}{\bm{\mathcal{H}}}_{c}\overrightarrow{\rho}\right)$
(30)

For a pure state, the rate of purity loss for the amplitude noise model
becomes [34]:

$\frac{d}{dt}\textrm{tr}\\{\hat{\bm{\rho}}^{2}\\}~{}=~{}-2\gamma\varepsilon^{2}(t)\Delta\hat{\mathbf{H}}_{c}$
(31)

where $\Delta\hat{\mathbf{H}}_{c}$ is the variance of the control operator
[18]. Examining Eq. (31), the optimal solutions will therefore minimize the
control amplitude as the amplitude of noise increases (Cf. Fig. 2 and Fig. 3
). For a gate transformation of the term
$\left(\overrightarrow{\rho}\cdot{\bm{\mathcal{H}}}_{c}{\bm{\mathcal{H}}}_{c}\overrightarrow{\rho}\right)$
in Eq. (30) was found to be almost constant during the propagation, when
averaged over a complete set of initial states.

###  III.3 Entangling gate-SU(4) space

The ultimate test of optimal control mitigation employs the full pallet of the
16 operator bases of SU(4). The dynamical generators connect all elements of
the algebra. Obtaining converged results was difficult. The entangling gate
set as a target is:

$\displaystyle\hat{U}=\left(\begin{array}[]{llll}1&0&0&0\\\ 0&1&0&0\\\
0&0&0&i\\\ 0&0&-i&0\end{array}\right)$ (36)

Employing many random initial guesses, the maximum fidelity obtained using a
single control field was only $\sim 80\%$. We, therefore, modified the drift
Hamiltonian such that the two qubits are not in resonance:

$\tilde{\bm{\mathcal{H}}}_{0}=a\tilde{\cal{I}}^{1}\otimes\tilde{\cal{I}}^{2}+\omega_{1}{\tilde{S}}_{Z}^{1}\otimes\tilde{\cal{I}}^{2}$
(37)

where $a$ is a phase factor and $\omega_{1}$ is the qubit frequency. In
addition, we increased the number of control fields so that the two control
generators are:

${{{}^{Z}}\tilde{\bm{\mathcal{H}}}_{c}}=\varepsilon_{Z}(t)\sum_{i=X,Y}a_{i}(\tilde{\cal
I}-{\tilde{S}}_{Z})^{1}\otimes{\tilde{S}}_{i}^{2}$ (38)
${{}^{E}}\tilde{\bm{\mathcal{H}}}_{c}=\varepsilon_{E}(t)({\tilde{S}}_{X}^{1}\otimes{\tilde{S}}_{Z}^{2})$
(39)

We tried several random guesses to imitate the control for this model, but
none produced high fidelity. We, therefore, employed as a guess field
$\varepsilon_{Z}$ obtained from the model of the previous sub-section. The
initial guess led to $IF_{U}\sim 0.05$ on the noiseless system. OC could
improve the infidelity to $IF_{U}\sim 1\times 10^{-6}$ when the Krotov
iterative process reached stagnation.

The CNOT+Phase gate can be part of the universal set. We, therefore, employ it
to test the influence of noise.

Figure 6: The noise cancellation $NC=\frac{IF_{n}}{IF_{F}}$ (log scale) Vs.
the noise rate $\gamma$ for the entangling gate under phase noise. The insert
focuses on the low noise regime.

The infidelity changes due to phase noise for the entangling SU(4) gate are
shown in figure 6. The noise significantly degrades the fidelity. Optimal
control can recover high fidelity. When the noise rate is low, we overcome the
saturation in the Krotov method and obtain infidelities better than the
unitary case (Cf. the inset) improving with noise rate. The reason is that
Krotov’s iterative procedure for systems with noise stagnates much later. New
control channels are introduced to the control protocol as suggested in [35].
For higher noise values, the mitigation due to OCT declines as expected.

The influence of amplitude noise on infidelity is shown in Figure 7. Optimal
control managed to restore $\sim 4$ orders of magnitude. As expected, this
mitigation decreases with the growth of the dephasing rate.

Figure 7: The noise cancellation $NC=\frac{IF_{n}}{IF_{F}}$ (log scale) Vs.
the noise rate $\gamma$ for entangling gate under amplitude noise.

Additional insight was obtained by using the optimal control field with noise
to a noiseless propagation. The resulting infidelity was $IF_{U}\sim 10^{-3}$,
3 orders of magnitude lower that the optimized $IF_{U}\sim 10^{-6}$. This
means that the optimal field with the noise represents a new active solution
that emerges due to the control under noise.

##  IV Conclusions

Quantum devices employ interference and entanglement as crucial resources.
Dissipation is the primary limiting factor to achieving quantum technology [1,
2]. Significant effort is therefore devoted to isolating the quantum device
from the environment. However, the noise originating from the controller is
typically ignored. Since such noise is unavoidable, we addressed the question:
Can optimal theory achieves the desired quantum gate while mitigating the harm
due to noise?

What is unique about controller noise is that it is ”hands-off; the damage is
avoided. Our Markovian noise model assumes the controller is fast. A GKLS
master equation describes the quantum dynamics of both amplitude and phase
noise. Extension of this study to non-Markovian noise model [24], and
Poissonian noise is feasible.

Employing advanced techniques of OCT for open systems, we first find a control
that generates a high-fidelity unitary gate. This reference gate studies the
degradation of fidelity due to the controller noise. Typically, the system is
more sensitive to phase noise. We then apply OCT again to mitigate the
fidelity loss due to noise. We have found that quantum control could find
control fields that restore fidelity or even surpass it.

The present study could not be possible without a highly accurate propagator
addressing explicit time dependence and non-hermitian operators. The gate
infidelity we demand is extremely small but significant. The differences
between the unitary and noisy solutions are minute.

Our findings emphasize the significance of error mitigation techniques
tailored to specific quantum gates and system characteristics. The varying
quality of fidelity improvement based on the target gate and system
underscores the need for targeted and adaptable strategies in the face of
noise.

Notably, the integration of Optimal Control Theory into the framework of open
quantum systems led us to discover novel control solutions distinct from
noiseless optimal solutions.

In conclusion, our study showcases the potential of Optimal Control Theory in
mitigating errors caused by phase and amplitude noise, elevating gate fidelity
to new heights within the domain of quantum computing.

## Acknowledgement

We thank Ido Schaefer for his help in implementing the semi-global propagator.
In addition we thank Raam Uzdin and Roie Dann for helpful discussions. This
study was supported by the Israel Science Foundation (grant nos. 510/17 and
526/21).

##  Appendix A Vectorizing Liouville space

The most common and expensive numerical operation in any propagation method is
the matrix multiplication of a state vector. The high precision propagators
presented in the following appendix contain an even higher order of
operations: super-operators acting of operators. Thus, it is more convenient
to vectorize the Liouville space to stay closer to the regular matrix-vector
operation. A description of the vectorizing procedure employed in both
analytical and numerical aspects. A Hilbert space composed of operators can be
generated by defining a scalar product between operators. This is equivalent
to a linear space of matrices, converting the matrices effectively into
vectors $(\rho\rightarrow|\rho\rangle\rangle)$. This is the Fock-Liouville
space (FLS) [36]. The usual definition of the scalar product of matrices
$\phi$ and $\rho$ is defined as
$\langle\langle\phi\mid\rho\rangle\rangle\equiv\operatorname{Tr}\left[\phi^{\dagger}\rho\right]$.
The Liouville superoperator from Eq. (1) is now an operator acting on the
Hilbert space of density matrices. The main utility of the FLS is to allow the
matrix representation of the evolution operator. For example, the evolution of
a two-level system density matrix can be expressed in the FLS as

$|\rho\rangle\rangle=\left(\begin{array}[]{c}\rho_{00}\\\ \rho_{01}\\\
\rho_{10}\\\ \rho_{11}\end{array}\right).$ (40)

The Liouvolle von Neumann equation describes the time evolution of a mixed
state Eq. (1). In vector notation, the Liouvillian superoperator is expressed
as a matrix,

$\tilde{\bm{\mathcal{L}}}=\left(\begin{array}[]{cccc}0&i\Omega&-i\Omega&0\\\
i\Omega&iE&0&-i\Omega\\\ -i\Omega&0&-iE&i\Omega\\\
0&-i\Omega&i\Omega&0\end{array}\right),$ (41)

Each row is calculated just by observing the output of the operation
$-i[H,\rho]$ in the computational basis of the density matrices space. The
system’s time evolution corresponds to the matrix equation
$\left.\frac{d|\rho\rangle\rangle}{dt}=\tilde{\mathcal{L}}|\rho\rangle\right\rangle$,
which in matrix notation would be

$\left(\begin{array}[]{l}\dot{\rho}_{00}\\\ \dot{\rho}_{01}\\\
\dot{\rho}_{10}\\\
\dot{\rho}_{11}\end{array}\right)=\left(\begin{array}[]{cccc}0&i\Omega&-i\Omega&0\\\
i\Omega&iE&0&-i\Omega\\\ -i\Omega&0&-iE&i\Omega\\\
0&-i\Omega&i\Omega&0\end{array}\right)\left(\begin{array}[]{l}\rho_{00}\\\
\rho_{01}\\\ \rho_{10}\\\ \rho_{11}\end{array}\right).$ (42)

A similar approach is used for the dissipative part $\tilde{\bm{\mathcal{D}}}$

##  Appendix B Semi-global Propagation method

Obtaining high-fidelity results for controlling quantum gates requires highly
accurate and efficient numerical propagators to solve the Liouville von
Neumann equation. For this task, we have modified the semi-global propagator
[31] to perform in the Liouville vector space.

The Liouville von Neumann equation describes the dynamics of an open quantum
system:

$\frac{d}{dt}\hat{\bm{\rho}}={\cal L}\hat{\bm{\rho}}$ (43)

where $\hat{\bm{\rho}}$ is the density operator and ${\cal L}$ the
Liouvillian. For a driven open system, the Liuouvillian can be partitioned to:

$\displaystyle\frac{d}{dt}\hat{\bm{\rho}}={\cal
L}_{t}\hat{\bm{\rho}}=\left({\cal L}_{H}(t)+{\cal
L}_{D}(t)\right)\hat{\bm{\rho}}$ (44) $\displaystyle{\cal L}_{H}={\cal
L}_{H_{0}}+{\cal L}_{H_{t}}$

$\mathcal{L}_{H}(t)$ is the generator of the unitary part of the dynamics and
can be decomposed into time-independent and time-dependent components. The
dissipative part $\mathcal{L}_{D}(t)$ implicitly describes the influence of
the environment and can be time-dependent.

For a time-independent Lindblidian ${\cal L}_{0}$, a formal solution becomes:

$\hat{\bm{\rho}}{(t)}=e^{{\cal L}_{0}t}\hat{\bm{\rho}}(0)$ (45)

with the initial conditions $\hat{\bm{\rho}}(0)$.

When the Liouviilan can be partitioned into a time-dependent and time-
independent part ${\mathcal{L}}={\mathcal{L}}_{0}+{\mathcal{L}}_{t}$ a formal
solution of Eq. (44) can be written as an integral equation:

$\hat{\bm{\rho}}(t)=e^{{\cal L}_{0}t}\hat{\bm{\rho}}(0)+\int_{0}^{t}e^{{\cal
L}_{0}(t-\tau)}{\cal{L}}_{t}\hat{\bm{\rho}}(\tau)d\tau$ (46)

Eq. (46) will form the basis for the numerical approximation.

In typical control problems $\mathcal{L}$ varies considerably with the time.
Therefore, the total evolution is practically broken into finite time steps
$\Delta t$. Then, one can concatenate the propagators and obtain the total
evolution from $t=0$ to $t=T$ by

$\overrightarrow{\mathbf{\rho}}(T)\approx\prod_{j=1}^{N_{t}}{\mathcal{G}}_{j}(\Delta
t)\overrightarrow{\mathbf{\rho}}(0)$ (47)

where ${\mathcal{G}}_{j}(\Delta t)$ is the propagator for $t$ to $t+\Delta t$
where $t=j\Delta t$.

A direct approximation assumes that ${\mathcal{L}}_{t}$ is time independent
within a time step, then

${\mathcal{G}_{j}}\approx e^{{\mathcal{L}}_{j}\Delta t}$ (48)

where ${\mathcal{L}}_{j}={\mathcal{L}}(t+\Delta t/2)$. Sampling
${\mathcal{L}}$ in the middle of the time step leads to second-order accuracy
in $\Delta t$.

In what follows, we formulate the density operator as a vector
$\overrightarrow{\mathbf{\rho}}$ and the Liouvillian superoperator ${\cal L}$
as a matrix ( see above App. A [28]).

A numerical method to solve Eq.(48) is based on expanding the exponent or any
analytic function in a polynomial series in the matrix ${\cal L}_{j}$:

$\overrightarrow{\mathbf{\rho}}(t+\Delta
t)\approx\sum_{n=0}^{K-1}a_{n}(t+\Delta
t)P_{n}\left(\mathcal{L}_{j}\right)\overrightarrow{\mathbf{\rho}}(t)$ (49)

where $P_{n}(\mathcal{L}_{j})$ is a polynomial of degree $n$, and
$a_{n}(t+\Delta t)$ is the corresponding expansion coefficient in the interval
$t,t+\Delta t$. This requires choosing the set of expansion polynomials
$P_{n}(\mathcal{L}_{j})$ and the corresponding coefficients $a_{n}$ [37].

The expansion (49) has to be accurate in the eigenvalue domain of
$\mathcal{L}_{j}$ so that the form (49) will be helpful to represent
${\mathcal{G}_{j}}$. Successive matrix-vector multiplications can compute this
series of polynomials at Eq. (49). It scales as $O\left(N^{2}\right)$, and the
computational effort can be reduced even further. For sparse superoperators,
the matrix-vector operation can be replaced by an equivalent semi-linear
scaling, with $\sim O(N)$ [38].

An immediate question is how to choose the set of expansion polynomials
$P_{n}(x)$. As a thumb rule, we search for a polynomial basis that converges
the fastest. An orthogonal set of polynomials is the first step for fast
convergence [39].

An efficient implementation can be done recursively. The Chebyshev polynomial
and Newton interpolation polynomials are two orthogonal expansion series using
$P_{n}(\mathcal{L}_{j})$. When the Hamiltonian is non-Hermitian, the
eigenvalue domain becomes complex, and the Chebychev approach is not
appropriate anymore. Then, the Newtonian or Arnoldi approach should be used
instead [39, 40, 41].

Note that in Eq. (45), only the coefficients $a_{n}(t)$ are time-dependent.
The solution at intermediate time points can be obtained by calculating the
coefficients for intermediate points with negligible additional computational
effort.

Quantum control of gates requires exceptionally high accuracy. The convergence
rate of Eq. (47) with a piecewise constant ${\mathcal{L}}_{j}$ is slow,
leading to extensive numerical effort. To obtain faster convergence, we must
consider the problem of time ordering within the time step $\Delta t$.

To overcome the problem of time ordering, we will combine the polynomial
solution of Eq. (45) and the integral equation formal solution (46). In Eq.
(46), the free propagator appears both as a free term and at the integrand.
The solution of the integral equation requires an iterative approach since
$\hat{\bm{\rho}}(\tau)$ appears also in the integrand. The iteration is
carried out by extrapolating the solution from one time step to the next, $t$
to $t+dt$.

The integral in the formal solution Eq. (46) is reformulated as an
inhomogeneous source term:

$\frac{d\overrightarrow{\mathbf{\rho}}(t)}{dt}=\mathcal{L}_{j}\overrightarrow{\mathbf{\rho}}(t)+\overrightarrow{\mathbf{s}}(t)$
(50)

The source term will represent the time-dependent/nonlinear part of the
dynamics.

Treating Eq. (50) as an inhomogeneous ODE will give rise to a formal solution
to the time-dependent problem. We can write the solution to Eq. (45):

$\displaystyle\overrightarrow{\mathbf{\rho}}(t+\Delta t)$
$\displaystyle=\tilde{\bm{\mathcal{G}}}_{j}(t)\overrightarrow{\mathbf{\rho}}(t)+\int_{t}^{t+\Delta
t}\tilde{\bm{\mathcal{G}}}_{j}(t-\tau)\overrightarrow{\mathbf{s}}(\tau)d\tau$
(51)
$\displaystyle=\exp\left(\mathcal{L}_{j}t\right)\overrightarrow{\mathbf{\rho}}+\int_{t}^{t+\Delta
t}\exp\left[\mathcal{L}_{j}(t-\tau)\right]\overrightarrow{\mathbf{s}}(\tau)d\tau$
$\displaystyle=\exp\left(\mathcal{L}_{j}t\right)\overrightarrow{\mathbf{\rho}}_{0}+\exp\left(\mathcal{L}_{j}t\right)\int_{t}^{t+\Delta
t}\exp\left(-\mathcal{L}_{j}\tau\right)\overrightarrow{\mathbf{s}}(\tau)d\tau$

Where $\tilde{\bm{\mathcal{G}}}_{j}$ is defined as the time-independent
propagator by the vec-ing procedure: subsection II.2.

The source term is expanded as a time-dependent polynomial to solve for the
integral analytically.

$\overrightarrow{\mathbf{s}}(t)=\sum_{n=0}^{M-1}\frac{t^{n}}{n!}\overrightarrow{\mathbf{s}}_{n}$
(52)

This expansion allows us to solve formally the integral in Eq. (51)
$\int e^{at}t^{m}/m!=\sum_{n=1}^{m}e^{at}t^{n-m}/a^{n}(n-m)!$.

The problem is now shifted to obtaining the expansion coefficients
$\overrightarrow{\mathbf{s}}_{m}$. This task is obtained by approximating
$\overrightarrow{\bm{s}}(t)$ by an orthogonal polynomial in the time interval.
We chose a Chebychev expansion

$\overrightarrow{\mathbf{s}}(t)\approx\sum_{n=0}^{M-1}\overrightarrow{\mathbf{b}}_{n}T_{n}(t)$
(53)

where the coefficients $\overrightarrow{\mathbf{b}}_{n}$ are computed by a
scalar product of the $T_{n}(t)$ with $\overrightarrow{\mathbf{s}}(t)$.
Approximating the coefficients using Chebychev sampling points in the time
interval $\Delta t$.

The coefficients $\overrightarrow{\mathbf{s}}_{m}$ are calculated by relating
the polynomial Eq. (53) to the Chebychev expansion. This source term is
inserted into the integral Eq. (51), leading to a numerical approximation to
the solution of the TDLE. The addition of the source term into the dynamics
gives rise to an analytical solution for the last term in Eq. (51), presented
here on the RHS of Eq. (54)

$J_{m+1}(\mathcal{L}_{j},t)\equiv\int_{t}^{t+\Delta
t}\exp(-\mathcal{L}_{j}\tau)\tau^{m}d\tau,\quad m=0,1,\ldots$ (54)

with the recursion relations:

$J_{m}(\mathcal{L}_{j},t)=-\frac{\exp(-\mathcal{L}_{j}t)t^{m-1}}{\mathcal{L}_{j}}+\frac{m-1}{\mathcal{L}_{j}}J_{m-1}(\mathcal{L}_{j},t),\quad
m=2,3,\ldots$ (55)

where

$\displaystyle J_{1}(\mathcal{L}_{j},t)\equiv\int_{t}^{t+\Delta
t}\exp(-\mathcal{L}_{j}\tau)d\tau=\frac{1-\exp(-\mathcal{L}_{j}t)}{\mathcal{L}_{j}}.$

Plugging Eq. (52) into this formulation leads to the following:

$\exp(\mathcal{L}_{j},t)\sum_{n=0}^{M-1}\frac{1}{n!}\int_{0}^{t}\exp(-\mathcal{L}_{j}\tau)t^{n}d\tau
s_{n}=\exp(\mathcal{L}_{j}t)\sum_{n=0}^{M-1}\frac{1}{n!}J_{n+1}(\mathcal{L}_{j},t)s_{n}=\sum_{n=0}^{M-1}f_{n+1}(\mathcal{L}_{j},t)s_{n}$
(56)

In Eq. (44), the Louivilian was split into explicit time-dependent and
approximated time-independent parts; this partition is used to re-write the
TDLE

$\frac{d\overrightarrow{\mathbf{\rho}}(t)}{dt}={\mathcal{L}}_{j}\overrightarrow{\mathbf{\rho}}(t)+\overrightarrow{\mathbf{s}}(\overrightarrow{\mathbf{\rho}}(t),t)$
(57)

Which is a reformulation of Eq. (50). The same analysis leads to:

$\overrightarrow{\mathbf{\rho}}(t,t+\Delta
t)=\exp(\mathcal{L}_{j}t)\overrightarrow{\mathbf{\rho}}(t)+\exp({\mathcal{L}_{j}}t)\int_{t}^{t+\Delta
t}\exp(-{\mathcal{L}_{j}}\tau)\overrightarrow{\mathbf{s}}(\overrightarrow{\mathbf{\rho}}(\tau),\tau)d\tau$
(58)

Now, we can use these formulations to approximate Eq. (58):

$\overrightarrow{\mathbf{\rho}}(t,t+\Delta t)\approx
P_{M}({\mathcal{L}_{j}},(t,t+\Delta t))\overrightarrow{S}(t,t+\Delta
t)_{M}+\sum_{n=0}^{M-1}\frac{t^{n}}{n!}\overrightarrow{S}(t,t+\Delta t)_{n}$
(59)

$P_{M}({\mathcal{L}_{j}},(t,t+\Delta t)\overrightarrow{S}(t,t+\Delta t)_{M}$,
is approximated by the Arnoldi method (the eigenvalue spectrum of
$\mathcal{L}$ is distributed on the complex plane) where

$\displaystyle\overrightarrow{{S(t)}}_{M}\equiv\overrightarrow{\mathbf{s}}(t)+\mathcal{L}_{t}\overrightarrow{\mathbf{\rho}(t)}$

are computed by expanding it by time in the same form of (52). We have used
here the fact that
$P_{n}(\mathcal{L_{H}},t)=\mathcal{L_{H}}^{k-n}P_{k}(\mathcal{L_{H}},t)+\sum_{j=n}^{k-1}\frac{t^{j}}{j!}\mathcal{L_{H}}^{j-n}$
Eq. (59) and $\overrightarrow{\mathbf{s}}_{j}$ includes dependence on
$\overrightarrow{\mathbf{\rho}}(t)$. It would seem we are back to the same
problem. Nevertheless, it can be overcome by an iterative procedure. First, we
guess a solution $\overrightarrow{\mathbf{\rho}}_{g}(t)$, within a time step
$\Delta t$ and use it in Eq. (59) to obtain a new approximate solution. This
procedure can be continued until the solution converges with the desired
accuracy. The initial guess is extrapolated from the previous time step to
accelerate convergence.

Three numerical parameters determine the precision of the propagation and the
convergence rate:

  * •

The size of the time step $\Delta t$.

  * •

Number of Chebyshev sampling points in each time step $M$.

  * •

The size of the Krylov space $K$ corresponds to the basis of the Arnoldi
algorithm. It is important to note that this parameter is limited by
$Dim\\{{\mathcal{L}\\}-1}$

Each of those is adjustable by the user to fit their needs best. For example,
in the Hadamard transformation system, we have used the following parameters:
$\Delta t=0.1$,$M=7$, and $K=3$. With these parameters, we got an accuracy of
$10^{-8}$ for the propagator, two orders of magnitude higher than the fidelity
of the target transformation. For the entangling gate, we adjusted the
parameters for a higher resolution that would fit $10^{-8}$ for the fidelity
of the target transformation ($M=K=9$).

## References

  * Preskill [2018] J. Preskill, Quantum 2, 79 (2018).
  * Schlosshauer [2019] M. Schlosshauer, Physics Reports 831, 1 (2019).
  * Kraus [1971] K. Kraus, Annals of Physics 64, 311 (1971).
  * Gorini and Kossakowski [1976] V. Gorini and A. Kossakowski, Journal of Mathematical Physics 17, 1298 (1976).
  * Lindblad [1976] G. Lindblad, Communications in Mathematical Physics 48, 119 (1976).
  * Bartana  _et al._ [1997] A. Bartana, R. Kosloff, and D. J. Tannor, The Journal of chemical physics 106, 1435 (1997).
  * Ohtsuki  _et al._ [1999] Y. Ohtsuki, W. Zhu, and H. Rabitz, The Journal of chemical physics 110, 9825 (1999).
  * Bartana  _et al._ [2001] A. Bartana, R. Kosloff, and D. J. Tannor, Chemical Physics 267, 195 (2001).
  * Schulte-Herbrüggen  _et al._ [2011] T. Schulte-Herbrüggen, A. Spörl, N. Khaneja, and S. Glaser, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 154013 (2011).
  * Koch [2016] C. P. Koch, Journal of Physics: Condensed Matter 28, 213001 (2016).
  * Pokharel  _et al._ [2018] B. Pokharel, N. Anand, B. Fortman, and D. A. Lidar, Physical review letters 121, 220502 (2018).
  * Petruhanov and Pechen [2023a] V. Petruhanov and A. Pechen, Journal of Physics A: Mathematical and Theoretical 56, 305303 (2023a).
  * Delben  _et al._ [2023] G. J. Delben, M. W. Beims, and M. G. E. da Luz, Phys. Rev. A 108, 012620 (2023).
  * Petruhanov and Pechen [2023b] V. Petruhanov and A. Pechen, arXiv preprint arXiv:2309.02063 (2023b).
  * Dann and Kosloff [2021] R. Dann and R. Kosloff, Quantum 5, 590 (2021).
  * Kallush  _et al._ [2022] S. Kallush, R. Dann, and R. Kosloff, Science Advances 8, eadd0828 (2022).
  * Viola  _et al._ [1999] L. Viola, E. Knill, and S. Lloyd, Physical Review Letters 82, 2417 (1999).
  * Khasin and Kosloff [2011] M. Khasin and R. Kosloff, Physical review letters 106, 123002 (2011).
  * Kallush  _et al._ [2014] S. Kallush, M. Khasin, and R. Kosloff, New Journal of Physics 16, 015008 (2014).
  * Jiang  _et al._ [2023] X. Jiang, J. Scott, M. Friesen, and M. Saffman, Physical Review A 107, 042611 (2023).
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976).
  * Feldmann and Kosloff [2010] T. Feldmann and R. Kosloff, Europhysics Letters 89, 20004 (2010).
  * Kiely [2021] A. Kiely, Europhysics Letters 134, 10001 (2021).
  * Kosloff and Feldmann [2010] R. Kosloff and T. Feldmann, Physical Review E 82, 011134 (2010).
  * Xuereb  _et al._ [2023] J. Xuereb, P. Erker, F. Meier, M. T. Mitchison, and M. Huber, arXiv preprint arXiv:2301.10767 (2023).
  * Agarwal [1976] G. S. Agarwal, Phys. Rev. Lett. 37, 1383 (1976).
  * Am-Shallem  _et al._ [2015] M. Am-Shallem, A. Levy, I. Schaefer, and R. Kosloff, arXiv preprint arXiv:1510.08634 (2015).
  * Palao and Kosloff [2003] J. P. Palao and R. Kosloff, Physical Review A 68, 062308 (2003).
  * Krotov [1995] V. Krotov,  _Global methods in optimal control theory_ , Vol. 195 (CRC Press, 1995).
  * Schaefer  _et al._ [2017] I. Schaefer, H. Tal-Ezer, and R. Kosloff, Journal of Computational Physics 343, 368 (2017).
  * Ramakrishna and Rabitz [1996] V. Ramakrishna and H. Rabitz, Physical Review A 54, 1715 (1996).
  * Levy  _et al._ [2018] A. Levy, A. Kiely, J. G. Muga, R. Kosloff, and E. Torrontegui, New Journal of Physics 20, 025006 (2018).
  * Boixo  _et al._ [2007] S. Boixo, L. Viola, and G. Ortiz, Europhysics Letters 79, 40003 (2007).
  * Wu  _et al._ [2007] R. Wu, A. Pechen, C. Brif, and H. Rabitz, Journal of Physics A: Mathematical and Theoretical 40, 5681 (2007).
  * Manzano [2020] D. Manzano, Aip Advances 10 (2020).
  * Berman and Kosloff [1991] M. Berman and R. Kosloff, Computer physics communications 63, 1 (1991).
  * Kosloff and Kosloff [1983] D. Kosloff and R. Kosloff, Journal of Computational Physics 52, 35 (1983).
  * Kosloff [1994] R. Kosloff, Annual review of physical chemistry 45, 145 (1994).
  * Arnoldi [1951] W. E. Arnoldi, Quarterly of applied mathematics 9, 17 (1951).
  * Lehoucq and Sorensen [1996] R. B. Lehoucq and D. C. Sorensen, SIAM Journal on Matrix Analysis and Applications 17, 789 (1996).