1 Introduction↩︎

Quantum optimal control (QOC) has emerged as a pivotal technique for manipulating quantum systems efficiently and precisely, essential in quantum computing, quantum sensing, and quantum communication. The primary objective is to design control pulses that drive a quantum system from an initial state to a target state or implement a specific quantum operation with high fidelity, often under stringent constraints on resources such as energy and time.

Two prominent methods in QOC are Gradient Ascent Pulse Engineering (GRAPE)[1] and Pontryagin’s Maximum Principle (PMP)[2]. GRAPE optimizes discretized control pulses by iteratively adjusting them in the direction of increasing fidelity. Conversely, PMP offers an analytical characterization of optimal solutions but often lacks practical computational implementation, especially for complex quantum systems.

In this work, we introduce MAGICARP (Method for Adjoint- and Gradient-based self-Iterative Construction And Refinement of Pulses), a hybrid method leveraging both the analytical insights from PMP and the numerical efficiency of gradient ascent algorithms. MAGICARP functions as a shooting method in control theory, determining the correct initial adjoint momentum (or costate) required to achieve a desired final state or gate operation. This initial momentum guides the generation of control pulses consistent with the PMP’s necessary conditions for optimality, explicitly accommodating either time-optimal or energy-optimal criteria. The time required to implement a gate is a crucial factor, it impacts the coherence time of the system, and the energetic cost of the control Hamiltonians is also an important factor to consider, as it impacts the energy dissipated in the system. Moreover, in the presence of control and environmental errors, the robustness of the gate to these errors is also a key factor to consider.

In the context of this work we consider two general types of objective or cost functionals [3], [4], known as problems of

(i) Mayer type:

\[\label{eq:MayerType} \mathcal{J}[u] = \phi(\rho(T), T),\tag{1}\]

and (ii) Lagrange type:

\[\label{eq:LagrangeType} \mathcal{J}[u] = \int_0^T \mathcal{L}[\rho(t), u(t), t] \, \mathrm{d}t.\tag{2}\]

As shown in 1 , Mayer problems focus on quantities dependent solely on the final state at time \(T\), such as maximizing the fidelity between the evolved state \(\rho(T)\) and a desired target state. Lagrange problems, represented in 2 , account for costs accumulated over the entire time interval \([0, T]\), including factors like the total energy used by control fields \(u(t)\) or penalties for deviations during the system’s evolution. One often considers hybrid cost functionals that combine both Mayer and Lagrange terms, known as problems of Bolza type:

\[\label{eq:BolzaType} \mathcal{J}[u] = \alpha \phi(\rho(T), T) + (1 - \alpha) \int_0^T \mathcal{L}[\rho(t), u(t), t] \, \mathrm{d}t,\tag{3}\]

where \(\alpha \in [0, 1]\) is a weighting parameter that balances the optimization priorities. By varying \(\alpha\), one can smoothly transition between focusing on the final-state fidelity and minimizing resource costs. For instance, starting with \(\alpha = 1\) emphasizes reaching the target state with maximum fidelity. Gradually decreasing \(\alpha\) introduces considerations like energy efficiency, allowing the optimization to find control fields that are both effective and resource-conscious [5].

The choice of a good cost functional is crucial for the success of the optimization process, in particular if implemented numerically and it is not a global optimization, the wrong choice of cost functional can lead to suboptimal solutions due to the optimization algorithm getting stuck in local minima.

1.1 Pontryagin’s Maximum Principle↩︎

There are numerous numerical optimization algorithms that can be used to solve the optimal control problem, such as the GRAPE algorithm, the Krotov method, or the CRAB algorithm, among others[6], [7]. However, these algorithms are often computationally expensive and can be sensitive to the choice of initial conditions and parameters. Moreover, they may not always provide the global optimum of the cost functional, and they may not be able to handle constraints on the control fields.

The cost problem can also be tackled analytically, and, while it is probably not solvable analytically, it can give us insights that could then be used to design better numerical optimization algorithms, with less free parameters for example.

The optimal control problem is a well-studied field in control theory, and one of the most powerful tools to solve it is Pontryagin’s Maximum Principle (PMP). The PMP is a necessary condition for optimality in optimal control problems, and it provides a set of differential equations, known as the Pontryagin equations, that the optimal control fields must satisfy. The PMP is a generalization of the Euler-Lagrange equations in classical mechanics, and it is a cornerstone of optimal control theory.

Let us consider a general optimal control problem, where we aim to minimize a cost functional of the form 2 :

\[\mathcal{J}[x, \vec{u}] = \int_0^T \mathcal{L}[x(t), \vec{u}(t), t] \, \mathrm{d}t,\]

where \(x(t)\) is the state of the system at time \(t\), \(\vec{u}(t)\) are the control fields, and \(\mathcal{L}[x(t), \vec{u}(t), t]\) is then the Lagrangian of the system.

