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., Reference Awschalom, Bassett, Dzurak, Hu and Petta2013). 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., Reference Glaser, Boscain, Calarco, Koch, Köckenberger, Kosloff, Kuprov, Luy, Schirmer and Schulte-Herbrüggen2015). 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., Reference Fu, Lee, Chen, Lim, Wu, Lin, Wei, Tsao, Chang and Fann2007; Koch et al., Reference Koch, Boscain, Calarco, Dirr, Filipp, Glaser, Kosloff, Montangero, Schulte-Herbrüggen, Sugny and Wilhelm2022; de Arquer et al., Reference de Arquer, Talapin, Klimov, Arakawa, Bayer and Sargent2021).
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., Reference Schirmer, Jonckheere and Langbein2018a, Reference Schirmer, Jonckheere, O.’Neil and Langbein2018b; Langbein et al., Reference Langbein, Schirmer and Jonckheere2015). The controllers D(|IN⟩,|OUT⟩) for this scheme are designed to maximize the fidelity |⟨OUT|U(T)|IN⟩|2 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., Reference Langbein, Schirmer and Jonckheere2015, Reference Langbein, O’Neil and Shermer2022b). 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. (Reference Schirmer, Jonckheere and Langbein2018a), 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., Reference Jonckheere, Schirmer and Langbein2018).
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, Reference Schirmer and Wang2010) and exploited to design backaction-based stabilization schemes (Ticozzi et al., Reference Ticozzi, Schirmer and Wang2010; Motzoi et al., Reference Motzoi, Halperin, Wang, Whaley and Schirmer2016) has been considered in other work (Schirmer et al., Reference Schirmer, Langbein, Weidner and Jonckheere2022).
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., Reference Langbein, Schirmer and Jonckheere2015; Schirmer et al., Reference Schirmer, Jonckheere and Langbein2018a; O’Neil et al., Reference O’Neil, Langbein, Jonckheere and Shermer2023) the Hamiltonian for this spintronic network in the single-excitation subspace is represented as
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., Reference Joel, Kollmar and Santos2013).
Assuming weak interaction with the environment, the dynamics of the system are described by the Lindblad differential equation:
where ℏ is the reduced Planck constant, H D the Hamiltonian defined above and L D is a Lindblad superoperator
By setting V D = 0 we recover the usual Hamiltonian dynamics considered in previous work (Schirmer et al., Reference Schirmer, Jonckheere and Langbein2018a; O’Neil et al., Reference O’Neil, Langbein, Jonckheere and Shermer2023). 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,
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., Reference D’Alessandro, Jonckheere and Romano2014, Reference D’Alessandro, Jonckheere and Romano2015). 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., Reference D’Alessandro, Jonckheere and Romano2014) as it appears at first glance. For a Hermitian V D , it is easily verified that the Lindblad superoperator simplifies to
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
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
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
Since $\sum _k \Pi _k = I_{\mathbb{C}^N}$ , the full solution is
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., Reference Langbein, Schirmer and Jonckheere2015; Schirmer et al., Reference Schirmer, Jonckheere and Langbein2018a), 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
The maximum transfer fidelity, F max = 1, is attained when ρ(T) = ρOUT as
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) = ρ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, Reference Schirmer and Solomon2004). We thus consider three different scenarios for optimal controller synthesis:
-
1. Optimize the transfer fidelity under unitary dynamics.
-
2. Optimize a convex combination of unitary transfer and asymptotic transfer fidelity
-
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., Reference Langbein, Schirmer and Jonckheere2015; Schirmer et al., Reference Schirmer, Jonckheere and Langbein2018a; Langbein et al., Reference Langbein, O’Neil and Shermer2022b). 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., Reference Langbein, Schirmer and Jonckheere2015, Reference Langbein, O’Neil and Shermer2022a).
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,
where ρ(t) = U(t)ρ 0 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 ρ 0 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 ρ 0 and ρ OUT only involve the first two eigenvectors of H D , we can restrict ourselves to considering the representation on this subspace,
The general state evolves as
In terms of the first objective (fidelity under unitary dynamics), γ 12 = γ 21 = 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 γ 12 = γ 21 > 0,
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 ρ 0 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., Reference Burhenne, Jacob and Henze2011). 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., Reference Gorini, Frigerio, Verri, Kossakowski and Sudarshan1978; Oi and Schirmer, Reference Oi and Schirmer2012). We further normalize each dephasing matrix Γ kℓ (S) as
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 2 basis matrices for Hermitian operators on the space ${{\mathbb{C}}^{N}}$ (Altafini and Ticozzi, Reference Altafini and Ticozzi2012). The result of this vectorization process is to transform the master equation (2) to the LTI form
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
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 2 − 1 generalized Pauli matrices (Bertlmann and Kramer, Reference Bertlmann and Philipp2008) 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, Reference Floether, de Fouquieres and Schirmer2012):
Here [⋅,⋅] is the matrix commutator and {⋅,⋅} is the anti-commutator. The solution to (15) is
where r 0 ∈ ${{\mathbb{R}}^{{N^2}}}$ is the vectorized version of ρ 0 with components given by r 0k = Tr(ρ 0 σ 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 2 = 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., Reference Schirmer, Jonckheere, O.’Neil and Langbein2018b; Jonckheere et al., Reference Jonckheere, Schirmer and Langbein2018; O’Neil et al., Reference O’Neil, Schirmer, Langbein, Weidner and Jonckheere2022), 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
or in the LTI formulation as
The fidelity error for the perturbed density matrix of (19) is then given as
or
where c is the N 2 × 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., Reference O’Neil, Schirmer, Langbein, Weidner and Jonckheere2022). We rely on the log-sensitivity based on its widespread use and relation to fundamental limitations of classical control (Jonckheere et al., Reference Jonckheere, Schirmer and Langbein2018; Dorf and Bishop, Reference Dorf and Bishop2000). 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
Note that we depart from the definition of log-sensitivity used in (O’Neil et al., Reference O’Neil, Schirmer, Langbein, Weidner and Jonckheere2022), 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., Reference O’Neil, Schirmer, Langbein, Weidner and Jonckheere2022) and includes the uncertain couplings J as in (Jonckheere et al., Reference Jonckheere, Schirmer and Langbein2018). 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 = 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., Reference Schirmer, Jonckheere, O.’Neil and Langbein2018b), 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, Reference Silverman1986). The bandwidth used in the calculation is h = 3.5σn −1/3 as derived in (Scott, Reference Scott1979) 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., Reference Perperoglou, Sauerbrei, Abrahamowicz and Schmid2019). A numeric differentiation of ê(T; S,δ), evaluated at the point where δ = 0, then provides an estimate of the differential sensitivity. We then have
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., Reference Langbein, O’Neil and Shermer2022c) contains the entire collection of these figures for all controllers and optimization options.
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, Reference Dorf and Bishop2000)
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
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 0: no correlation between s a (S, T) and s k (S, T);
-
H 1: 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 0: no correlation between e(T) and {s a (S, T), s k (S, T)};
-
H 1: 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, Reference Abdi and Salkind2007) where n = 100 denotes the number of samples. We then determine the statistical significance of the test by evaluating
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
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 0 if p τ, r ≥ α,
-
reject H 0, 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 0 in support of H 1 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.
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 0 in favor of H 1 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−4 in the figure) we still observe a nearly linear trend. However, in the right-hand side of the plot (e(T) > 10−4 in the figure), we observe a large number of controllers with varying log-sensitivity for the same error.
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.
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.
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.
Data availability statement
The data is available at (Langbein and others, Reference Langbein, O’Neil and Shermer2022a) and (Langbein and others, Reference Langbein, O’Neil and Shermer2022c).
Financial support
Sean O’Neil acknowledges PhD funding through the US Army Advanced Civil Schooling program.
Competing interests
The authors report no conflicts of interest.
Comments
No accompanying comment.