Robustness of energy landscape control to dephasing

Abstract

As shown in previous work, in some cases closed quantum systems exhibit a non-conventional absence of trade-off between performance and robustness in the sense that controllers with the highest fidelity can also provide the best robustness to parameter uncertainty. As the dephasing induced by the interaction of the system with the environment guides the evolution to a more classically mixed state, it is worth investigating what effect the introduction of dephasing has on the relationship between performance and robustness. In this paper we analyze the robustness of the fidelity error, as measured by the logarithmic sensitivity function, to dephasing processes. We show that introduction of dephasing as a perturbation to the nominal unitary dynamics requires a modification of the log-sensitivity formulation used to measure robustness about an uncertain parameter with nonzero nominal value used in previous work. We consider controllers optimized for a number of target objectives ranging from fidelity under coherent evolution to fidelity under dephasing dynamics to determine the extent to which optimizing for a specific regime has desirable effects in terms of robustness. Our analysis is based on two independent computations of the log-sensitivity: a statistical Monte Carlo approach and an analytic calculation. We show that despite the different log-sensitivity calculations employed in this study, both demonstrate that the log-sensitivity of the fidelity error to dephasing results in a conventional trade-off between performance and robustness.

Introduction

The advent of quantum technology and promise of applications ranging from quantum computing to quantum sensing has resulted in strong interest in a range of quantum systems. In particular coupled spin systems, or spin networks for short, show potential as simple prototypes on the path to scaling to more complex systems (Awschalom et al., 2013). As control plays a fundamental role in the translation of physical phenomena into technology, the development and implementation of effective control schemes for quantum systems are essential to harness their technological potential (Glaser et al., 2015). Coupled with the design of controllers for quantum systems, tools to assess and guarantee robustness of these controllers to the effect of the environment are essential to realizing the benefits of quantum technology for high fidelity medical imaging, operation of quantum gates and quantum computing (Fu et al., 2007; Koch et al., 2022; de Arquer et al., 2021).

A paradigm for quantum control based on energy landscape shaping has been proposed and applied to derive feedback control laws for selective transfer of excitations between different nodes in a spin network (Schirmer et al., 2018a, 2018b; Langbein et al., 2015). The controllers D(|IN⟩,|OUT⟩) for this scheme are designed to maximize the fidelity |⟨OUT|U(T)|IN⟩|² of transfer from an input state |IN⟩ to an output state |OUT⟩ at a specified readout time T. U(T) is the unitary time evolution operator of the system, which depends on static fields to shift the energy levels of the system (Langbein et al., 2015, 2022b). These optimal controllers D are selective in that no input other than |IN⟩ can drive the system to |OUT⟩ at time T.

Although this quantum control problem can be formulated as a linear time-invariant (LTI) control system with state feedback, there are numerous differences between this quantum control problem and a classical tracking problem. The unitary evolution of a closed quantum system is characterized by persistent oscillations. As such, the system is not stable in a classical sense, and the target states are not attractive or asymptotically stable steady states. Lack of stability might be expected to bode poorly for robustness but this is not necessarily always the case. In Schirmer et al. (2018a), it was shown that controllers achieving perfect state transfer also have vanishing sensitivity with respect to perturbations to the coupling strengths of the drift Hamiltonian. At the same time, statistical analysis of a set of optimal energy landscape controllers for uniformly coupled spin rings, ranging in size from 3 to 20 spins, found that under certain conditions, a concordant relationship between the error and the log-sensitivity is possible for controllers that achieve low but nonvanishing errors (Jonckheere et al., 2018).

Our goal in this paper is to investigate whether a nonconventional lack of trade-off between robustness and performance is observed in simple spin rings evolving under dephasing dynamics or whether a conventional trade-off exists. We assess robustness through the log-sensitivity, calculated both analytically and via kernel density estimation (KDE). We also consider controllers optimized to provide maximum fidelity under different assumptions on the system–environment interaction. In particular, we consider controllers optimized for maximum fidelity under unitary evolution, those optimized for fidelity with dephasing introduced and those optimized for a linear combination of fidelity under unitary dynamics and steady-state fidelity. To understand this effect of decoherence on the controller design we will focus on the intermediate regime where coherent dynamics play a significant role but are modified by dephasing as a result of weak interactions with the environment. The strongly dissipative regime where asymptotic stability can be recovered (Schirmer and Wang, 2010) and exploited to design backaction-based stabilization schemes (Ticozzi et al., 2010; Motzoi et al., 2016) has been considered in other work (Schirmer et al., 2022).

While the non-conventional absence of a trade-off between robustness and performance observed in some cases under coherent dynamics may carry over to systems subject to decoherence, the addition of decoherence alters the dynamics significantly. Pure dephasing, in particular, results in quantum superposition states converging to classical mixed states. It is thus reasonable to expect a more classical, in a control-theoretic sense, behavior for systems subject to decoherence.

In Section “Coupled spin dynamics – single-excitation subspace,” we introduce the theory of coupled spin systems and their evolution under decoherence in the single-excitation subspace. In Section “Control objectives and controller design,” we introduce the control objectives in terms of maximization of the transfer fidelity and the objective functions for optimal controller design in the different regimes. Next, in Section “LTI formulation of the state equation,” we introduce an LTI form of the dynamical equations amenable to robustness analysis. In Section “Robustness assessment,” we provide the pair of methods (statistical and analytic) used to gauge the sensitivity and robustness of the controllers. In Section “Analysis,” we present results detailing the level of concordance between the two log-sensitivity calculations, the degree to which the robustness properties of the controllers agree with the trade-offs from classical control theory, and explore the effect of specific controller types on the fidelity and observed robustness properties. Finally, we conclude with Section “Conclusion.”

Methods

In this section, we outline the equations governing excitation transfer for spin rings in the single-excitation subspace, describe the control objectives and optimization scheme used to develop the controllers, and detail the methods used to assess the robustness of these controllers to perturbations in the form of dephasing.