Assuming the system is governed by a set of differential equations, known as the state equations, that describe the evolution of the system in the absence of control fields:

\[\dot{x}(t) = f[x(t), t],\]

Note: If \(\lambda^0\) is not zero, then it is possible to rescale the equations by defining \(\lambda \leftarrow \lambda/|\lambda^0|\), and \(\mathcal{H} \rightarrow \mathcal{H}/|\lambda^0|\), such that, effectively \(\lambda^0=-1\).

Note: When one applies the PMP where \(x\) are the center of mass coordinates of an object in classical mechanics and the Langrangian is the traditional kinetic energy minus the potential energy, the PMP gives the equations of motion of the object, with \(\lambda\) being the momentum \(p\) of the object, and the pseudo-Hamiltonian is then the real Hamiltonian, i.e. total energy of the object. Therefore, the PMP is a generalization of the Hamilton equations in classical mechanics.

In our particular case, the state is not the density matrix (i.e. the litteral state of the system), but the unitary operator \(U(t)\), and the control fields are \(u_k(t)\) in the control Hamiltonian

\[H_c(t) = \sum_k u_k(t) H_k(t).\]

The previous set of equations can be written as:

\[\begin{align} f[U] &= \dot{U} = -i H U, \tag{4}\\ \mathcal{H}[U, \lambda, \vec{u}, t] &= \lambda^\dagger(t) \cdot f[U] - \mathcal{L}[U,\vec{u},t], \tag{5}\\ \dot{\lambda} &= -i H \lambda + \left(\frac{\partial \mathcal{L}}{\partial U}\right)^\dagger, \tag{6}\\ \vec{0} &= \left.\frac{\partial \mathcal{H}}{\partial \vec{u}}\right\vert_{u=u^*}, \tag{7}\\ \mathcal{H}[U^*(t), \lambda^*(t), \vec{u}^*(t), t] &= C \geq 0,\quad\forall t\in [0,T].\tag{8} \end{align}\]

where \(\lambda\) is the adjoint matrix, and \(\mathcal{L}\) is the Lagrangian to be determined. The scalar product for two matrices \(A\) and \(B\) is defined as \(A\cdot B = \operatorname{ReTr}(A^\dagger B)\).

1.1.1 Multiple control fields for a qudit in the interaction picture↩︎

We apply the Pontryagin Maximum Principle (PMP) to the problem of identifying optimal control fields for implementing a quantum gate on a \(d\)‑level qudit. Working in the interaction picture, the free evolution is removed and the Hamiltonian is \[H(t)=\sum_k u_k(t)H_k,\] where \(H_k\) are control Hamiltonians and \(u_k(t)\) their corresponding control fields. Our goal is to realise a target unitary \(U_{\text{target}}\) at a normalised final time \(T=1\).

It will be convienient to define the cumulative envelope of the control fields as \[c(t)\equiv\sqrt{\sum_k u_k^2(t)}.\]

The choice of the final time \(T\) is arbitrary, as, for a driven quantum system in the interaction picture, for any continuous integrable scalar function \(S(t)\), the following time-dependent rescaling of the control fields

\[\Bigl\{u_k(t)\Bigr\}_{k=1}^m \;\longrightarrow\; \Bigl\{\tfrac{u_k(t)}{S(t)}\Bigr\}_{k=1}^m,\]

accompanied by the time-reparametrisation

\[t'=\!\!\int_{0}^{t}\!|S(s)|\,ds, \qquad dt'=|S(t)|\,dt, \qquad \tilde{T}=\int_0^T |S(s)|\,ds,\]

