Noise Resistant Quantum Control Using Dynamical Invariants

A systematic approach to design robust control protocols against the influence of different types of noise is introduced. We present control schemes which protect the decay of the populations avoiding dissipation in the adiabatic and non-adiabatic regimes and minimize the effect of dephasing. The effectiveness of the protocols is demonstrated in two different systems. Firstly we present the case of population inversion of a two level system in the presence of either one or two simultaneous noise sources. Secondly, we present an example of the expansion of coherent and thermal states in harmonic traps, subject to noise arising from monitoring and modulation of the control respectively.


Introduction
A major obstacle to manipulate quantum systems and develop quantum technologies is the unavoidable presence of noise. There are different types of noise sources that disturb control protocols, including noise induced by the environment, noise caused from interaction with the measurement apparatus, or errors in the implementation of the control protocol. Different approaches proposed to suppress or mitigate the effects of noise include (for a recent review see [1]): the use of decoherence-free subspaces which are immune to noise [2], correction of errors using quantum feedback controls [3], performing sudden interactions on the system on timescales for which the noise only slightly interferes with the process such as in dynamical decoupling [4], and applying shortcuts to adiabaticity (STA) [5,6]. The noise has also been proposed as a resource to achieve the desired control in specific processes [7,8,9].
A large family of control protocols are based on adiabatic following of instantaneous eigenstates of a time dependent Hamiltonian by smoothly and slowly changing control parameters. These methods are widespread as they are in principle robust against control imperfections. However, they are also prone to suffer the effects of noise due to long operation times. As a result the fidelity of the final state with respect to the target state is reduced [1,10].
STA provide a strategy to combat the effects of noise thanks to two different mechanisms: (i) In principle faster than adiabatic processes are desirable to avoid pervasive, long interactions with the environment. In practice the fidelities may present maxima at specific times [8], and STA can be set for these optimally short process times. Many studies have tested the achieved robustness with master equations including the noise, see e.g. [35,36]; (ii) in addition, the parameter paths leading to STA are typically not unique, so this freedom may be used to choose the most robust ones with respect to specific perturbations, noise or control imperfections. This optimization has been performed so far by minimizing the excitation energy or maximizing the fidelity in perturbative schemes, e.g. in two-level systems [37,38,39,40] or for ion transport [41], using decoherence free subspaces [42,43], super-operator [44] and non-Hermitian invariants [45], and effective Hamiltonians [46] .
In this work we introduce an alternative systematic method for smooth control under the influence of noise which is applicable for both adiabatic and nonadiabatic time scales. The technique proposed intends to go beyond the perturbative regime and can be applied to the strong noise regime [47]. The central idea is to inverse engineer the noiseless Hamiltonian by designing its dynamical invariants (i.e., the dynamics of the noiseless system) such that the noise has a minimal effect. The control functions to be minimized are state independent, and measure the deviation of the actual invariant, i.e., the one for the full dynamics including noise, from the noiseless invariant. This technique does not require one to solve the full dynamics iteratively as is often done in optimal control methods [48,49]. This property makes the method appealing and simple for implementation. Furthermore, it is not restricted to very fast operations, where tipycally very short control time is limited by experimental constraints.
The dynamics of the system including the dissipative term resulting from the noise takes the form (for = 1) where HereĤ(t) is the total Hamiltonian of the system including the control Hamiltonian and the noise term given by Eq. (2). TheX k (t) represent Hermitian operators acting on the Hilbert space of the system and can be explicitly time dependent. The pre-factors η k are scaling factors representing the strength of the noise and may have different dimension depending onX k (t). The sum over k includes the possibility of independent types of noise simultaneously affecting the dynamics. This equation was derived in different contexts, including the singular coupling limit [50], phase noise [51], action noise [8], amplitude noise [52], noise from monitoring weakly some quadrature of the system [3,53], Gaussian noise and Poisson noise for SU(2) algebra [54,55], and more [56,57].
In Sec. 2 we present the main results of this work. First we construct the dynamical invariant method in the density operator formalism which can also be then applied to the study of noise and naturally extend the treatment from pure states to general mixed states. We derive two measures to quantify the effect of noise introducing constraints on the noiseless dynamical invariant. In Sec. 3 we study the example of the two-level system with single and multiple noise terms in the dynamics. Section 4 studies the control of thermal and coherent states of the harmonic oscillator. We conclude with a discussion and outlook on future work in Sec. 5, plus some technical appendices.

