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

(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., Reference Joel, Kollmar and Santos2013).

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., 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,

(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., 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

(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., 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

(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|.

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., 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

(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, Reference Altafini and Ticozzi2012). 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, 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):

(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., 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

(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., 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

(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., 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. 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

(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., Reference Langbein, O’Neil and Shermer2022c) 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, Reference Dorf and Bishop2000)

(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).