Coupled spin dynamics – single-excitation subspace

We consider a ring of N spin-1/2 particles with nearest-neighbor coupling in the subspace of the state space where the total number of excitations is one, which consists of states where one spin is in an excited state and N − 1 spins remain in the ground state, and superpositions of such states. As detailed in (Langbein et al., 2015; Schirmer et al., 2018a; O’Neil et al., 2023) the Hamiltonian for this spintronic network in the single-excitation subspace is represented as

(1) $${H_D} = \left( {\matrix{ {{D_1}} & {{J_{1,2}}} & 0 & \ldots & 0 & {{J_{1,N}}} \cr {{J_{1,2}}} & {{D_2}} & {{J_{2,3}}} & {} & 0 & 0 \cr 0 & {{J_{2,3}}} & {{D_3}} & {} & 0 & 0 \cr \vdots & {} & \ddots & \ddots & \ddots & \vdots \cr 0 & 0 & 0 & {} & {{D_{N - 1}}} & {{J_{N - 1,N}}} \cr {{J_{1,N}}} & 0 & 0 & \ldots & {{J_{N - 1,N}}} & {{D_N}} \cr } } \right)$$

in the basis where each natural basis vector represents the excitation localized at spin |n⟩. The terms J _{k, ℓ} represent the coupling between spins |k⟩ and |ℓ⟩ and are all assumed equal. The terms D _n are scalar values of the time-invariant control fields applied to shape the energy landscape. This single-excitation subspace model is a simplification of the model for a system with any number of excited spins up to N. Specifically, this is the subspace that results by retaining only those eigenvectors with eigenvalue 1 for the total spin operator $S_N={1 \over 2}\sum _{k = 1}^{N}(I+Z_k)$ (Joel et al., 2013).

Assuming weak interaction with the environment, the dynamics of the system are described by the Lindblad differential equation:

(2) $${\rm\dot{\rho}} (t) = - {{j \over \hbar }}[{H_D},{\rm{\rho}} (t)] + {L_D}({\rm{\rho}} (t)),$$

where ℏ is the reduced Planck constant, H _D the Hamiltonian defined above and L _D is a Lindblad superoperator

(3) $${L_D}({\rm{\rho}} ) = {V_D}^\dagger {\rm{\rho}} {V_D} - {{1 \over 2}}({V_D}^\dagger {V_D}{\rm{\rho}} + {\rm{\rho}} {V_D}^\dagger {V_D}).$$

By setting V _D = 0 we recover the usual Hamiltonian dynamics considered in previous work (Schirmer et al., 2018a; O’Neil et al., 2023). Here, we study systems subject to decoherence that can be modeled as dephasing in the Hamiltonian basis. This is a common model for weak decoherence, described by a Lindblad operator L _D of dephasing type, given by a Hermitian dephasing operator that commutes with the system Hamiltonian,

(4) $$ V_D ={V_D}^\dagger {\rm \ where\ }[H_D,V_D] = 0. $$

The subscript D in the Lindbladian indicates a dependence on the control as, strictly speaking, decoherence in the weak coupling limit depends on the total Hamiltonian, and hence on the control (D’Alessandro et al., 2014, 2015). Although this is a simple decoherence model, it is closer to the master equation in the weak coupling limit developed in (D’Alessandro et al., 2014) as it appears at first glance. For a Hermitian V _D , it is easily verified that the Lindblad superoperator simplifies to

(5) $${L_D}({\rm{\rho}} ) = - {{1 \over 2}}[{V_D},[{V_D},{\rm{\rho}} ]].$$

As H _D and V _D commute for the dephasing model, they are simultaneously diagonalizable and there exists a set of projectors {Π _k } _k onto the (orthogonal) simultaneous eigenspaces of H _D and V _D such that $\sum _{k}\Pi _k=I_{{\mathbb{C}}^{N}}$ is a resolution of the identity on the Hilbert space ${{\mathbb{C}}^{N}}$ of the single-excitation subspace and

$$ H_D = \sum _{k} \lambda _k \Pi _k, \quad V_D = \sum _{k} c_k \Pi _k, $$

where λ _k and c _k are the real eigenvalues of H _D and V _D , respectively.

Pre- and post-multiplying the master equation (2) with Lindblad term (3) by Π _k and Π _ℓ yields

(6) $$ \Pi _k\dot{{\rm{\rho}} }(t) \Pi _\ell = (-j \omega _{k \ell } - \gamma _{k \ell }) \Pi _k {\rm{\rho}} (t) \Pi _\ell, $$

where ${\omega _{k\ell }} = {{1 \over \hbar }}({\lambda _k} - {\lambda _\ell })$ and ${\gamma _{k\ell }} = {{1 \over {2\hbar }}}{({c_k} - {c_\ell })^2} \ge 0$ . The solution for the projection Π _k ρ(t)Π _ℓ is

$$ \Pi _k{\rm{\rho}} (t)\Pi _\ell = e^{-t(j \omega _{k\ell } + \gamma _{k\ell })} \Pi _k {\rm{\rho}} _0 \Pi _\ell. $$

Since $\sum _k \Pi _k = I_{\mathbb{C}^N}$ , the full solution is

(7) $$ {\rm{\rho}} (t) = \sum _{k,\ell } e^{-t(j \omega _{k\ell } + \gamma _{k\ell })} \Pi _k {\rm{\rho}} _0 \Pi _\ell. $$

Control objectives and controller design

In this section, we define the control objective as the transfer fidelity and discuss differing conditions under which we seek to maximize this measure. These varying conditions manifest as distinct sets of controllers aimed at optimizing the fidelity under differing conditions.

Transfer fidelity