Dynamical Invariant for Unitary Dynamics
We refer the reader who is unfamiliar with the dynamical invariant method to Appendix A where we present the method in the wave function formalism.
For noiseless, unitary dynamics, the evolution of the density operator is described by A dynamical invariant satisfies the equation [58] i and can be expressed in diagonal form, Here λ k are the real time independent eigenvalues and |φ k ≡ |φ k (t) are the time dependent eigenvectors of the invariant. In this basis, the density matrix elements ρ lk ≡ φ l |ρ(t) |φ k can be calculated froṁ where the dot represents the time derivative. The off diagonal terms of the density matrix depend on the difference of time derivatives of two Lewis-Riesenfeld phases (compare with Eq.(A.5)), while the populations remain constant with time [58]. As the system is driven through the instantaneous eigenstates of the invariant, imposing [Î(0),Ĥ(0)] = [Î(t f ),Ĥ(t f )] = 0 we ensure that the system starts and ends in an energy eigenstate ofĤ without unwanted excitations. The state transfer is designed by choosinĝ I(t) and then determiningĤ(t). (See Appendix A for more details.)

Dynamical invariant under the influence of noise
We now consider the influence of Eq. (2) on the control process. In order to demonstrate the effect of noise we consider a single operatorX ≡X 1 (t) and η ≡ η 1 . In a later example, we will also consider the case for simultaneous noise sources. The dynamics for an arbitrary observableÂ including the noise effect in the Heisenberg representation reads Assuming the structure of the invariant (5) for the unitary dynamics we insert it in Eq. (7) to account for the noise. The eigenvalues of the invariant λ l are now no longer constant in time and evolve according tȯ Note that if {|φ k } are eigenstates ofX thenλ l = 0 as required in the unitary noiseless method. In this case the invariant is not affected by the noise. Although the requirement forÎ andX having common eigenstates cannot generally be achieved for all times during the process, the effect of noise can be significantly reduced by constructing the Hamiltonian from an invariant which shares common eigenvectors with those of X during most of the process. Since at final time we impose that the invariant and the Hamiltonian share common eigenstates, protecting the invariant from the noise will drive the system to the desired target state.
To express the density matrix elements, Eq. (6) is now modified by adding to the r.h.s. the additional termρ d lk ≡ φ l |Lρ|φ k which accounts for dissipation and decoherence resulting from the noise term and is given bẏ − 2 nm ρ nm φ l |X|φ n φ m |X|φ k .
In the limit thatÎ andX share common eigenstates the contribution due to noise to the change in population and the off diagonal terms iṡ where {x k } are the eigenvalues ofX. Note that in this limit the decay of the populations in the invariant eigenbasis is suppressed, however, a decay of the coherences is still present, although it can be minimized for sufficiently fast processes. The strategy proposed here relies on the ample freedom provided by STA. By adding constraints on the unitary invariant we can design a control Hamiltonian that optimizes the fidelity under the influence of the noise.
To identify the amount of overlap between the two bases sets ofÎ andX we define the overlap matrix S with the entries Here {|ψ j } are the eigenvectors ofX. The sum of the overlap matrix can be bounded by: n n ij |S ij | < n 2 , where n is the dimension of the space. The upper bound is not tight, and obtaining a tight bound, typically, becomes a difficult optimization problem for high dimension. However, we are only interested in minimizing S. In the scenario of n = 2, the tight upper bound is given by 2 √ 2. Next, we define the measure for the overlap along the process as the time average of the distance between the overlap matrix and its minimal value.
The measure O is zero if and only if the eigenbasis ofX andÎ are identical. A different measure which stems from similar considerations and in some cases can be easier to compute is In the above expressions we use the Frobenius norm defined as M ≡ tr (MM † ). The normalization factor z = 2 guarantees that A is dimensionless and equal or smaller than 1 (this is an immediate consequence of the sub-additivity of the Frobenius norm). The measure A is zero if and only ifX andÎ commute at all times during the processÎX =XÎ, and takes the maximal value A = 1 whenÎX = −XÎ.
For unbounded operators extra care is needed. The norm should be calculated on a finite domain or using other techniques as will be demonstrated in a later section. In order to improve the fidelity of the evolved state with respect to the target by minimizing the effect of noise, the controls that drive the system are inverse engineered through the invariantÎ(t) of the unitary dynamics subject to the minimization of the measures O or A.