leaves the overall evolution operator \(\,U=\mathcal{T}\,\exp\!\bigl[-\tfrac{i}{\hbar}\!\int_0^{T}\!H(t)\,dt\bigr]=\mathcal{T}\,\exp\!\bigl[-\tfrac{i}{\hbar}\!\int_0^{\tilde{T}}\!H(t')\,dt'\bigr]\) unchanged, i.e. the total pulse area is conserved; only the clock that measures it is stretched.

This means that for a given set of pulses applied over a time interval \([0,1]\), there are easily computable rescaled pulses \(u_k(t')=\Omega_\text{max}u_k(t)/c(t)\) that yield the same evolution operator \(U(1)\), but with a different duration

\[\widetilde{T}=\int_0^1 c(t)/\Omega_\text{max}\,\mathrm dt=\frac{1}{\Omega_\text{max}}\int_0^1 \sqrt{\sum_k u_k^2(t)}\,\mathrm dt, \label{eq::T95max}\tag{9}\]

Moreover, this rescaling assures that the cumulative envelope \(c(t)\) is normalised to a given maximum attainable amplitude constraint on the control fields, \(\Omega_\text{max}\), the duration of the gate in units of \(\Omega_{\max}^{-1}\) is then simply the integral of the cumulative envelope over the time interval \([0,1]\).

1.1.1.1 Simple quadratic cost functional

We first treat the elementary quadratic cost \[\mathcal{J}[\vec{u}]=\int_0^1\sum_k u_k^2(t)\,\mathrm dt = \int_0^1 c^2(t)\,\mathrm dt, \label{eq::costsimple}\tag{10}\] which penalises the integrated pulse energy. PMP gives \[\begin{align} \mathcal{H}[U(t), \vec{\lambda}(t), \vec{u}(t)] &= \operatorname{ImTr}\left(\lambda^\dagger(t) H(t) U(t)\right) - \sum_k u_k^2(t), \tag{11}\\ \dot{\lambda}(t) &= -i H(t) \lambda(t), \tag{12}\\ \vec{0} &= \operatorname{ImTr}\left(\lambda^\dagger(t) \frac{\partial H(t)}{\partial \vec{u}^*} U(t)\right) - 2 \vec{u}^*(t),\tag{13}\\ \mathcal{H}[U^*(t), \vec{\lambda}^*(t), \vec{u}^*(t)] &= C \equiv c^2 \geq 0,\quad\forall t\in [0,1].\tag{14} \end{align}\]

Let us now analyze the adjoint equation 12 . Since \(U(t) \in SU(d)\), the adjoint state \(\lambda(t)\) belongs to the tangent space at \(U(t)\) and can be written as:

\[\lambda(t) = U(t) \lambda_0,\]

where \(\lambda_0 \in \mathfrak{su}(d)\) is a constant, skew-Hermitian matrix. Inserting into the adjoint equation, we obtain:

\[\dot{\lambda}_0 = 0 \quad \Rightarrow \quad \lambda(t) = U(t) \lambda_0. \label{eq::adjoint95matrix950}\tag{15}\]

Using this in 13 and the identity \(\frac{\partial H}{\partial u_k} = H_k\), we find that for all \(k\):

\[\operatorname{ImTr}\left( U(t) \lambda_0 U^\dagger(t) H_k \right) = 2 u_k^*(t). \label{eq::optimal95control95l0952}\tag{16}\]

Introducing \(g = i \lambda_0\), we find that \(g\) is a traceless Hermitian matrix, and the optimal control fields satisfy:

\[u_k^*(t) = \frac{1}{2} \operatorname{ReTr}\left(U(t) g U^\dagger(t) H_k \right) = \frac{1}{2} g(t) \cdot H_k = \frac{1}{2} g \cdot H_k(-t), \label{eq::optimal95control2}\tag{17}\]

with \(g(t) = U(t) g U^\dagger(t)\) the time-evolved adjoint, or equivalently \(H_k(-t) = U^\dagger(t) H_k U(t)\) the time-reversed Hamiltonians.

\(u_k^*(t)\) is uniquely defined by \(g(t)\), the latter being a solution of the first-order differential equation:

\[\dot{g}(t) = -i[H(t), g(t)],\]

implies \(g(t)\) is uniquely determined by its initial condition \(g(0) = g = i \lambda_0\).

Any set of control fields satisfying 17 at all times is optimal in the sense of minimizing the cost in 10 , and the pulses – and a fortiori the final gate \(U(1)\) – can be constructed a priori, by choosing an initial adjoint momentum \(\lambda_0\) or \(g\) and then solving the differential equation for \(g(t)\), which in turn gives the optimal control fields \(u_k^*(t)\).

Note: The trivial solution \(u_k(t) = 0\) corresponds to the identity gate and arises when \(g \propto \mathbb{1}_d\), noting that \(H_k\) are traceless.

Moreover, 14 leads to useful insights about the envelope of the optimal control fields, since, by recalling that \(H(t) = \sum_k u_k(t) H_k\), it gives us:

\[\sum_k u_k^*(t) \operatorname{ImTr}\left(U(t) \lambda_0 U^\dagger(t)H_k\right) - \sum_k (u_k^*(t))^2 = c^2,\]

which, using 16 , simplifies to:

\[\sum_k u_k^2(t) = c^2,\]

where \(c\) is a constant. This means that the cumulative envelope of the control fields is constant over the time interval \([0,1]\), and this implies that for any given maximal driving amplitude \(\Omega_{\max}\), as per 9 , the duration of the pulse is fixed to \[\widetilde{T} = \frac{1}{\Omega_{\max}}\int_0^1 \sqrt{\sum_k u_k^2(t)} \, dt = |c|/\Omega_{\max}. \label{eq::T95max2}\tag{18}\]

Meanwhile, the energy of the pulses is simply

\[\widetilde{E} = \int_0^{\widetilde{T}} \Omega_{\max}^2 \, dt = |c|\Omega_{\max}. \label{eq::E95max}\tag{19}\]

Conversely, if the time duration \(\widetilde{T}\) is fixed, the energetic cost of the pulses is

\[\widetilde{E} = c^2 / \widetilde{T} \label{eq::E95max2}\tag{20}\]

In both cases, going from the time interval \([0,1]\) to the physical time interval \([0,\widetilde{T}]\) is done by rescaling the time by a factor \(c/\Omega_{\max} = \widetilde{T}\) and the control fields by a factor \(\Omega_{\max}/c=\widetilde{T}^{-1}\). The latter can be absorbed into \(g\) in 17 (\(\overline{g}_E = g\Omega_{\max}/c = g/T\)), which allows to apply the definition of the optimal controls to any time interval \([0,T]\) without loss of generality.

This means that if one wants to implement a gate with a fixed duration \(T\), but a minimal energetical cost, a procedure is to find a constant traceless Hermitian matrix \(\overline{g}\) such that \(U(t)\) is the desired gate at time \(t=T\) – up to a global phase and up to a tolerance on the fidelity of the gate. The energetic cost of the gate, which is simply the integral of the squared cumulative envelope over the time interval \([0,T]\), can be computed from \(\overline{g}_E\) as \[\widetilde{E} = \int_0^T \sum_k u_k^2(t)\, dt = \sum_k (u_k^*(0))^2 T = \sum_k \left(\frac{1}{2} \operatorname{ReTr}\left(\overline{g}_E H_k\right)\right)^2 T.\]

In the case where the amplitude \(\Omega_{\max}\) is fixed, finding a constant traceless Hermitian matrix \(g\) such that \(U(t)\) is the desired gate at time \(t=1\) allows to compute the time necessary to implement the gate while trying to minimize the energy, c.f. 18 , i.e. \[\widetilde{T} = \frac{1}{2\Omega_{\max}} \sqrt{\sum_k \left( \operatorname{ReTr}\left(g H_k\right) \right)^2}.\]

Note: The PMP does not provide a constructive algorithm to find the optimal control fields, as the PMP only gives necessary conditions for optimality, not sufficient ones. However, it does provide a structure that can be used to either construct good candidate pulses, or that cen be incorporated into a numerical optimization algorithm, to improve its convergence toward the optimal solution.

1.1.1.2 Time‑optimal cost functional

To minimise the gate duration itself, we refer to 9 , and we take \[\mathcal{J}[\vec{u}]=\widetilde{T}=\int_0^1\sqrt{\sum_k u_k^2(t)}\,\mathrm dt = \int_0^1 c(t)\,\mathrm dt, \label{eq::Lagrangian}\tag{21}\] PMP now yields

\[\begin{align} \mathcal{H}[U(t), \vec{\lambda}(t), \vec{u}(t)] &= \operatorname{ImTr}\left(\lambda^\dagger(t) H(t) U(t)\right) - \sqrt{\sum_k u_k^2(t)}, \tag{22}\\ \dot{\lambda}(t) &= -i H(t) \lambda(t), \tag{23}\\ \vec{0} &= \operatorname{ImTr}\left(\lambda^\dagger(t) \frac{\partial H(t)}{\partial \vec{u}^*} U(t)\right) - \vec{u}^*(t)/c(t),\tag{24}\\ \mathcal{H}[U^*(t), \vec{\lambda}^*(t), \vec{u}^*(t)] &= C \geq 0,\quad\forall t\in [0,1].\tag{25} \end{align}\] With \(\lambda(t)=U(t)\lambda_0=-iU(t)g\), similarily to the previous cost function under interest, the optimal control fields satisfy \[u_k^*(t)=\operatorname{ReTr}\bigl(U(t)gU^\dagger(t)H_k\bigr), \label{eq::optimal95control}\tag{26}\] where \(g\) is a constant traceless Hermitian matrix, and \(c(t)=\sqrt{\sum_k u_k^2(t)}\) is the cumulative envelope of the control fields.

Substituting \(H(t)=\sum_k u_kH_k\) and (26 ) into \(\mathcal{H}=C\) (25 ) gives \[\sum_k u_k^*(t)\operatorname{ReTr}\bigl(U(t)gU^\dagger(t)H_k\bigr)-\sqrt{\sum_k u_k^{*2}(t)}=C.\] This leads to the following condition on the cumulative envelope of the control fields: \[c(t)(c(t)-1) = C,\] which, assuming c(t) is cotinuous over the time interval \([0,1]\), implies that the cumulative envelope \(c(t)\) is also constant over the time interval \([0,1]\).

Note: An intersting consequence, is that for a given maximal amplitude \(\Omega_{\max}\), and with an envelope \(c(t)\) fixed to the constant \(\Omega_{\max}\) over the time interval \([0,\widetilde{T}]\), the structure of the time-optimal control fields is the same as the energy-optimal pulse in ?? , with \(\Omega_{\max} \overline{g}_T = \frac{1}{2}\overline{g}_E\).

Combining the two necessary conditions for both the time-optimal and energy-optimal control fields, we find that an time-optimal family of control fields \(\{u_k(t)\}\) is also energy-optimal over the time interval \([0,\widetilde{T}]\).

1.1.1.3 Open questions

Although the PMP reveals the structure of optimal fields, it does not yet provide

1. a constructive computational algorithm, nor

2. a mechanism to incorporate gate fidelity.

These topics are addressed in the next section.

1.1.2 Single control field for a qudit in the laboratory frame↩︎

Let us now consider the case of a single control field \(u(t)\) acting on a qudit system, in the laboratory frame however, where the free evolution is not removed. The system Hamiltonian reads:

\[H(t) = H_0 + u(t) H_\text{c},\]

The free evolution is given by the drift Hamiltonian \(H_0\), and the control Hamiltonian is \(H_\text{c}\). The goal is to implement a target unitary operator \(U_\text{target}\) at time \(T\).

The presence of the drift Hamiltonian \(H_0\) does not change the form of the PMP equations 11 12 13 14 , nor does it change the form of the adjoint equation 15 , however the evolution operator \(U(t)\) incorporates now the drift Hamiltonian.

?? as a necessary condition for energy-optimality is therefore still valid with the drift Hamiltonian included in the evolution operator \(U(t)\), and the optimal control fields are given by: \[u^*(t) = {\Omega_\text{max}} \operatorname{ReTr}\left(U(t) \overline{g}_E U^\dagger(t) H_\text{c} \right),\]

where \(\overline{g}_E\) is a constant, traceless, Hermitian matrix.

The first major difference with the interaction picture is that rescaling the duration of the gate to \(T=1\) does not imply that the duration of the maximally driven control field is equal to its energetic cost as in 9 in units of \(\Omega_\text{max}^{-1}\). Therefore, a pulse of the aforementioned would only be energy-optimal for a fixed duration e.g. \(T=1\), and not time-optimal. However, the case of energy-optimal control fields is still interesting, and the last PMP equation

\[\mathcal{H}[U^*(t), \vec{\lambda}^*(t), u^*(t)] = C \geq 0,\quad\forall t\in [0,1],\]

allows us to determine that the envelope of the control field in the presence of a drift Hamiltonian is given by:

\[c(t) = |u^*(t)| = \sqrt{C-\operatorname{ReTr}(U(t) \overline{g}_E U^\dagger(t) H_0)}.\]

This leads to the following constrained differential equation for \(\overline{g}_E(t) = U(t) \overline{g}_E U^\dagger(t)\) to satisfy:

\[\begin{cases} \dot{\overline{g}_E }(t) = -i[H_0 + \operatorname{ReTr}(\overline{g}_E (t) H_\text{c})H_\text{c}, \overline{g}(t)], \\ \left(C - \operatorname{ReTr}(\overline{g}_E (t) H_0)\right)^2 = \operatorname{ReTr}(\overline{g}_E (t) H_\text{c})^2. \end{cases}\]

The case of time-optimal control fields with a drift Hamiltonian is more complex, as no rescaling of the time interval is possible, the time cost cannot be expressed in terms of the cumulative envelope, the correct cost functional is actually given by the control-independent time integral:

\[\mathcal{J}[T]=\int_0^T \mathrm dt=T, \label{eq::cost95time95optimal}\tag{27}\]

this means that the Langrangian \(\mathcal{L}\) is simply \(1\), moreover since \(T\) is not fixed anymore, the PMP equations are given by Boscain2021?:

\[\begin{align} \mathcal{H}[U(t), \lambda(t), u(t)] &= \operatorname{ImTr}\left(\lambda^\dagger(t) H(t) U(t)\right) - 1, \tag{28}\\ \dot{\lambda}(t) &= -i H(t) \lambda(t), \tag{29}\\ 0 &= \operatorname{ImTr}\left(\lambda^\dagger(t) \frac{\partial H(t)}{\partial u^*} U(t)\right),\tag{30}\\ \mathcal{H}[U^*(t), \lambda^*(t), u^*(t)] &= 0,\quad\forall t\in [0,T].\tag{31} \end{align}\]

The adjoint equation is unchanged, and we can rewrite the transversality condition 31 as: \[\operatorname{ReTr}\left(U(t) g U^\dagger(t) H_0\right) + u^*(t) \operatorname{ReTr}\left(U(t) g U^\dagger(t) H_\text{c}\right) = 1,\]

where \(g = i \lambda_0\) is a constant, traceless, Hermitian matrix. We find that the time-optimal control fields satisfy:

\[u^*(t) = \frac{1 - \operatorname{ReTr}\left(U(t) g U^\dagger(t) H_0\right)}{\operatorname{ReTr}\left(U(t) g U^\dagger(t) H_\text{c}\right)}.\]

Or, equivalently, by defining \(\Phi(t) = \operatorname{ReTr}\left(U(t) g U^\dagger(t) H_\text{c}\right)\), we have \[u^*(t) = \frac{1 - \operatorname{ReTr}\left(U(t) g U^\dagger(t) H_0\right)}{\Phi(t)} = \frac{1 - \operatorname{ReTr}\left(U(t) g U^\dagger(t) H_0\right)}{\|\Phi(t)\|} \frac{\Phi(t)}{\|\Phi(t)\|},\]

If constraints are imposed on the maximal amplitude of the control field, the optimal control field is bang-bang, i.e. it switches between -\(\Omega_{\max}\) and \(\Omega_{\max}\) Boscain2021?, except on singular arcs where \(\operatorname{ReTr}\left(U(t) g U^\dagger(t) H_\text{c}\right)=0\), which do generally not persist, and

\[\left|\frac{1 - \operatorname{ReTr}\left(U(t) g U^\dagger(t) H_0\right)}{\operatorname{ReTr}\left(U(t) g U^\dagger(t) H_\text{c}\right)}\right| = \Omega_{\max}.\]

Moreover if \([H_0, H_\text{c}] = 0\), \(\dot{\Phi}(t) = 0\), and the control field is constant, which is a trivial case.

However, in the unbounded case, the only physically acceptable solution is on singular arcs, otherwise the control field diverges.

This imposes the following constraints on \(g(t) = U(t) g U^\dagger(t)\): \[\begin{cases} \dot{g}(t) = -i[H_0 + u^*(t)H_\text{c}, g(t)], \\ \operatorname{ReTr}(g(t) H_\text{c}) = 0, \\ \operatorname{ReTr}(g(t) H_0) = 1. \end{cases}\]

On a singular arc we have \(\Phi(t)\equiv 0\). Differentiating once gives \[\begin{align} \dot{\Phi}(t) &= \operatorname{ReTr}\!\big(\dot{g}(t)\,H_{\mathrm c}\big) = \operatorname{ReTr}\!\big(-i[H_0+uH_{\mathrm c},\,g(t)]\,H_{\mathrm c}\big) \nonumber\\ &= \operatorname{ReTr}\!\big(g(t)\,i[H_0+uH_{\mathrm c},\,H_{\mathrm c}]\big) = \operatorname{ReTr}\!\big(g(t)\,i[H_0,\,H_{\mathrm c}]\big), \end{align}\] so the singularity persists only if \(\dot{\Phi}(t)=0\), i.e. \[\operatorname{ReTr}\!\big(g(t)\,i[H_0,\,H_{\mathrm c}]\big)=0.\] Differentiating again (note that \(H_0\) and \(H_{\mathrm c}\) are time-independent), \[\begin{align} \ddot\Phi(t) &= \operatorname{ReTr}\!\big(\dot{g}(t)\,i[H_0,\,H_{\mathrm c}]\big) = \operatorname{ReTr}\!\big(-i[H_0+uH_{\mathrm c},\,g(t)]\,i[H_0,\,H_{\mathrm c}]\big) \nonumber\\ &= \operatorname{ReTr}\!\big(g(t)\,i[H_0+uH_{\mathrm c},\,i[H_0,\,H_{\mathrm c}]]\big) \nonumber\\ &= \operatorname{ReTr}\!\big(g(t)\,i[H_0,\,i[H_0,\,H_{\mathrm c}]]\big) \;+\;u(t)\,\operatorname{ReTr}\!\big(g(t)\,i[H_{\mathrm c},\,i[H_0,\,H_{\mathrm c}]]\big). \end{align}\] Enforcing \(\ddot\Phi(t)=0\) on the singular arc yields the (second-order) singular control: \[u^*(t) = -\,\frac{\operatorname{ReTr}\!\big(g(t)\,[H_0,\,[H_0,\,H_{\mathrm c}]]\big)}{\operatorname{ReTr}\!\big(g(t)\,[H_{\mathrm c},\,[H_0,\,H_{\mathrm c}]]\big)} \label{eq:singular95u95second95order}\tag{32}\] provided the denominator is nonzero. If \([H_{\mathrm c},\,[H_0,\,H_{\mathrm c}]] = 0\), one must differentiate \(\Phi(t)\) further until the control appears (higher-order singular arc).

2 MAGICARP algorithm↩︎

MAGICARP iteratively refines an adjoint matrix \(g\) to optimize the target quantum operation in the interaction picture. The algorithm proceeds as follows:

1. Initialize adjoint matrix \(g\).

2. Compute initial controls using: \[u_k(0) = \frac{1}{2}\operatorname{ReTr}(g H_k).\]

3. Compute the unitary evolution using discretized pulses: \[U(\delta t) = \exp\left(-i \delta t \sum_k u_k(0) H_k\right).\]

4. Update control fields at subsequent steps using adjoint propagation: \[\tilde{u}_k(\delta t) = \frac{1}{2} \operatorname{ReTr}\left(U(\delta t) g U^\dagger(\delta t) H_k \right).\]

5. Iterate until final time \(T\).

6. Evaluate fidelity and update \(g\) using gradient ascent until convergence.

[fig::flowMAGICARP] shows the flowchart of the MAGICARP algorithm.

Let us compare the numerical requirements and advantages of MAGICARP and GRAPE, the two algorithms that we have discussed in this chapter.

2.0.0.1 Number of optimization parameters

GRAPE requires \(N_{\text{steps}}\times N_{\text{controls}}\) real optimization parameters, where \(N_{\text{steps}}\) is the number of time steps for the discretization of the evolution and \(N_{\text{controls}}\) is the number of control fields. MAGICARP on the other hand requires \(d^2-1\) real optimization parameters, where \(d\) is the dimension of the system. This is because MAGICARP optimizes the adjoint matrix \(g\), which is a \(d\times d\) traceless hermitian matrix, and the control fields are then computed from \(g\) through the self-iterative method. The dimension of the optimization space can then be very different for the two methods, and is summarised in 1.

MAGICARP can be very useful for systems with a large number of control fields, and/or a large number of time steps, it also yields, by construction, more continuous-looking pulses. However, MAGICARP scales poorly with the dimension of the system, as the number of optimization parameters grows quadratically with the dimension of the system, on the other hand GRAPE scales linearly with the number of control fields and time steps, which, for systems with a small number of control fields relatively to the dimension of the system, can yield a smaller number of optimization parameters than MAGICARP, in particular as the number of time steps can always be adjusted and therefore allows for a trade-off between the number of optimization parameters and the accuracy of the optimization for example.

2.0.0.2 Numerical complexity and stability

The numerical complexity of the two algorithms is also different. GRAPE relies on analytically computed gradients of the cost functional with respect to the control fields, which can be evaluated efficiently. MAGICARP on the other hand relies on the numerical computation of the gradients of the cost functional with respect to the adjoint matrix \(g\), which can be computationally expensive.

An interesting point can be made about the numerical stability of the MAGICARP method, as the adjoint matrix \(g\) is a traceless hermitian matrix, it is possible to use the Lie algebra \(\mathfrak{su}(d)\) to represent \(g\) as a linear combination of the generators (up to a global prefactor \(i\)) of \(SU(d)\), which could maybe be used to compute the gradients of the cost functional with respect to \(g\) efficiently.

Moreover, the dependence of the final evolution operator \(U(N\delta)\) with \(N\) the number of time steps, on the initial guess for the adjoint matrix \(g\) can be a source of numerical instability. \(U(N\delta) \in SU(d)\), therefore it has \(d^2-1\) real degrees of freedom, and the adjoint matrix \(g\) has also \(d^2-1\) real degrees of freedom, which means that the optimization space is the same for the two matrices. However, the mapping between the two spaces is not straightforward and highly non-linear, and an in-depth study of the stability of \(U(N\delta t)\) or at least \(\require{physics} \Tr(U^\dagger_\text{targ}U(N\delta t))\) with respect to small variations in \(g\) is necessary. For example, can two different adjoint matrices \(g_1\) and \(g_2\) that are very close in the optimization space lead to two very different final evolution operators \(U(N\delta)\)?

Finally a crucial point is the stability and convergence of the solution with respect to \(\delta t\), the time step for the discretization of the evolution. Does the solution converge to the optimal continuous solution for \(\delta t\rightarrow 0\)? A possible improvement to the algorithm could be to use a variable time step, where the time step is decreased as the optimization progresses, in order to refine the solution for example. Otherwise, the discretization of the evolution can be accounted for in the PMP in the first place for example, as in the work of Dionis et al. [8].

2.0.0.3 Optimality of the solution

A straightforward advantage of MAGICARP is that it combines the insights from the PMP to restrict the mathematical structure of the optimal control fields with respect to a Lagrange type cost functional. The gradient-ascent is then used to adjust this mathematical form in order to minimize the Mayer-type cost functional that is the fidelity of the obtained gate. Therefore, MAGICARP aims for an time or energy-optimal solution, while basic GRAPE only aims for a solution that minimizes the fidelity, moreover if one includes other costs in GRAPE, the gradient of the new cost functional with respect to the control fields amplitude at each time step is not straightforward to compute, and the optimization process can be very slow.

2.0.1 Example of results↩︎

As an example, consider the two-qubit system and the Hadamard gate \(H = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\) as the target gate, the control Hamiltonians are chosen to be \(\sigma_x\) and \(\sigma_y\). The results of the optimization process for the Hadamard gate using MAGICARP are shown in [fig::bloch95sphere95control95pulses95MAGICARP] where it is compared with the constrained control field from [fig::bloch_sphere_control_pulses]. The MAGICARP constrained control field is more continuous and has a duration of \(\sim1.25\tau_{QSL}\) compared to the constrained control field from [fig::bloch_sphere_control_pulses] that has a duration of \(\sim1.33\tau_{QSL}\). The latter was actually calculated using an augmented GRAPE algorithm with 2 time-steps, with the term \[\left(1-\int_0^T \left(\sum_k u_k^2\right)dt\tau_\text{QSL}^{-1}\right)^2\] added to the cost functional \(\require{physics} 1-\Tr(U_{\text{targ}}^\dagger U(T))\) to ensure that the duration of the gate is minimized.

As an example, let us consider the optimization of the Hadamard gate for systems of dimension \(d = 2, 3, 4, 5, 6\). For each dimension, the control Hamiltonians in the interaction picture are chosen to be the \(2(d-1)\) generalized Pauli matrices \(\sigma^x_{k,k+1}\) and \(\sigma^y_{k,k+1}\) ¹. The adjoint matrix \(g\) is initialized randomly, and the optimization process is repeated 300 times for each dimension.

The results of this optimization process for the QFT gate using MAGICARP are shown in [fig::cost95duration95qu2-6it-Chrestenson], where the cost functional \(\require{physics} 1-\frac{1}{d}\Tr(U_{\text{targ}}^\dagger U(T))\) is plotted as a function of the gate duration \(T\) in units of \(\Omega_{\text{max}}^{-1}\), where \(\Omega_{\text{max}}\) represents the maximum Rabi field driving the system. The results are also presented in [fig::cost95duration95qu2-6it-Chrestenson95LG], where the durations are given in units of \(\tau_{QSL}=\pi/\Omega_\text{max}\left(1-\frac{1}{d}\right)\), the quantum speed limit time for the QFT(\(d\)) gate.

The results indicate that as the dimension increases, the fidelity of achieving the target gate decreases, and the required duration increases. This is expected because the optimization space grows quadratically with the system’s dimension. Consequently, the numerical optimization is more likely to converge to local minima, rather than the global minimum, as the dimension increases. This also explains why, at higher dimensions, the optimization results are more scattered and the process becomes less stable.

It is noteworthy however that, for every dimension, there appears to be a minimal achievable duration for the gate. The dashed lines in [fig::cost95duration95qu2-6it-Chrestenson] [fig::cost95duration95qu2-6it-Chrestenson95LG] indicate, for each dimension, the minimal duration at which a cost functional of \(\require{physics} 1-\frac{1}{d}\Tr(U_{\text{targ}}^\dagger U(T))\) of \(10^{-7}\) or lower is reached. A significant fraction of the runs achieve durations close to this minimal duration, particularly for lower dimensions. This suggests that the optimization process is often able to find solutions that are near-optimal and possibly close to the quantum speed limit time with the chosen constrained set of control fields. The latter is different than the previously discussed \(\tau_\text{QSL}\), the quantum speed limit with unconstrained set of controls.

If verified, it would be interesting to analyze how this minimal duration evolves with the system’s dimension, especially in comparison to the quantum speed limit time with unconstrained control fields. As it stands, the minimal duration appears to deviate further from the quantum speed limit time with unconstrained control fields as the system’s dimension increases. This can be expected since, while the set of control Hamiltonians satisfies the Lie Algerba Rank Condition (LARC), the ratio of these control Hamiltonians to the total number of \(SU(d)\) generators decreases as the system’s dimension increases. Specifically, this ratio is \(2(d-1)/(d^2-1) = 2/(d+1)\).

2.0.2 Conclusive remarks↩︎

The MAGICARP algorithm seems to provide some insights on the quantum speed limit with a constrained set of control pulses. In particular, in a system where the control hamiltonians are the \(2(d-1)\) generalized Pauli-X and Pauli-Y matrices between two adjacent levels of the system, the minimal achievable durations seem to increase with the dimension, and to stray from optimality relative to an unconstrained set of linearly independent \(d^2-1\) control Hamiltonians as the dimension increases. This would imply that in order to minimize the duration of a gate with increased dimension, a higher connectivity is required, i.e. more allowed and controllable transitions.

The results also suggest that the optimization process is often able to find solutions that are near-optimal and possibly close to the quantum speed limit time with the chosen constrained set of control fields however. The MAGICARP algorithm proves quite promising, and, provided the necessary improvements and studies f the numerical stability and convergence are made, it could be a valuable tool for the optimization of quantum gates in the future.

References↩︎