Following the framework adopted in earlier work (Langbein et al., 2015; Schirmer et al., 2018a), we seek static controllers that map an input state to a desired output state by shaping the energy landscape of the system. Specifically, our design objective is to find a controller that steers the dynamics to maximize the transfer fidelity of an excitation at an initial node of the network, |IN⟩, to an output node |OUT⟩ at a specific readout time T. If the output state is a pure state, this target state can be represented as ρ_OUT = |OUT⟩⟨OUT|. We then evaluate the fidelity of the state ρ(t) at time T in terms of the overlap with ρ_OUT as

(8) $$ F[{\rm{\rho}} (T)] = {\rm Tr} [{\rm{\rho}} _{\rm{OUT}} {\rm{\rho}} (T)]. $$

The maximum transfer fidelity, F _max = 1, is attained when ρ(T) = ρ_OUT as

(9) $$ F[{\rm{\rho}} (t)] \leq F[{\rm{\rho}} _{\rm{OUT}}] = {\rm Tr}[{\rm{\rho}} _{\rm{OUT}}^2] = 1. $$

The fidelity error e(T) at the readout time T is therefore given by e(T) = 1 − F[ρ(T)]. We thus seek controllers that maximize this transfer fidelity (equivalently minimize the fidelity error) defined in (8), where ρ(T) is the solution of Eq. (2) with ρ(0) = ρ₀ = |IN⟩⟨IN|.

Optimal controller design

For the energy landscape control paradigm, finding a controller is equivalent to finding an ordered N-tuple of control parameters ${\left\{ {{D_1},\;{D_2}, \ldots, \;{D_N}} \right\}^T}$ and a time T that maximizes the transfer fidelity for a system evolving under the total Hamiltonian H _D , according to Eq. (2). If the decoherence process is known precisely then it is straightforward to optimize the control objective by evolving the system according to Eq. (2) and evaluating the fidelity. However, the exact dephasing rates for a given system are often not precisely known (Schirmer and Solomon, 2004). We thus consider three different scenarios for optimal controller synthesis:

1. 1. Optimize the transfer fidelity under unitary dynamics.

2. 2. Optimize a convex combination of unitary transfer and asymptotic transfer fidelity

3. 3. Optimize the fidelity averaged over a sampling of decoherence processes.

Option 1 is a reasonable choice if decoherence is a weak perturbation to the Hamiltonian dynamics. Optimizing solely for asymptotic transfer fidelity may be a reasonable choice if the decoherence is so strong that the system is likely to reach a steady state before the transfer is complete. However, in the intermediate regime, when we are unsure of the dephasing rates, optimizing for Options 2 and 3 is more practical.

Optimization under unitary dynamics

Optimization of Option 1 has been considered in previous work (Langbein et al., 2015; Schirmer et al., 2018a; Langbein et al., 2022b). The optimization problem in the other cases can be solved similarly, using standard optimization algorithms with suitable modification of the objective functional. Despite the complex optimization landscape, we have found that the L-BFGS (limited memory Broyden–Fletcher–Goldfarb–Shanno) quasi-Newton algorithm with restarts using randomly selected initial values in a sufficiently large domain based on stratified sampling works well for all options (Langbein et al., 2015, 2022a).

Optimization of coherent and asymptotic transfer fidelity

To simultaneously optimize coherent and asymptotic transfer, we define an objective function that is a weighted average of both, for example,

(10) $$ \alpha {\rm Tr}[{\rm{\rho}} (T) {\rm{\rho}} _{\rm{OUT}}] + \left (1-\alpha \right ) {\rm Tr}[{\rm{\rho}} _{\infty }{\rm{\rho}} _{\rm{OUT}}], $$

where ρ(t) = U(t)ρ ₀ U(t)^† is the initial state propagated by unitary evolution and ρ _∞ is the steady state of decoherent evolution. Both ρ(T) and ρ _∞ are efficiently calculated from (7) by setting t = T, γ _kℓ = 0 for coherent evolution and t = ∞, k = ℓ for the decoherent steady state, respectively.

To maximize this weighted average fidelity, a controller must be superoptimal. Specifically, it must enable perfect state transfer from ρ ₀ to ρ _OUT at time t = T. Simultaneously, it must maximize the overlap of the steady state ρ _∞ with the target state ρ _OUT. To see that such controllers exist in principle, consider a controller that achieves maximum asymptotic transfer by rendering the input and output states orthogonal superposition states of the form ${\left | \textrm{IN} \right \rangle } = {1 \over \sqrt{2}} \left ({\left | e_1 \right \rangle } +{\left | e_2 \right \rangle } \right )$ and ${\left | \rm{OUT} \right \rangle } = {1 \over \sqrt{2}} \left ({\left | e_1 \right \rangle } -{\left | e_2 \right \rangle } \right )$ in the eigenbasis of the Hamiltonian, $H_D = \sum _{k=1}^N E_k{\left | e_k \right \rangle }{\left \langle e_k \right |}$ . We refer to these states as orthogonal pairs. Here, E _k and |e _k ⟩ represent the energy eigenvectors and eigenvalues of H _D . Since ρ ₀ and ρ _OUT only involve the first two eigenvectors of H _D , we can restrict ourselves to considering the representation on this subspace,

(11) $${{\rm{\rho}} _0} = {1 \over {\sqrt 2 }}\left[ {\matrix{ 1 & 1 \cr 1 & 1 \cr } } \right],\;\;\;{\kern 1pt} {{\rm{\rho}} _{{\rm{OUT}}}} = {1 \over {\sqrt 2 }}\left[ {\matrix{ 1 & { - 1} \cr { - 1} & 1 \cr } } \right].$$

The general state evolves as

(12) $${\rm{\rho}} (t) = {1 \over {\sqrt 2 }}\left[ {\matrix{ 1 & {{e^{ - t(j{\omega _{12}} + {\gamma _{12}})}}} \cr {{e^{ - t(j{\omega _{21}} + {\gamma _{21}})}}} & 1 \cr } } \right].$$