Two-level system
As a first example we consider the control problem of a full population inversion in a two level system (TLS) [16,59,60,61,62]. The Hamiltonian takes the form: where ∆(t) and Ω(t) are real, time-dependent functions resulting from an interaction with some external field, andσ z andσ x are the Pauli matrices. Initially the system is set to the ground state,ρ 0 = |0 0|, with the initial Hamiltonian corresponding to ∆(0) = ∆ 0 and Ω(0) = 0. The desired target stateρ tar = |1 1| corresponds to the ground state of the final Hamiltonian ∆(t f ) = −∆ 0 and Ω(t f ) = 0. To evaluate the success of the control protocol we will use the fidelity that measures the overlap between the final state and the target stateρ tar . To connect the statesρ(0) andρ tar we engineer the controls ∆(t) and Ω(t) from the dynamical invariant. Associated with the Hamiltonian (14) there is a dynamical invariant expressed as (see Appendix B), where G ≡ G(t) and B ≡ B(t) are auxiliary real time dependent functions of the invariant obeying The frictionless conditions [ . At intermediate times these two functions are totally free. In particular interpolating G(t) = 3 i=0 g i t i and B(t) = 3 i=0 b i t i by polynomials with at least the same degree as the number of boundary conditions lets us deduce from Eq. (17) the desired controls ∆(t) and Ω(t). However, extra-coefficients can be added to the interpolation, for example, Here g 4 can be used to also control the values of the measures O or A.

Single noise source
We first consider amplitude noise of the formX 1 ≡σ z and η 1 ≡ η z . For this particular type of noise the measure O is given explicitly by which is independent of the free function B of the invariant. It takes its minimal value O z → 0 when G(t) → nπ with n ∈ Z and the maximal 2 √ 2 − 2 when G(t) → π/2 + nπ. Similarly we can write explicitly (13) after some simple algebra: As for the measure O z the minimal value A z → 0 occurs when G → nπ and A z → 1 for G → π/2 + nπ.
In Figs. 1 we plot the fidelity as function of the measures O z and A z for a given final time t f for the full population inversion problem. Both measures show a similar behavior, when O z → 0 and A z → 0 the fidelity is improved significantly and monotonically decreasing as O z and A z increase. Thus, by adding constraints on these measures when constructing the invariant we obtain a control field which minimizes the effect of the noise. When the dynamics is subject to noise from a single source, i.e., a singleX 1 , a control protocol which leads to fidelity 1 can be found. Generally, when the dynamics is subject to several independent sources of noise,X j , obtaining fidelity 1 is not guaranteed. Nevertheless, the influence of the overall noise can still be minimized and the final fidelity is improved.

Multiple noise sources
When multiple noise terms (differentX j ) are present in the dynamics, the measures O j or A j cannot always be minimized simultaneously. In the next example we study the worst case scenario where the two noise terms have mutually unbiased bases [63]. In this case minimization of one of the noise terms will lead to maximization of the other. In particular we consider amplitude noise both in theσ z andσ x fields, i.e.,X 1 =σ z , η 1 = η z andX 2 =σ x , η 2 = η x . For the TLS we employ O z(x) to quantify the effect of the noise in the dynamics. As for O z we can write explicitly O x as By examining the integrands of Eqs. (18) and (20) we observe that (i) when G(t) → nπ, For multiple noise terms we suggest to minimize the averageŌ of the single noise measures weighted according to their relative strength. In the example above this average reads, with a minimum value (see Appendix C for more details) This implies that in order to optimize the fidelity, protocol (i) or (ii) are chosen depending on the relation of the noise strength η z and η x . Thus, minimizing the influence of the stronger noise term will lead to higher fidelity as is demonstrated in Fig. 2. In this figure we plot the fidelity against O z and O x for different η z and η x ratios. Maximal fidelity is obtained when the averageŌ is minimal and given by Eq. (22). We remark that equivalently, optimization can be performed using the measure A by the replacement of O → A in Eq.(21).

Quantum Harmonic Oscillator
In this section we study the quantum harmonic oscillator which for example can describe a particle with reduced mass m (for the simulations the mass of 100 ions of 40 Ca + is used) in a harmonic trap with time dependent frequency ω(t). The Hamiltonian takes the form,Ĥ Note that this problem can be mapped to a general control problem of the SU(1,1) algebra [5,14,58,64,65]. Thus, Eq.(23) can be written as, where we define a = 1/m, b(t) = mω 2 (t), and identifyT 1 =p 2 /2,T 2 =q 2 /2 and T 3 = (pq +qp)/2 that satisfy the commutation relations Associated with the harmonic oscillator Hamiltonian (23) there is a dynamical invariant of the form (see Appendix D) where [x,π] = i withx ≡q/ρ,π ≡ ρp − mρq, and ρ is an auxiliary scaling function satisfying Ermakov's equation with imposed by the frictionless conditions [Ĥ(t b ),Î(t b )] = 0 and continuity. As in the example of the TLS we use the freedom to interpolate the free function ρ at intermediate times. We choose functions of polynomials with sufficient parameters to satisfy the previous six boundary conditions. As we showed in the previous section, extra coefficients can be incorporated with higher order polynomials to impose other constraints such as the minimization of O or A.
In the next two examples, we study the expansion control of coherent and thermal states. In these cases the success of the control protocol is evaluated according to the previous fidelity definition, Eq. (15), for Gaussian states [66].

Coherent states
We assume that the initial coherent state |α with the initial frequency ω 0 = ω(0) is driven to the final target state |α with ω f = ω(t f ), whereα = αe −igω 0 and g = t f 0 dt /ρ 2 . For this end we interpolate ρ(t) = ( 5 i=0 r i t i ) −1/2 and deduce ω(t) from Eq. (27) (see Appendix D). This noise arises from weakly and continuously measuring (monitoring) the position of the particle in the trap leading toX =q [3,53]. As was discussed above, for unbounded operators the calculation of the overlap between the bases to compute O should be carried on a finite domain or as we will see next it can be evaluated using This overlap can be written explicitly as (see Appendix D) where H n are the Hermite polynomials. In principle, to compute O we should consider the sum over n from 0 to ∞ of the elements S n in Eq. (28). Nevertheless, we find that minimizing Eq. (28) for a certain n will necessarily minimizes all the different n terms. We prove this by showing that the spatial integration over q is independent of the function ρ. Proof: We preform the following coordinate substitution, u 1 = q/ρ and u 2 = 1/ √ ρ.
The determinant of the Jacobian is given by Then, Eq. (28) takes the form The integration over u 1 depends on n, but it is independent of ρ. Thus, different designs of ρ influence only the integration over u 2 which is independent of n, implying that it is sufficient to minimize Eq. (28) for an arbitrary n when constructing the invariant. In Fig. 3a we design different protocols and plot the fidelity against S 0 normalized by the maximal S 0 value out of the protocols considered in the figure. This is done by adding two extra-coefficients in the invariant interpolation ρ(t) = ( 7 i=0 r i t i ) −1/2 , where r 6 and r 7 control the values of S 0 in Eq. (28) and g that let us fix the final target coherent state independently of the ρ interpolation (see Eq. (D.11).) As S 0 becomes smaller the fidelity is enhanced. Figure 3b presents two control protocols corresponding to the green and red points of Fig. 3a. The green dashed line represents the standard STA protocol [5] (standard refers to those protocols where the free functions in the invariant are only constrained by the boundary frictionless conditions). This protocol can be improved using the method we presented, minimizing S 0 to achieve higher fidelities. We remark that higher fidelity than those shown in Fig. 3a can be achieved just if a higher order polynomial is incorporated when interpolating ρ.

Thermal states
Consider again the harmonic oscillator Hamiltonian (23), we now choose different states to protect against noise. The initial state is assumed to be the thermal state,ρ 0 = exp(−β 0Ĥ (0))/Z, with the normalization factor Z, and the initial inverse temperature β 0 and frequency ω 0 ≡ ω(0). The final Hamiltonian corresponds to the frequency ω f ≡ ω(t f ) and the target state is the thermal stateρ tar = exp(−β fĤ (t f ))/Z with the final inverse temperature β f = β 0 ω 0 /ω f . The noise considered in this example is noise in the modulation of the frequency described by the noise operatorX = 2T 2 =q 2 and constant η. Since this noise is more problematic for long operation times, a natural way to avoid it is to have short operation times. However, very short expansion times are typically not feasible experimentally. Designing protocols protected against amplitude noise improve the final fidelities even at longer times. In Fig. 4a we plot the fidelity against final times t f for three different control protocols. In blue we plot the fast adiabatic protocol of constant µ ≡ω/ω 2 [8], in red the standard STA protocol, and in green the improved STA protocol (for both STA protocols ω(t) is deduced from Eq. (27) using the following ansatzes: ρ(t) = ( 5 i=0 r i t i ) 1/2 for the standard and ρ(t) = ( 6 i=0 r i t i ) 1/2 for the improved protocols, respectively). We see that for the optimized STA protocol higher fidelity for all final times t f is obtained. This introduces high flexibility for controlling the final time and the average instantaneous power consumption/production, In Fig. 4b we plot the fidelity vs. the absolute value of the averaged instantaneous powerP.

Discussion
In this work we introduce a method to construct a control protocol which is robust against dissipation of the population and minimizes the effect of dephasing. Doing so, we optimize the fidelity of the final state with respect to the target state. As is shown in Eq.(10) the diagonal terms of the density matrix will remain constant at the end of the process while the off diagonal will be affected by dephasing at a rate proportional to the square of the distance between the eigenvalues of the noise operator. This is a clear indication that Markovian noise cannot be completely suppressed without adding an auxiliary system which will store the information about the coherence. The idea of the method presented is based on the fact that the dynamical invariant provides a family of infinite solutions from which the control Hamiltonian can be constructed for a particular state transfer problem. By imposing additional constraints on the invariant, namely, minimization of the measure O or equivalently A, we protect the invariant from the noise during the process. Since at the final time the invariant and the Hamiltonian share common eigenvalues the target state is achieved with high fidelity. The main advantages of this method compared to other control optimization methods is its simple implementation. It does not require calculations by iteration and does not involve perturbation methods. Since the structure of the invariant for the unitary dynamics is already known in many cases, imposing additional constraints on this invariant is not a difficult task. Moreover, the method is applicable to different time scales, from the nonadiabatic to the adiabatic regime. This implies that in order to suppress the noise we are not limited to frequent sudden operations which in many cases are not feasible experimentally and will typically be costly in terms of power. Formulation of the method in terms of density operator is necessary to treat noise but it also makes controlling mixed states possible as in the example of the thermal state.
The idea of protecting the invariant from the noise during the process and by that designing an optimal noise resistant control can in principle be applied to all types of noise including thermal and non-Markovian noise. If the noise operators are not Hermitian, we suggest (like in the procedure above) to find the invariant for the unitary dynamics with additional degrees of freedom which can be set later to minimize the effect of noise. Next, instead of considering the overlap between the bases which now might not be computable, we can use the rate of change of the eigenvalues subject to the full dynamics as a measure for minimization. In the limit l |λ l | → 0 the noise will not affect the invariant and a noise resistant control can be found.

Appendix A. Invariant inverse engineering based on Lie algebras
We summarize [67], a systematic approach to inverse engineering the controls from the dynamical invariants of a system when it is described by a closed Lie algebra. Let us assume that the time-dependent HamiltonianĤ(t) describing a quantum system is given by a linear combination of Hermitian generatorsT a , where the h a (t) are real time-dependent functions and theT a span a Lie algebra [68] [ with α abc the structure constants. Associated with the Hamiltonian there are timedependent Hermitian invariants of motionÎ(t) that satisfy [69] A wave function |Ψ(t) which evolves withĤ(t) can be expressed as a linear superposition of the instantaneous invariant modes [69] |Ψ(t) = n c n e iαn |φ n (t) , where the c n are constants, the phases α n fulfill and the eigenvectors |φ n (t) ofÎ(t) where λ n are the constant eigenvalues.
If the invariant is also a member of the dynamical algebra, it can be written aŝ where where the kets are defined in terms of the component of each generator [67]. Note that the relation between the Hamiltonian and the invariant is a property of the algebra, i.e. the structure constants, and is independent of the representation. Usually these coupled equations are interpreted as a linear system of ordinary differential equations for f a (t) when the h a (t) components of the Hamiltonian are known [68,70,71,72]. Here we consider a different perspective taking them as an algebraic system to be solved for the h a (t), when the f a (t) are given. As there are many Hamiltonians for a given invariant [73] we cannot generally invert Eq. (A.8) as |h = G −1 |ḟ to get |h . This means that det(G)= 0, so at least one of the eigenvalues a (i) (t) of the G matrix vanishes. Different approaches, such as Gauss elimination or projector techniques [67], can be used to find the pseudo-inverse matrix of G and deduce the Hamiltonian component |h in terms of the invariant |f .
When inverse engineering shortcuts to adiabaticity [5,6], the Hamiltonian is usually given at initial and final times. In general the invariantÎ (equivalently |f (t) ) is chosen to drive, through its eigenvectors, the initial states of the Hamiltonian H(0) to the states of the finalĤ(t f ) [5,14,69] according to Eq. (A.4). This is ensured by imposing at the boundary times t b = 0, t f , the frictionless conditions [Ĥ(t b ),Î(t b )] = 0 [5]. Equivalently, using Eqs. (A.1), (A.7), and since theT a generators are independent this condition implies As we pointed before for this algebra det(G) = 0, so Eq. (B.2) is not directly invertible.
After some simple algebra we find the h a (t) components in terms of f a (t) if the constrainṫ with all indices i, j, k different, E ijk is the Levy-Civita symbol (1 for even permutations of (123) and -1 for odd permutations), c is a constant, and h k (t) is considered a Hamiltonian free component chosen for convenience. The frictionless conditions (A.10) for this algebra is To be more specific note that the TLS in Sec. 3 is governed by this algebra with the following representation of generators, where h 1 (t) = Ω(t), h 3 (t) = ∆(t), and the boundary Hamiltonians Ω(0) = 0, ∆(0) = ∆ 0 at t = 0 and Ω(t f ) = 0, ∆(t f ) = −∆ 0 at t = t f to produce the population inversion among the |0 and |1 states. The objective is to design Ω(t) and ∆(t) to connect these two states by imposing partially the structure Eq. (14) ofĤ(t), i.e. h 2 (t) = 0 ∀t, as it is not always experimentally feasible to implement T 2 . Imposing h 2 (t) we chose to interpolate f 1 and f 2 satisfying the boundary conditions We use polynomial interpolations with at least the same degree as the number of boundary conditions, nevertheless, higher order polynomials can be considered to impose even more constraints. Then f 3 is given by Once the f a are fixed Ω(t) and ∆(t) are deduced from Eq. (B.4), An alternative and sometimes convenient choice to express the invariant is using the angles on the Bloch sphere G(t) and B(t), parametrazing f 1 (t) = Ω R sin(G(t)) cos(B(t)), f 2 (t) = Ω R sin(G(t)) sin(B(t)), and choosing Ω R = Ω 2 (0) + ∆ 2 (0). The constraint (B. 3) imposes f 3 (t) = Ω R cos(G(t)) with c = Ω 2 R , so the invariantÎ is expressed as in Eq. (16). According to Eq. (B.4) the Hamiltonian coefficients h a (t) are given in terms of these polar angles as [73] ∆ = −Ḃ +Ġ tan(G) tan(B) , Ω =Ġ sin(B) , (B.8) with the boundary conditions for population inversion G(0) = π, G(t f ) = 0,Ġ(t b ) = 0, remaining B(t b ) andḂ(t b ) as free parameters.

Appendix C. Overlap Matrix for TLS
For two arbitrary bases in the Hilbert space C 2 the overlap matrix S is bounded by In this scenario there are three bases that maximize the overlap ij |S ij |. These bases are the mutually unbiased bases given by the eigenvectors of the Pauli matrices. In the main text we study the simultaneous overlap of two mutually unbiased bases with the basis of the invariant. In particular the bases of interest are the eigenvectors ofσ z andσ x which reads {(1, 0) t , (0, 1) t } and (  The weighted average of the overlaps of the two bases with C.2 is given by 2p sin θ 2 + cos θ 2 + 2(p − 1) 1 − cos(ϕ) sin(θ) 2 + 1 + cos(ϕ) sin(θ) 2 .