In terms of the first objective (fidelity under unitary dynamics), γ ₁₂ = γ ₂₁ = 0, and at a time $T = { (2n + 1) \pi \over \omega _{12}}$ for n ∈ ${\mathbb{Z}}$ the controller achieves perfect state transfer or ρ(t) = ρ _OUT, maximizing the first half of the objective. In terms of the asymptotic component of the objective, with γ ₁₂ = γ ₂₁ > 0,

(13) $$\mathop {\lim }\limits_{t \to \infty } {\rm{\rho}} (t) = {1 \over {\sqrt 2 }}\left[ {\matrix{ 1 & 0 \cr 0 & 1 \cr } } \right].$$

And so ${\rm Tr}\left [ {\rm{\rho}} _{\infty }{\rm{\rho}} _{\textrm{OUT}} \right ] = {1 \over 2}$ , maximizing the possible overlap.

Decoherence-averaged optimization

Optimization of the transfer fidelity, averaged over many dephasing processes, is in principle also straightforward. For a given initial state ρ ₀ and controller D, the output state ρ ^(D,S) (T) subject to a dephasing process S is calculated according to (7) and transfer fidelity from (8). From this the average transfer fidelity can be computed by taking the mean of the transfer fidelity over all decoherence processes.

The computational overhead of the average fidelity evaluation, and thus the optimization as a whole, depends linearly on the number of decoherence processes averaged over. Efficient sampling of the possible decoherence processes to minimize the number of required decoherence processes and avoid sampling bias is therefore important.

Decoherence in the form of dephasing in the Hamiltonian basis is modeled by sampling the space of pure dephasing processes. We generate a large set of N × N lower triangular matrices with entries Γ _kℓ ^(S) in [0,1]. Here, N is the system dimension given by the number of qubits in the network. To ensure even sampling of the whole space, the entries of the triangular matrices are drawn from a Sobol sequence for low-discrepancy sampling, thus allowing an even covering of the sample space (Burhenne et al., 2011). A set of at least 10, 000 dephasing operators is then generated by eliminating all trial dephasing matrices that violate the complete-positivity physical constraints for evolution of an open quantum system (Gorini et al., 1978; Oi and Schirmer, 2012). We further normalize each dephasing matrix Γ _kℓ ^(S) as

(14) $$ \bar{\Gamma }_{k\ell }^{(S)} = \Gamma _{k\ell }^{(S)} \left/ \sum\limits _{\substack{1 \lt k \le N 1 \le \ell \lt k}} |\Gamma _{k\ell }^{(S)}|\right .. $$

We then use 1, 000 of these dephasing operators in each optimization run for fidelity with dephasing dynamics.

LTI formulation of the state equation

To facilitate the following analysis we reformulate the Lindblad equation (2) and its solution (7) through expansion by a set of N ² basis matrices for Hermitian operators on the space ${{\mathbb{C}}^{N}}$ (Altafini and Ticozzi, 2012). The result of this vectorization process is to transform the master equation (2) to the LTI form

(15) $$ \dot{r}(t) = (A + L)r(t) $$

where r(t) ∈ ${{{\mathbb{R}}^{N^2}}}$ is the vectorized version of the density matrix in the Hamiltonian basis, A ∈ ${{\mathbb{R}}^{{N^2\times}{N^2}}}$ is the matrix representation of the Liouville superoperator which determines the unitary evolution of the system, and L ∈ ${{\mathbb{R}}^{{N^2\times}{N^2}}}$ is the matrix representation of the Lindblad superoperator. To see this, consider the decomposition of the controlled Hamiltonian H _D = UΛU ^† where U ∈ U(N) and Λ ∈ ${{\mathbb{R}}^{{N\times}{N}}}$ is the diagonal matrix of real eigenvalues of H _D . Pre-multiplying (2) by U ^†, post-multiplying by U and noting that UU ^† = U ^† U = I, we have

(16) $$\hbar \dot{\tilde{{\rm{\rho}}}} (t) = j \Lambda {\rm{\tilde{\rho}}} (t) + j {\rm{\tilde{\rho}}} (t)\Lambda + C {\rm{\tilde{\rho}}} (t)C - {{1 \over 2}}{C^2} {\rm{\tilde{\rho}}} (t) - {{1 \over 2}} {\rm{\tilde{\rho}}} (t){C^2}$$

where ρ῀(t) = U ^† ρ(t)U is the representation of the density matrix in the Hamiltonian basis, and C = U ^† V _D U is a diagonal matrix of the N scalar c _k ’s from the decomposition $V_D = \sum _{k}c_k \Pi _k$ .

We choose the N ² − 1 generalized Pauli matrices (Bertlmann and Kramer, 2008) complemented by I _N to form an orthonormal basis for the Hermitian operators on ${{\mathbb{C}}^{N}}$ which we designate as {σ _n }. The orthonormality conditions are expressed as Tr(σ _m ,σ _n ) = δ _mn . Expansion of (16) in terms of the basis {σ _n }, yields the following (Floether and others, 2012):

(17a) $${r_k}(t) = {\rm{Tr}}\left( {\rm{\tilde{\rho}}} {(t){\sigma _k}} \right),$$

(17b) $${A_{k\ell }} = {\rm{Tr}}\left( {{j \over \hbar }\Lambda [{\sigma _k},{\sigma _\ell }]} \right),$$

(17c) $${L_{k\ell }} = {1 \over \hbar }{\rm{Tr}}(C{\sigma _k}C{\sigma _\ell }) - {1 \over {2\hbar }}{\rm{Tr}}({C^2}\{ {\sigma _k},{\sigma _\ell }\} ).$$

Here [⋅,⋅] is the matrix commutator and {⋅,⋅} is the anti-commutator. The solution to (15) is

(18) $$ r(t) = e^{t(A + L)}r_0 $$

where r ₀ ∈ ${{\mathbb{R}}^{{N^2}}}$ is the vectorized version of ρ ₀ with components given by r _0k = Tr(ρ ₀ σ _k ).

Before proceeding, we make a few observations. Firstly, as C and Λ are diagonal and commute, then so do A and L which simplifies calculation of the log-sensitivity. Secondly, given the requirement that Tr(ρ(t)) = 1 ∀t, we see that $r_{N^2}(t) = {1 \over \sqrt{N}}$ , which implies ṙ _{N ²} = 0. This guarantees the existence of at least one zero eigenvalue of the dynamical equation (15) and, similarly, the existence of a subspace along which the trajectory is constant (a steady state).

Additionally, if the eigenvalues of A + L not equal to zero are distinct, and the dephasing process is characterized by N distinct jump operators not equal to zero, then A + L has N zero eigenvalues corresponding to k = ℓ, consistent with the constant populations on the main diagonal of ρ(t) in (7). We return to this LTI formulation to assess robustness in Section “Robustness assessment.”

Robustness assessment

Following (Schirmer et al., 2018b; Jonckheere et al., 2018; O’Neil et al., 2022), we assess the robustness of the implemented controllers to perturbations to the nominal system through the logarithmic sensitivity of the fidelity error, or log-sensitivity for short.

We consider our nominal system as the closed system evolving under unitary dynamics with a nominal trajectory given by (7) with γ _kℓ = 0 or the vectorized form (15) with $L = {\bf 0}_{N^2 \times N^2}$ . We then consider perturbations of this nominal trajectory due to the introduction of dephasing. We represent the dephasing process as the matrix S _μ ∈ ${{\mathbb{C}}^{{N^2\times}{N^2}}}$ where μ is used to index distinct dephasing processes. Even though S _μ is an operator on the Hilbert space over ${\mathbb{C}}^{{N^2}}$ , it consists of strictly real elements Γ _kℓ ^{(S _μ )} = γ _kℓ ^{(S _μ )} as defined in (14). We will consider a set of 1000 dephasing operators in each trial, indexed as μ ∈ {1, 2, …, 1000}. Furthermore, we introduce a parameter δ ∈ [0,1] to modulate the strength of the dephasing process for each S _μ . We then denote the perturbed trajectory in the density matrix formalism as

(19) $$ \tilde{{\rm{\rho}} }(t;S_{\mu },\delta ) = \sum _{k,\ell }e^{-t(j\omega _{kl} - \delta \gamma ^{(S_{\mu })}_{k\ell })} \Pi _{k}{\rm{\rho}} _0 \Pi _{\ell } $$

or in the LTI formulation as

(20) $$ \tilde{r}(t;S_{\mu },\delta ) = e^{t(A + \delta S_{\mu })}r_{0}. $$

The fidelity error for the perturbed density matrix of (19) is then given as

(21) $$ \tilde{e}(T;S_{\mu },\delta ) = 1 - {\rm Tr}[\tilde{{\rm{\rho}} }(T;S_{\mu },\delta ){\rm{\rho}} _{\rm{OUT}}] $$

or

(22) $$ \tilde{e}(T;S_{\mu },\delta ) = 1 - \boldsymbol ce^{(A + \delta S_{\mu })T}r_0 $$

where c is the N ² × 1 row vector corresponding to the transpose vectorized representation of ρ _OUT .

We are now in a position to assess the robustness of $\tilde{e}(T;S_\mu, \delta)$ by measuring the effect of the dephasing perturbation through the log-sensitivity (O’Neil et al., 2022). We rely on the log-sensitivity based on its widespread use and relation to fundamental limitations of classical control (Jonckheere et al., 2018; Dorf and Bishop, 2000). Specific to the excitation transfer problem, we quantify performance as high fidelity F[ρ(T)] (equivalently low error 1 − F[ρ(T)] = e(T)) and measure robustness as the size of the differential change in the performance induced by an external perturbation or ${\partial e(T) \over \partial \delta }$ . In terms of fundamental limitations, we expect that those controllers which exhibit the best performance will display the greatest differential log-sensitivity (the least robustness) and vice versa for those controllers exhibiting low performance.

Next, in order to better gauge a normalized percentage change in $\tilde{e}(T;S_\mu, \delta)$ for a given S _μ , we consider a logarithmic sensitivity of the form ${\partial \ln (\tilde{e}(T) \over \partial \delta }$ . We calculate the log-sensitivity of the fidelity error to the perturbation S _μ as

(23) $$ s(S_{\mu },T) =\left. {\partial \ln (\tilde{e}(T;S_\mu,\delta )) \over \partial \delta } \right |_{\delta = 0} = \left. {1 \over e(T)}{\partial \tilde{e}(T;S_{\mu },\delta ) \over \partial \delta } \right |_{\delta = 0}. $$

Note that we depart from the definition of log-sensitivity used in (O’Neil et al., 2022), as we see two distinct cases when applying the log-sensitivity. In the first case there are no changes in the inertia of the system matrix when parameters drift about their nominal values. By inertia of a system matrix, we mean the number of eigenvalues λ _k with ℜλ _k = 0, ℜλ _k < 0 and ℜλ _k > 0, characterizing the response of the system. This first case is the basis for the analysis in (O’Neil et al., 2022) and includes the uncertain couplings J as in (Jonckheere et al., 2018). In such cases, the genuine log-sensitivity $$\partial \ln(e(T)) / \partial \ln(\xi) \rvert_{\xi = \xi_0}$$ provides a meaningful, dimensionless measure of sensitivity to the uncertain parameter $$\xi$$ with nominal value $$\xi_0$$ . In the more complicated second case there are changes of inertia around the nominal parameter values. The present case falls in this category as the nominal decoherence rate is δ ₀ = 0. Introduction of the perturbation induces a bifurcation in the system dynamics from unitary evolution to decoherent evolution. In this case, evaluating the log-sensitivity as ∂ln (e(T))/∂ln (δ)| _{δ = δ 0} yields a zero value for all controllers. This requires a revision of the log-sensitivity as (∂e(T)/∂δ) ⋅ (1/e(T))| _{δ = 0} to obtain a meaningful log-sensitivity. Regardless, s(S _μ ,T) as defined in (23) provides a percentage change in the error with respect to the introduction of dephasing. Put differently, the value of s(S _μ ,T)δ for nonvanishing δ provides a means to compare the effect of decoherence process S _μ of strength δ across different controllers.

Kernel density estimate approach

As shown in (Schirmer et al., 2018b), evaluating the log-sensitivity of the error numerically over a large number of perturbations provides one method of assessing robustness. In the current study, we use this as the first approach to calculating the log-sensitivity. Specifically, we consider spin rings of size N = 5 and N = 6. For each ring, we consider transfers from spin |IN⟩ = |1⟩ to |OUT⟩ = |2,3⟩ for N = 5 and |OUT⟩ = |2,3,4⟩ for N = 6. For each transfer we select the best, as measured by highest nominal fidelity, 100 controllers from the three optimization categories described in Section “Optimal controller design”: fidelity under unitary dynamics, fidelity under dephasing, and unitary transfer combined with asymptotic fidelity. For brevity, we refer to these optimization options as fidelity, dephasing and overlap, respectively. For each controller, we select 1000 dephasing operators generated by the process described in Section “Decoherence-averaged optimization.” For each dephasing process, we consider a perturbation δ ∈ [0,1] quantized into 1001 points with a uniform interval of 0.001.

With this setup, we calculate $$\tilde{e}(T;S_\mu, \delta)$$ for each controller across the entire population of S _μ and δ by (19) and (21). The result is an 1000 × 1001 element array of error results arranged by dephasing process along the rows and perturbation strength along the columns. This array serves as the kernel for the MATLAB function ksdensity to compute a kernel density estimate of the error to dephasing for a given strength (Silverman, 1986). The bandwidth used in the calculation is h = 3.5σn ^−1/3 as derived in (Scott, 1979) where σ is the standard deviation of the samples in each error array noted above and n is the number of samples. For each step in δ we also calculate the mean and variance of the error across the 1000 dephasing processes. These 1001 samples of the mean error over all dephasing operators serves as input to the MATLAB function fit with option “smoothingspline” to produce a functional representation of the mean error designated as ê(T; S,δ) where we drop the subscript μ on the dephasing operator to indicate that the averaging process in the density estimation has already been taken into account. We choose a smoothing spline over a cubic fit to provide the greatest degree of freedom while providing a functional fit that minimizes any distortion in the data (Perperoglou et al., 2019). A numeric differentiation of ê(T; S,δ), evaluated at the point where δ = 0, then provides an estimate of the differential sensitivity. We then have

(24) $$ s_{k}(S,T) = {1 \over e(T)}\left. {\partial \hat{e}(T;S,\delta ) \over \partial \delta } \right |_{\delta = 0} $$

where e(T) = ê(T; S,0) is the nominal error. The subscript k indicates the value of the log-sensitivity is calculated from the density estimate of the mean error.

As an example of the output of this KDE, Figure 1 displays a heatmap visualization of the error versus decoherence strength in the upper pane. The slope of the green line at the point where δ = 0 provides the estimate of the differential sensitivity used in the log-sensitivity calculation. The repository located at (Langbein et al., 2022c) contains the entire collection of these figures for all controllers and optimization options.

Figure 1 Top: heat map of KDE-based fidelity error distribution as a function of the dephasing strength for a 0 → 1 transfer in a 5-ring with dephasing-optimized controller. The green lines indicate the mean and standard deviation of the distribution. The slope of the mean error as a function of the decoherence strength δ is used to estimate the sensitivity in the limit δ → 0, which provides the numerical estimate s _k (S,T) of the log-sensitivity. Bottom left: values of the biases of the energy landscape controller, with transfer time T and steady-state overlap indicated in the title. Bottom center and right: heat map of the initial state ρ ₀ and final state ρ(T) with respect to an eigenbasis of the controlled Hamiltonian. Notice that in this example, the initial and final state form an orthogonal pair in this basis.

Analytic calculation

The structure of the matrices A and S _μ greatly simplifies the calculation of the log-sensitivity. Based on the fidelity error defined in (20) and noting that A and S _μ commute we have that $$\tilde{e}(T;S_\mu, \delta) = {c}e^{AT}e^{\delta S_\mu T}r_0$$ . The log-sensitivity is then easily calculated as (Dorf and Bishop, 2000)

(25) $$ s(S_{\mu },T) = \left. {1 \over e(T)}{\partial \tilde{e}(T;S_{\mu },\delta ) \over \partial \delta } \right |_{\delta = 0} = \left \vert {-1 \over e(T)}\boldsymbol {\boldsymbol c}e^{AT}(S_{\mu }T)r_0 \right \vert. $$

For each controller, we calculate s(S _μ , T) for the same 1000 dephasing operators used in the approach of Section “Kernel density estimate approach” for μ ∈ {1, 2, …, 1000}. To arrive at a single value of the log-sensitivity for each controller and to maintain consistency with the KDE approach, we take the arithmetic mean over all dephasing operators and define

(26) $$ s_{a}(S,T) = {1 \over 1000} \sum _{\mu = 1}^{1000}s(S_{\mu },T) $$

where the subscript a indicates the log-sensitivity is derived from an analytic calculation of the perturbed trajectory of (20).

Analysis

We focus our analysis on three topics: the level of concordance between the two log-sensitivity measures s _k (S, T) and s _a (S, T), the degree to which the controllers in the data set exhibit robustness properties that align with the trade-off induced by the S + T = I identity of classical feedback control, and the role played by input–output state orthogonal pairs in the fidelity and robustness of the controllers. For the first two categories we execute a pair of hypothesis tests based on the Kendall τ as a nonparametric measure of correlation and the Pearson r as a measure of linear correlation. The analysis of the role of orthogonal pairs is based on a visual interpretation of the fidelity and robustness plots for controllers that render input–output states as orthogonal pairs and those that do not.

Hypothesis test

In Section “Comparison of robustness assessments” below, we test the concordance between s _a (S, T) and s _k (S, T). We expect a positive correlation between the two metrics and test the level of concordance through a one-tailed hypothesis test with the right tail. We establish

• H ₀: no correlation between s _a (S, T) and s _k (S, T);

• H ₁: positive correlation between s _a (S, T) and s _k (S, T).

In Section “Robustness trend analysis,” we compare the trend between each measure of log-sensitivity (s _a (S, T) and s _k (S, T)) with the nominal error e(T). The classical control trade-offs require a negative correlation between e(T) and the log-sensitivity, so we establish a one-tailed hypothesis test on the left tail. For this test:

• H ₀: no correlation between e(T) and {s _a (S, T), s _k (S, T)};

• H ₁: negative correlation between e(T) and {s _a (S, T), s _k (S, T)}.

In both cases we consider the N = 5 transfers from |IN=1⟩ to |OUT⟩ = |2,3⟩ and N = 6 transfers from |IN⟩ = |1⟩ to |OUT⟩ = |2,3,4⟩. For each of these five transfers we consider the controllers optimized for fidelity, dephasing and overlap. With 100 controllers within each transfer-optimization target combination we thus have 15 tests for each hypothesis above, each with 100 samples.

To execute the computation of the Kendall τ and Pearson r we leverage the MATLAB function corr(⋅,⋅) with the option “Kendall” or “Pearson” as appropriate. In the following discussion of hypothesis tests, use of τ refers to the Kendall rank correlation coefficient and r to the Pearson correlation coefficient. We compute the test statistic for the Kendall τ as $Z_{\tau } = \tau \left ( \sqrt{{2(2n+5) \over 9n(n-1)}} \right )^{-1}$ (Abdi, 2007) where n = 100 denotes the number of samples. We then determine the statistical significance of the test by evaluating

(27) $${p_\tau } = \left\{ {\matrix{ {1 - \Phi ({Z_\tau }),} \hfill & {{\rm{ for }}\ {s_k}(S,T)\,{\rm{ vs}}{\rm{. }}\,{s_a}(S,T),} \hfill \cr {\Phi ({Z_\tau }),} \hfill & {{\rm{ for }}\ \{ {s_a}(S,T),\,{s_k}(S,T)\} \,{\rm{ vs}}{\rm{. }}\,e(T),} \hfill \cr } } \right.$$

where Φ(⋅) is the normal cumulative distribution function. For the Pearson r-based test, we calculate the test statistic as $t_{r} = r \left ( \sqrt{{1-r^{2} \over n-2}}\right )^{-1}$ . Here again, n = 100. We then quantify the statistical significance of the test for a given value of r as

(28) $${p_r} = \left\{ {\matrix{ {1 - {\cal S}({t_r}),} \hfill & {{s_k}(S,T)\,{\rm{ vs}}{\rm{. }}\,{s_a}(S,T),} \hfill \cr {{\cal S}({t_r}),} \hfill & {\{ {s_a}(S,T),{s_k}(S,T)\} \,{\rm{ vs}}{\rm{. }}\,e(T),} \hfill \cr } } \right.$$

where $\mathcal{S}$ represents the cumulative Student’s t-distribution.

Finally, we establish the level of significance at α = 0.02 so that the hypothesis test itself is

• accept H ₀ if p _{τ, r} ≥ α,

• reject H ₀, if p _{τ, r} < α,

for both tests.

Comparison of robustness assessments

The hypothesis test reveals strong agreement between s _a (S, T) and s _k (S, T) for all 15 tests with a Kendall τ or Pearson r of near unity in all cases. Table 1 summarizes the results for the Kendall τ-based hypothesis test. As depicted, all 15 test cases result in rejection of H ₀ in support of H ₁ with p-values near zero. The Pearson r-based test provides similarly strong evidence of linear correlation. Though not indicative of equality, the hypothesis test confirms strong concordance. Visually, we can see this strong correlation in Figure 2, where both measures lie nearly on-top of the other. While the analytic calculation is the preferred method of assessing robustness, evidence of this correlation supports the efficacy of the KDE measure when the equations of motion are so complex that available computing power makes this Monte Carlo approach more efficient than an analytic calculation.

Figure 2 Comparison s _a (S, T) (solid blue line) vs. s _k (S, T) (dashed red line) for fidelity-optimized controllers in an N = 6 ring with transfer to spin 3. Note the near perfect agreement of the two log-sensitivity measures. We can also observe here the conventional trend between both log-sensitivity measures and the fidelity error indicating a trade-off between performance (low e(T)) and robustness (small s _{{a, k}}(S, T)).

Table 1. Results of Kendall τ-based hypothesis test for the concordance of s _a (S, T) and s _k (S, T). The hypothesis test provides strong confirmation that both the analytic and kernel density estimation are consistent in the evaluation of robustness to dephasing

Robustness trend analysis

The results of the hypothesis test to evaluate concordance of the log-sensitivity with the fidelity error are summarized in Table 2. The table presents the results of the Pearson r-based test, which evaluates the level of linear correlation between the metrics on a log − log scale. For all 15 test cases we see rejection of H ₀ in favor of H ₁ with p-values near zero. This suggests a strong negative correlation between the two metrics. Furthermore, the Pearson r provides the slope of the best linear fit through the data point, valuable for assessing the impact of a given change in error on robustness. Figure 3(a) shows a plot of s _a (S, T) and s _k (S, T) versus e(T) on a log-log scale for a 5-ring, nearest-neighbor transfer and fidelity-optimized controller. The near unity slope of the plot and r value of − 0.972 indicates a nearly uniform cost in robustness, measured by the log-sensitivity, for a given increase in performance, quantified by the fidelity error. Conversely, Figure 3(b) depicts the correlation for controllers optimized for fidelity in a 6-ring for nearest-neighbor transfer. The linear correlation coefficient is less strong than that in Figure 3(a) with r = − 0.5882, indicating a less stringent adherence to a uniform cost in robustness for increased performance. The visual plot confirms this. For the highest-fidelity controllers (e(T) < 10⁻⁴ in the figure) we still observe a nearly linear trend. However, in the right-hand side of the plot (e(T) > 10⁻⁴ in the figure), we observe a large number of controllers with varying log-sensitivity for the same error.

Figure 3 s _a (S, T) and s _k (S, T) versus e(T) on a log − log scale. Plot (a) suggests a strong linear correlation, which is confirmed by a Pearson r of − 0.972, while (b) shows a much weaker linear correlation with r = − 0.589.

Table 2. Table summarizing the Pearson r-based hypothesis test results for log-sensitivity versus e(T). The test results in rejection of H ₀ in favor of H ₁ in all cases, indicating a trend in agreement with the conventional trade-off between performance and robustness

Orthogonal pair states, robustness and fidelity

Given the suitability of orthogonal pairs to maximize the combined objective of fidelity under unitary transfer with asymptotic transfer fidelity, we expect that controllers which render the input and output states orthogonal pairs will dominate the “overlap” controller data set. Of more interest, however, are what fidelity and robustness properties input–output orthogonal pairs present across the breadth of the controllers. As Table 3 confirms, input–output states that form orthogonal pairs produce the best (highest-fidelity) controllers for the case of overlap optimization. For both 5 and 6 rings, all overlap-optimized controllers for nearest-neighbor and next-nearest-neighbor (i.e., |OUT⟩ = |3⟩) transfers create an orthogonal-pair input and output state. However, for the N = 6 and 1 → 4 transfer, the population of overlap-optimized controllers is almost evenly split. For both fidelity-optimized and dephasing-optimized controllers, the controllers that render input–output orthogonal pairs dominate the nearest-neighbor transfers for both N = 5 and N = 6. However, for all other transfers, orthogonal pairs constitute a minority of the controllers. Given that the controllers under consideration were filtered for the best performance (fidelity), this suggests that controllers facilitating the highest levels of fidelity for nearest-neighbor transfer, regardless of optimization choice, exhibit this orthogonal-pair property. This suggests optimizing to produce input–output eigenstructures that replicate orthogonal pairs as a means to generate high fidelity controllers under a range of conditions for nearest-neighbor transfer.

Table 3. Table depicting the percentage of input–output states rendered as orthogonal pairs by each type of controller and each transfer

The robustness properties of the orthogonal-state controllers are less clear. For fidelity-optimized and dephasing-optimized controllers and nearest-neighbor transfer, the lower sensitivity controllers exhibit the orthogonal pair property as seen in Figure 4(a). Whereas for non-nearest-neighbor transfers, the more robust controllers appear to not render the input–output states orthogonal pairs. This behavior is evident in Figure 4(b) where the most robust fidelity-optimized controllers in a 6-ring for the |1⟩ → |4⟩ transfer are not of the orthogonal pair variety. Whether the robustness properties of the controllers in these cases are a characteristic of the orthogonal-like controllers or simply due to the dominance of orthogonal-type controllers in nearest-neighbor transfers still requires further investigation.

Figure 4 Plot of s _a (S, T) vs. e(T) on a log − log scale showing controllers that yield orthogonal-pair input–output states in the eigenbasis of the Hamiltonian (blue crosses) and those that do not (red boxes). Orthogonal pairs dominate, show the greatest fidelity and also exhibit the smallest log-sensitivity in (b) while in (b) both the highest-fidelity controllers and most robust controllers are not of the orthogonal-pair variety.

Conclusion

In this paper we applied two distinct approaches to evaluate the log-sensitivity of the fidelity error to a perturbation of the system dynamics in the form of dephasing in the Hamiltonian basis. The KDE approach based on the error measurements of 1000 dephasing operators introduced at varying strength and with the bandwidth input as described in Section “Kernel density estimate approach” produced log-sensitivity values nearly identical to those of the analytic calculation based on the same dephasing operators. As distinct from previous work, we also considered controllers optimized for not only fidelity under unitary dynamics but fidelity under dephasing and fidelity in the asymptotic regime. We showed that in all cases the relationship between the fidelity error and the log-sensitivity adheres to the trade-off between performance and robustness expected from classical control theory. Of note, though we examined controllers optimized to maximize fidelity under dephasing, these controllers exhibited no better robustness to dephasing than the other optimization choices. This suggests that optimizing for fidelity under decoherence does not accrue robustness benefits, as measured by the log-sensitivity, over optimizing for fidelity under unitary transfer without exact knowledge of the dephasing process.

Several items still require further investigation. Firstly, we only considered perturbations in the form of dephasing. It is necessary to generalize the analytic formula to allow for perturbations in the form of dephasing simultaneously with uncertainty in the Hamiltonian or control parameters as well as to extend the analysis to consider general decoherence processes that include dissipation with dephasing. This would facilitate a better understanding of robustness for general open quantum systems, a persistent challenge within quantum control. Next, an investigation of the robustness properties of the controllers that produce input–output state orthogonal pairs is in order. In particular, though these types of controllers tend to provide the best fidelity for nearest-neighbor transfer, a comparison of their robustness properties with controllers that do not provide this orthogonal-pair property is missing, mainly due to the paucity of non-orthogonal-type controllers for nearest-neighbor transfers for the controllers considered in this study. An investigation into this matter would provide a justification for either pursuing controllers of another sort, should they provide greater robustness for the same fidelity, or relying on orthogonal-type controllers to provide the best fidelity and robustness for nearest-neighbor transfer and optimizing with that target as the goal.