T Schulte-Herbrüggen

A Spörl

N Khaneja

S J Glaser

B, special issue Quantum Control Theory for Coherence and Information Dynamics Optimal Control for Generating Quantum Gates in Open Dissipative Systems

Keywords: numbers: 03.67.-a, 03.67.Lx, 03.65.Yz, 03.67.Pp; 82.56.Jn I

Optimal control methods for implementing quantum modules with least amount of relaxative loss are devised to give best approximations to unitary gates under relaxation. The potential gain by optimal control fully exploiting known relaxation parameters against time-optimal control (the alternative for unknown relaxation parameters) is explored and exemplied in numerical and in algebraic terms: for instance, relaxation-based optimal control is the method of choice to govern quantum systems within subspaces of weak relaxation whenever the drift Hamiltonian would otherwise drive the system through fast decaying modes. In a standard model system generalising ideal decoherence-free subspaces to more realistic scenarios, opengrape-derived controls realise a cnot with delities beyond 95% instead of at most 15% for a standard Trotter expansion. As additional benet their control elds are orders of magnitude lower in power than bang-bang decouplings.

INTRODUCTION

Using experimentally controllable quantum systems to perform computational tasks or to simulate other quantum systems [1,[START_REF] Feynman | Feynman Lectures on Computation[END_REF] is promising: by exploiting quantum coherences, the complexity of a problem may reduce when changing the setting from classical to quantum. Protecting quantum systems against relaxation is therefore tantamount to using coherent superpositions as a resource. To this end, decoherence-free subspaces have been applied [START_REF] Zanardi | [END_REF], bang-bang controls [4] have been used for decoupling the system from dissipative interaction with the environment, while a quantum Zeno approach [5] may be taken to projectively keep the system within the desired subspace [6]. Controlling relaxation is both important and demanding [7,[START_REF] Lidar | Decoherence-Free Subspaces and Subsystems[END_REF]9,10,[START_REF] Viola | Quantum Communication, Computing, and Measurement 2[END_REF], also in view of fault-tolerant quantum computing [START_REF] Kempe | [END_REF] or dynamic error correction [13,14,15]. Implementing quantum gates or quantum modules experimentally is in fact a challenge: one has to ght relaxation while simultaneously steering the quantum system with all its basis states into a linear image of maximal overlap with the target gate. We also showed how near time-optimal control by grape [16] take pioneering realisations from their delity-limit to the decoherence-limit [17].

While decoupling techniques typically protect against certain types of relaxation for a broad array of scenarios whose specic parameters need not be known (or are unknown anyway), optimal control exploits the knowledge of the dynamic system parameters for low-cost protection against relaxation dressed to the specics of the experimental system given. For instance, in spectroscopy, optimal control has proven to keep the state in slowly relaxing modes of the Liouville space [START_REF] Khaneja | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Xu | [END_REF]20]. In quantum computing, however, the entire basis has to be trans- * Electronic address: tosh@ch.tum.de formed. For generic relaxation scenarios, this precludes simple adaptation to the entire Liouville space: the gain of going along protected dimensions is outweighed by losses in the orthocomplement. Yet embedding logical qubits as decoherence-protected subsystem into a larger Liouville space of the encoding physical system raises questions: is the target module reachable within the protected subspace by admissible controls?

In this category of setting, the extended gradient algorithm opengrape turns out to be particularly powerful to give best approximations to unitary target gates in relaxative quantum systems thus extending the toolbox of quantum control, see e.g. [16,21,22,23,24,25,[START_REF] Ganesan | Proc. 44th. IEEE CDC-ECC[END_REF][START_REF] Sklarz | [END_REF]28,[START_REF]Introduction to Quantum Control and Dynamics[END_REF][START_REF] Machnes | [END_REF]. Moreover, building upon Ref. [START_REF] Schulte-Herbrüggen | [END_REF], it has been shown in [START_REF] Rebentrost | [END_REF] that non-Markovian relaxation models can be treated likewise, provided there is a nite-dimensional embedding such that the embedded system itself ultimately interacts with the environment in a Markovian way. Time dependent Γ(t) have recently also been treated in the Markovian [33,34] and non-Markovian regime [35].

Here we study model systems that are fully controllable [36,37,38,39], i.e. those in whichneglecting relaxation for the momentto any initial density operator ρ, the entire unitary orbit U(ρ) := {U ρ U -1 | U unitary} can be reached [37] by evolutions under the system Hamiltonian (drift) and the experimentally admissible controls. Yet, certain tasks can be performed within a subspace, e.g., a subspace protected totally or partially against relaxation.

II. THEORY

Unitary modules for quantum computation require synthesising a simultaneous linear image of all the basis states spanning the Hilbert space or subspace on which the gates shall act. Thus gate synthesis generalises the (rank-1) spectroscopic task of transferring the state of a system from a given initial one into maximal overlap with a desired target state to a full-rank unitary image.

A. Preliminaries

Yet in either case, the underlying control problem of maximising an overlap subject to the dynamics being governed by an equation of motion may readily be addressed by our algorithm grape [16].

ρ = -i [H, ρ] . (2)
In quantum computation, however, the above have to be lifted to the corresponding operator equations, which is facilitated using the notations Ad U (

•) := U (•)U † and ad H (•) := [H, (•)] obeying e -it ad H (•) = Ad U (•) with U := e -itH (3)
and using `•' for the composition of maps in

U = -iH U (4)
d dt Ad U = -i ad H • Ad U . (5)
These operator equations of motion occur in two scenarios for realising quantum gates or modules U (T ) with maximum trace delities: Let g := 1 N tr{U † target U (T )} denote the normalised overlap (setting N := 2 n for an n-qubit system henceforth). Then the quality function

f := 1 N Re tr{U † target U (T )} = Re g (6)
covers the case where overall global phases shall be respected, whereas if a global phase is immaterial [25] (while the xed phase relation between the matrix columns is kept), the quality function

f := 1 N 2 Re tr{Ad † Utarget Ad U (T ) } = g 2 (7)
applies. The latter identity is most easily seen [25] in the so-called vec-representation [START_REF] Horn | Topics in Matrix Analysis[END_REF] of ρ, where vec(ρ) is the column vector stacking all colums of the matrix ρ.

In this representation, Ad U takes the form of the conjugation superoperator Ad U = Ū ⊗ U (with Ū denoting the complex conjugate) and ad H takes the form of the commutator superoperator ad H = 1l ⊗ H -H ⊗ 1l.

B.

Open grape

Likewise, under relaxation introduced by the operator Γ (which may, e.g., take GKS-Lindblad form), the respective master equations for state transfer [36] and its ρ0 = ρSE(0) ⊗ ρB(0)

Ad W (t) ---------→ (unitary) ρ(t) = W (t)ρ0W † (t) Π SE ? ? y tr B Π SE ? ? y tr B ρSE(0) F SE (t) ---------→ (Markovian) ρSE(t) Π S ? ? y tr E Π S ? ? y tr E ρS(0) F S (t) -----------→ (non-Markovian) ρS(t)
Figure 1: Time-evolution of a quantum system (S) embedded in some environment (E) and coupled to a bath (B). If the universal system evolves under the global unitary W (t), opengrape provides optimal controls for the category of scenarios for which there is a nite-dimensional embedding such that the embedded system (SE) follows a time evolution under a Markovian quantum map FSE(t), while the reduced system of concern (S) may evolve in a non-Markovian way by a generic quantum map FS(t). In the simplest of cases, ρSE(0) is of tensor product form with FS(t) being Markovian itself.

lift for gate synthesis read

ρ = -(i ad H + Γ) ρ (8) Ḟ = -(i ad H + Γ) • F . (9)
Again with N := 2 n in n-qubit systems, F denotes a quantum map in GL(N 2 ) as linear image over all basis states of the Liouville space representing the open system. The Lie-semigroup properties of F (t) have recently been elucidated in detail [START_REF] Dirr | [END_REF]: it is important to note that only in the special (and highly unusual) case of

[ad H , Γ ] = 0 the map F (t) boils down to a mere contraction of the unitary conjugation Ad U . In general, however, one is faced with an intricate interplay of the respective coherent (i ad H ) and incoherent (Γ) part of the time evolution: it explores a much richer set of quantum maps than contractions of Ad U , as expressed in [START_REF] Dirr | [END_REF] in terms of a k, p-decomposition of the generators in gl(N 2 , C) of quantum maps. It is this interplay that ultimately entails the need for relaxation-optimised control based on the full knowledge of the master Eqn. (9), while in the special case of mere contractions of Ad U , tracking maximum qualities against xed nal times (`top curves', vide infra, e.g. Fig. 3 (a) upper panel) obtained for Γ = 0 plus an estimate on the eigenvalues of Γ suce to come up with good guesses of controls. Now for a Markovian master equation to make sense in terms of physics, it is important that the quantum subsystem of concern is itself coupled to its environment in a way justifying to neglect any memory eects. This means the characteristic time scales under which the environment correlation functions decay have to be suciently smaller than the time scale for the quantum evolution of the subsystem (see, e.g., [START_REF] Breuer | The Theory of Open Quantum Systems[END_REF]). More precisely, as exemplied in Fig. 1, we will assume that either the quantum system itself (S) or a nite-dimensional embedding of the system (SE) can be separated from the environmental bath (B) such that (at least) one of the quantum maps of the reduced system F S (t) or F SE (t) is Markovian and allows for a description by a completely positive semigroup [START_REF] Gorini | [END_REF]44,[START_REF] Davies | Quantum Theory of Open Systems[END_REF], if the time evolution for the universal composite of (embedded) system plus bath is unitary. Examples where F S (t) is Markovian have been given in [START_REF] Schulte-Herbrüggen | [END_REF], while a concrete setting of a qubit (S) coupled on a non-Markovian scale to a two-level uctuator (E), which in turn interacts in a Markovian way with a bosonic bath (B) has been described in detail in [START_REF] Rebentrost | [END_REF].

Henceforth, we will drop the subscript to the quantum map F (t) and tacitly assume we refer to the smallest embedding such that the map is Markovian and governed by Eqn. (9).

Moreover, if the Hamiltonian is composed of the drift term H d and control terms H j with piecewise constant control amplitudes u j (t k

) for t k ∈ [0, T ] H(t k ) := H d + j u j (t k )H j with u j (t k ) ∈ U ⊆ R (10)
then Eqn. ( 9) denes a bilinear control system.

With these stipulations, the grape algorithm can be lifted to the superoperator level in order to cope with open systems by numerically optimising the trace delity

f tr := 1 N 2 Re tr{Ad † Utarget F (T )} (11)
for xed nal time T . For simplicity, we henceforth assume equal time spacing ∆t :

= t k -t k-1 for all time slots k = 1, 2, . . . , M , so T = M • ∆t. Therefore F (T ) = F M • F M -1 • • • F k • • • F 2 • F 1 with every map tak- ing the form F k = exp{-(i ad H(t k ) +Γ(t k ))∆t} leads to
the following approximated derivatives

∂ftr ∂uj (t k ) -1 N 2 Re tr Ad † U •F M • F M -1 • • • F k+1 × × i ad Hj + ∂Γ(uj (t k )) ∂uj (t k ) F k ∆t × F k-1 • • • F 2 • F 1 , (12)
which hold for suciently small time intervals obeying ∆t 1/||i ad H + Γ|| 2 . Otherwise one may employ, e.g., nite dierences, since with ad H and Γ not commuting, the generators -(i ad H +Γ) are generically non-normal and thus the eigendecomposition looses its numerical benets. In the simplest case, the gradients can be used for a recursive gradient-ascent scheme (with appropriate

step size α r > 0), which gives in each time slot t k u (r+1) j (t k ) = u (r) j (t k ) + α r ∂ftr ∂uj (t k ) (13)
as update from iteration r to r + 1 of the control amplitude u j to control H j . When needed for speeding up the algorithm at higher delities, one may readily resort to standard quasi-Newton methods [START_REF] Nocedal | Numerical Optimization[END_REF] (like limitedmemory Broyden-Fletcher-Goldfarb) as shown recently [START_REF] Machnes | [END_REF].

Numerical Setting

Numerical opengrape typically started from some 50 initial conditions to each xed nal time then taking the 2-3 fastest increasing runs through some r = 10-30×10 3 iterations (see Eqn. (13)) to arrive at one point in the top curve shown as upper trace in Fig. 3.

In contrast, for nding time-optimised controls in the closed reference system, we used grape for tracking top curves: this is done by performing optimisations with xed nal time, which is then successively decreased so as to give a top curve g(T ) of quality against duration of control, a standard procedure used in, e.g., Ref. [25].

Here, in closed systems, nding controls for each xed nal time started out from some 20 random initial control sequences. Convergence to one of the points in time (where Fig. 3 shows mean and extremes for a familiy of 15 dierent such optimised control sequences) required some r = 1000 recursive iterations each. Numerical experiments were carried out on single workstations with 512 MHz to 1.2 GHz tact rates and 512 MB RAM.

Conclusions with regard to algorithmic concepts for optimising high-delity quantum maps of open systems would be premature on the grounds of this study and thus will be given elsewhere as for reliable statistics, smaller (i.e. faster to-come-by) systems will have to be addressed. and `bang-bang decoupling', we choose two model systems that can partially be tracted by algebraic means in the sense that Lie closures can readily be evaluated (vide infra ). This allows for establishing controllability [37,38,39,[START_REF] Schulte-Herbrüggen | Aspects and Prospects of High-Resolution NMR[END_REF][START_REF] Schirmer | [END_REF]49] in the subspace of low decoherence thus ensuring the existence of solutions be they by paper-and-pen or numerical. Thus comparing numerical results with analytical ones will elucidate the pros of numerical optimal control over previous approaches.

In order to avoid misunderstandings, however, we should emphasize our algorithmic approach to controlling open systems (opengrape) is by no means limited to operating within such predesigned subspaces of weak decoherence: e.g., in Ref. [START_REF] Rebentrost | [END_REF] we have worked in the full Liouville space of a non-Markovian target system. Yet, not only are subspaces of weak decoherence practically important, they also lend themselves to demonstrate the advantages of relaxation-optimised control in the case of Markovian systems (here with time in dependent relaxation operator Γ).

The starting point is the usual encoding of one logical qubit in Bell states of two physical ones

|0 L := 1 √ 2 {|01 + |10 } = |ψ + |1 L := 1 √ 2 {|01 -|10 } = |ψ - (14)
Four elements then span a Hermitian operator subspace

B := span R {|ψ ± ψ ± |} , (15)
whose protection against T 2 -type relaxation can readily be seen, since for any ρ ∈ B

Γ 0 (ρ) := [zz, [zz, ρ]] = 0 , ( 16
)
where henceforth we use the short-hand zz := σ z ⊗ σ z /2 and likewise xx as well as 1lµν1l : 16) as perfect protection against T 2 -type decoherence is in line with the slow-tumbling limit of the Bloch-Redeld relaxation by the spin tensor A 2,(0,0

= 1 2 1l 2 ⊗ σ µ ⊗ σ ν ⊗ 1l 2 for µ, ν ∈ {x, y, z, 1l}. Interpreting Eqn. (
) := 1 √ 6 ( 3 2 zz -I 1 I 2 ) [50] Γ T2 (ρ) := [A † 2,(0,0) , [A 2,(0,0) , ρ]] = 9 24 [zz, [zz, ρ]] = 0 . ( 17
)
For the sake of being more realistic, the model relaxation superoperator mimicking dipole-dipole relaxation within the two spin pairs in the sense of Bloch-Redeld theory is extended from covering solely T 2 -type decoherence to mildly including T 1 dissipation by taking (for each basis state ρ) the sum [START_REF] Ernst | Principles of Nuclear Magnetic Resonance in One and Two Dimensions[END_REF] Γ(ρ) :=

1 m1,m2=-1 A † 2,(m1,m2) , A 2,(m1,m2) , ρ , (18)
in which the zeroth-order tensor A 2,(0,0) ∼ zz is then scaled 100 times stronger than the new terms. So the resulting model relaxation rate constants nally become

T -1 2 : T -1 1 = 4
.027 s -1 : 0.024 s -1 170 : 1.

A. Controllability Combined with Protectability against Relaxation

In practical applications to a given system, a central problem boils down to simultaneously solving two questions: (i) is the (sub)system fully controllable and (ii) can the (sub)system be decoupled from fast relaxing modes while being steered to the target.

It is for answering these questions in algebraic terms that we have chosen the following coupling interactions:

if the two physical qubits are coupled by a Heisenberg-XX interaction and the controls take the form of z-pulses acting jointly on the two qubits with opposite sign, one obtains the usual fully controllable logical single qubit over B, because

i(z1l -1lz), i(xx + yy) Lie rep = su(2) , (19)
where • Lie denotes the Lie closure under commutation, which here gives i(yx -xy) as third generator to su(2). In System II the Hamiltonian adH D1+D2 comprises o-diagonal blocks (empty boxes) that make the protected modes exchange with the fast decaying ones. [NB: for pure T2-relaxation (Eqn. 17), the relaxation-rate eigenvalues would further degenerate to 0 s -1 (`decoherence-free'), 4 s -1 (medium) and 8 s -1 (fast) while maintaining the same block structure].

Model System I By coupling two of the above qubit pairs with an Ising-ZZ interaction as in Refs. [START_REF] Lidar | [END_REF]52,53] one gets the standard logical two-spin system serving as our reference System I: it is dened by the drift Hamiltonian H D1 and the control Hamiltonians H C1 , H C2 +J xyz 1lxx1l + 1lyy1l + 1lzz1l (22) and it takes the system out of the decoherence-protected subspace due to the o-diagonal blocks in Fig. 2; so the dynamics nds its Lie closure in a much larger algebra isomorphic to so [START_REF] Kempe | [END_REF], dim (iH D1+D2 ), iH C1 , iH C2 Lie = 66 , (23) to which su(4) is but a subalgebra.

H D1 :=J xx (xx1l1l + 1l1lxx + yy1l1l + 1l1lyy) + J zz 1lzz1l H C1 :=z1l1l1l -1lz1l1l H C2 :=1l1lz1l -1l1l1lz ,
Note that e -iπH Cν (H D1+D2 )e iπH Cν = H D1-D2 for either ν = 1, 2. So invoking Trotter's formula lim n→∞ e -i(H D1+D2 )/(2n) e -i(H D1-D2 )/(2n) n = e -iH D1 (24) it is easy to see that the dynamics of System II may reduce to the subspace of System I in the limit of innitely many switchings of controls H C1 or H C2 and free evolution under H D1+D2 . It is in this decoupling limit that System II encodes a fully controllable logical two-qubit system over the then dynamically protectable basis states of B ⊗ B.

In the following paragraph we may now compare the numerical results of decoherence-protection by optimal control with alternative pulse sequences derived by paper and pen exploiting the Trotter limit. As an example we choose the CNOT gate in a logical two-qubit system encoded in the protected four-qubit physical basis B ⊗ B.

B.

Results on Performing Target Operations under

Simultaneous Decoupling

The model systems are completely parameterised by their respective master equations, i.e. by putting together the Hamiltonian parts of Eqns. (20) for System I or Eqn. (22) for System II and the relaxative part expressed in Eqn. [START_REF] Khaneja | Proc. Natl. Acad. Sci. USA[END_REF]. We will thus compare dierent scenarios of approximating the logical CNOT target gate (Ad UCNOT ) by the respective quantum map F (T ) while at the same time, the logical two-qubit subsystem has to be decoupled from the fast decaying modes in order to remain within a weakly relaxing subspace. This is what makes it a demanding simultaneous optimisation task. The numerical and analytical results are summarised in Fig. 3; they come about as follows.

1.

Comparison of Relaxation-Optimised and Near

Time-Optimal Controls

With decoherence-avoiding numerically optimised controls one obtains a delity beyond 95%, while near timeoptimal controls show a broad scattering as soon as relaxation is taken into account: among the family of 15 sequences generated, serendipity may help some of them to reach a quality of 85 to 90%, while others perform as bad as giving 65%. With opengrape performing about two standard deviations better than the mean obtained without taking relaxation into account, only 2.5%

of near time-optimal control sequences would roughly be expected to reach a delity beyond 95% just by chance.

Since relaxative losses enter the quality function [START_REF] Viola | Quantum Communication, Computing, and Measurement 2[END_REF], the algorithm intrinsically nds controls entailing decoherence avoiding evolutions: Fig. 4 elucidates that the new decoherence avoiding controls keep the system almost perfectly within the slowly-relaxing subspace, whereas conventional near time-optimal controls partly sweep the system through the fast-relaxing subspace thus leading to inferior quality.

Comparison to Paper-and-Pen Solutions

Algebraic alternatives to numerical methods of optimal control exploit Trotter's formula for remaining within the slowly-relaxing subspace when realising the target, see, e.g. [START_REF] Storcz | [END_REF]. Though straightforward, they soon become unhandy as shown in Fig. 5. Assuming for the moment that to any evolution under a drift H d the inverse evolution under -H d is directly available, the corresponding naive expansions take almost 3 times the length of the numerical results, yet requiring much stronger control elds (1 -17 kHz instead of 50 Hz) as shown in Fig. 3.

In practice, however, the inverse is often not immediately reachable, but will require waiting for periodicity.

For instance, in the Trotter decomposition of Fig. 5 (c), the Ising term H ZZ := 1lzz1l as part of the drift Hamiltonian H D1+D2 is also needed with negative sign so that all terms governed by J xyz in Eqn. (22) cancel and only the Heisenberg-XX terms governed by J XX survive. But H ZZ cannot be sign-reversed directly by the z-controls in the sense 1lzz1l → -1lzz1l since it clearly commutes with the z-controls. Thus one will have to choose evolution times (τ 2 in Fig. 5) long enough to exploit (quasi) periodicity. However, H D1+D2 shows eigenvalues lacking periodicity within practical ranges altogether. Moreover, the non-zero eigenvalues of H D1+D2 do not even occur in pairs of opposite sign, hence there is no unitary transform U : H → -H = U AU † to reverse them, and a forteriori there is no local control that could do so either [56].

Yet, when shifting the coupling to J xx = 2.23 Hz to introduce a favourable quasi-periodicity, one obtains almost perfect projection (f tr ≥ 1-10 -10 ) onto the inverse drift evolution of System II, to wit U -1 := e +i π 4 H D 1 +D 2 after 3.98 sec and onto -U -1 after 1.99 sec. Thus the identity Ad (-U -1 ) = Ad (U -1 ) may be exploited to cut the duration for implementing Ad U -1 to 1.99 sec. Yet, even with these facilitations, the total length required for a realistic Trotter decomposition (with an overall trace delity of f tr ≥ 94.1 % in the absence of decoherence) amounts to some 28.5 sec as shown in Fig. 3. Moreover, as soon as one includes very mild T 1 -type processes, the (see text). For reference, the top panel of (a) shows the top curve g0(T ) consisting of maximum obtainable delities against xed nal time T in the absence of relaxation; mean and rmsd are shown for families of 15 independent control sequences generated for each T . The lower panel of (a) shows their performance vastly scattering in the presence of relaxation (•) with the intervals giving mean ± rmsd for all the 15 control sequences tested (dots for best and worst values), while numerical optimal control under explicit relaxation (•) is far superior. In (b) these results are compared to (•) representing `naive' Trotter calculations assuming to every interaction the inverse is directly obtainable and to (•) depicting a realistic Trotter approach, where the inverse has to be explicitly generated. The numbers in brackets (n1, n2, p) give the expansion coecients n1, n2 of Fig. 5 and the max. control power p (counted as number of 2π-rotations per second). Note that numerical optimal control requires some four orders of magnitude less power than the bang-bang type decoupling from fast relaxing modes used in the Trotter expansions. The upper quality limit is imposed by slow T1-type relaxation (see Eqn. [START_REF] Khaneja | Proc. Natl. Acad. Sci. USA[END_REF] in the text). Without relaxation, all the Trotter sequences would achieve delities between 93 and 99 %, except ( ) the ones limited to control elds of powers p ≤ 50 Hz: they would fall below 5%. relaxation rate constants in the decoherence-protected subspace are no longer strictly zero (as for pure T 2 -type relaxation), but cover the interval [0 s -1 , 0.060 s -1 ]. Under these realistic conditions, a Trotter expansion gives no more than 15% delity, while the new numerical methods allow for realisations beyond 95% delity in the same setting (even with the original parameter J xx = 2.0 Hz). I.

Before going into them in more detail, recall (from the section on numerical setting) that the top curve g 0 (T ) shall denote the maximum delity against nal times T as obtained for the analogous closed quantum system (i.e. setting Γ = 0) by way of numerical optimal control. Moreover, dene T * as the smallest time such that g 0 (T * ) = 1 -ε, where ε denotes some error-correction threshold.

First (I), consider the simple case of a Markovian quantum system with no encoding between logical and physical qubits, and assume g 0 (T ) has already been deter- These pulses can be seen as bang-bang type controls; due to the high repetition rates (n1 = 2 but n2 = 64) the Trotter expansions accumulate eld strengths of p 500 kHz, while the pulse shapes from optimal control require p ≤ 50 Hz. Expanding the interior Hamiltonian in (b) is even more complicated (not shown here but documented in [START_REF] Spörl | Numerische und analytische Lösungen für quanteninformatisch relevante Probleme[END_REF]).

pick the set of controls leading to g 0 (T * ) calculated in the absence of relaxation for tracking g 0 (T ). In the simplest setting, they would already be `optimal' without ever having resorted to optimising an explicitly open system. More roughly, the time-optimal controls at T = T * already provide a good approximation to ghting relaxation if T * -T * ≥ 0 is small, i.e. if γ > 0 is small. Next consider a Markovian system without coding, where Γ = γ1l is not fully degenerate (I.b). Let {γ j } denote the set of (the real parts of the) eigenvalues of Γ. Then, by convexity of {e -γt | t, γ > 0}, the following rough [62] yet useful limits to the delity f (T ) obtainable in the open system expressed by g 0 (T ) and Γ apply

g 0 (T )•exp{- T N 2 N 2 j=1 γ j } f (T ) g 0 (T )• 1 N 2 N 2 j=1
e -γj T .

(

) 25
Hence the optimisation task in the open system amounts to approximating the target unitary gate (Ad Utarget ) by the quantum map F (T ) resulting from evolution under the controls subject to the condition that modes of different decay rate constants γ j = γ k are interchanged to the least possible amount during the entire duration 0 ≤ t ≤ T of the controls. An application of this strategy known in NMR spectroscopy as TROSY [START_REF] Pervushin | Proc. Natl. Acad. Sci[END_REF] makes use of dierential line broadening [START_REF] Griey | [END_REF] and partial cancellation of relaxative contributions. Clearly, unless the eigenvalues γ j do not signicantly disperse, the advantage by optimal control under explicit relaxation will be modest, since the potential gain in this scenario relates to the variance σ 2 ({γ j }). Yet also in the case of entirely unknown relaxation characteristics, where, e.g., model building and system identication of the relaxative part is precluded or too costly, we have demonstrated that guesses of timeoptimal control sequences as obtained from the analogous closed system mayby trial and errorcope with relaxation, though mostly in a suboptimal way. This comes at the cost of making sure a suciently large family of time-optimal controls is ultimately tested in the actual experiment for selecting among many such candidates.

Clearly, this `poor man's approach' is most cumbersome thus making a strong case for taking the eort to identify a The problem actually roots in nding a viable protected subspace rather than drawing prot from it.

or at least estimate the (pertinent) system parameters.

In the non-Markovian case, however, it becomes in general very dicult to nd a common weakly relaxing subspace for encoding (II.b): there is no master equation of GKS-Lindblad form, the Γ(t) of which could serve as a guideline to nding protected subspaces. Rather, one would have to analyse the corresponding non-Markovian Kraus maps for weakly contracted subspaces allowing for encodings. However, in non-Markovian scenarios, the pros of relaxation-optimised control already become signicant without encoding as demonstrated in [START_REF] Rebentrost | [END_REF].

Simultaneous Transfer in Spectroscopy

Finally, note that the presented algorithm also solves (as a by-product) the problem of simultaneous state-tostate transfer that may be of interest in coherent spectroscopy [59]. While Eqns. ( 9) and [START_REF] Viola | Quantum Communication, Computing, and Measurement 2[END_REF] In the limit of Π being a rank-1 projector it reproduces the scenario of Eqn. [START_REF] Lidar | Decoherence-Free Subspaces and Subsystems[END_REF]. While such rank-1 problems under relaxation were explored in [START_REF] Khaneja | Proc. Natl. Acad. Sci. USA[END_REF], the algorithmic setting of opengrape put forward here allows for projectors of arbitrary rank, e.g., 1 ≤ rk Π ≤ N for n spin-1 2 qubits

with N := 2 n . Thus the rank can readily be chosen such as to equal the number of orthogonal state-to-state optimisation problems to be solved simultaneously by the same controls.

V.

CONCLUSIONS AND OUTLOOK

We have provided numerical optimal-control tools to systematically nd near optimal approximations to unitary target modules in open quantum systems. The pros of relaxation-optimised controls over time-optimised ones depend on the specic experimental scenario. We have elucidated strategies for ghting relaxation in Markovian (and non-Markovian [START_REF] Rebentrost | [END_REF]) settings with and without encoding logical qubits in protected subspaces. Numerical results have been complemented by algebraic analysis of controllability in protected subspaces under simultaneous decoupling from fast relaxing modes. The progress is quantitatively exemplied in a typical Markovian model system of four physical qubits encoding two logical ones: when the master equation is known, the new method is systematic and signicantly superior to near time-optimal realisations, which in turn are but a guess when the relaxation process cannot be quantitatively characterised. In this case, testing a set of 10 -20 such near time-optimal control sequences empirically is required for getting acceptable results with more condence, yet on the basis of trial and error. Hence determining the pertinent relaxation parameters by compressedsensing assisted process tomography will be of increasing inportance.

Complementing analytical approaches in special cases [15], optimal control tools like opengrape are therefore the method of choice in generic systems with known relaxation parameters. They accomplish decoupling from fast relaxing modes with several orders of magnitude less decoupling power than by typical bang-bang controls. In order to fully exploit the power of optimal control of open systems the challenge is shifted to (i) thoroughly understanding the relaxation mechanisms pertinent to a concrete quantum hardware architecture and (ii) being able to determine its pertinent relaxation parameters to sucient accuracy. In spite of the diculty of the latter [60],

by the recent advances in quantum state and process tomography due to techniques of comressed sensing (see, e.g., [61]), we anticipate the optimal control strategies presented here will nd wide application.

  III. EXPLORING APPLICATIONS BY MODEL SYSTEMS It is the purpose of this section is to use illustrative examples to demonstrate the power of optimal control of open quantum systems as a realistic means for protecting systems from relaxation based on the knowledge of the full master equation. In order to compare the results with idealised scenarios of `decoherence-free subspaces'

  Figure 2:(Colour) In a physical four-qubit system for encoding two logical qubits, the Hamiltonians (in their superoperator representations of adH ) take the form of 256 × 256 matrices. In the eigenbasis of Γ (Eqn. 18), the drift Hamiltonian adH D1 of System I block diagonalises into slowly relaxing modes (blue) with relaxation rate constants in the interval [0 s -1 , 0.060 s -1 ], moderately relaxing modes (hatched) with [4.01 s -1 , 4.06 s -1 ], and fast relaxing modes (red) with [8.02 s -1 , 8.06 s-1 ]. In System II the Hamiltonian adH D1+D2 comprises o-diagonal blocks (empty boxes) that make the protected modes exchange with the fast decaying ones. [NB: for pure T2-relaxation (Eqn. 17), the relaxation-rate eigenvalues would further degenerate to 0 s -1 (`decoherence-free'), 4 s -1 (medium) and 8 s -1 (fast) while maintaining the same block structure].

( 20 )

 20 where the coupling constants are set to J xx = 2 Hz and J zz = 1 Hz. Hence, over the T 2 -decoherence protected subspace spanned by the four-qubit Bell basis B ⊗ B one obtains a fully controllable logical two-qubit system iH D1 , iH C1 , iH C2 Lie B⊗B rep = su(4) .

( 21 )

 21 As illustrated in Fig.2, in the eigenbasis of Γ of Eqn.[START_REF] Khaneja | Proc. Natl. Acad. Sci. USA[END_REF] the Hamiltonian superoperators ad H take block diagonal form, where the rst block acts on the Liouville subspace B ⊗ B spanning the states protected against T 2 - type relaxation. Thus in more abstract terms and recalling Eqn. (3), the Hamiltonians of System I restricted to the T 2 -protected block, {i ad H D1 , i ad H D1 , i ad H D1 } B⊗B , generate Ad SU (4) as group of inner automorphisms over the protected states. Model System II Now, by extending the Ising-ZZ coupling between the two qubit pairs to an isotropic Heisenberg-XXX interaction, one gets what we dene as System II. Its drift term with the coupling constants being set to J xx = 2 Hz and J xyz = 1 Hz (= J zz of Eqn.(20)) reads H D1+D2 := J xx xx1l1l + 1l1lxx + yy1l1l + 1l1lyy

Figure 3 :

 3 Figure 3: (Colour) Fidelity of a CNOT gate encoded in an open system of four physical qubits in dierent scenarios of System II(see text). For reference, the top panel of (a) shows the top curve g0(T ) consisting of maximum obtainable delities against xed nal time T in the absence of relaxation; mean and rmsd are shown for families of 15 independent control sequences generated for each T . The lower panel of (a) shows their performance vastly scattering in the presence of relaxation (•) with the intervals giving mean ± rmsd for all the 15 control sequences tested (dots for best and worst values), while numerical optimal control under explicit relaxation (•) is far superior. In (b) these results are compared to (•) representing `naive' Trotter calculations assuming to every interaction the inverse is directly obtainable and to (•) depicting a realistic Trotter approach, where the inverse has to be explicitly generated. The numbers in brackets (n1, n2, p) give the expansion coecients n1, n2 of Fig.5and the max. control power p (counted as number of 2π-rotations per second). Note that numerical optimal control requires some four orders of magnitude less power than the bang-bang type decoupling from fast relaxing modes used in the Trotter expansions. The upper quality limit is imposed by slow T1-type relaxation (see Eqn.[START_REF] Khaneja | Proc. Natl. Acad. Sci. USA[END_REF] in the text). Without relaxation, all the Trotter sequences would achieve delities between 93 and 99 %, except ( ) the ones limited to control elds of powers p ≤ 50 Hz: they would fall below 5%.

Figure 4 :

 4 Figure 4: (Colour) (a) Time evolution of all the protected basis states under a typical time-optimised control of Fig.3. Projections into the slowly-relaxing and fast-relaxing parts of the Liouville space are shown. (b) Same for the new decoherenceavoiding controls. System II (see text) then stays almost entirely within the T2-protected subspace.

  extract strategies of how to ght relaxation by means of optimal control, we classify open quantum systems (i) by their dynamics being Markovian or non-Markovian and (ii) by the (Liouville) state space directly representing logical qubits either directly without encoding or indirectly with one logical qubit being encoded by several physical ones. So the subsequent discussion will lead to assigning dierent potential gains to dierent scenarios as summarised in Table

  mined. If as a trivial instance (I.a) one had a uniform decay rate constant γ so Γ = γ1l, then the delity in the presence of relaxation would simply boil down to f (T ) = g 0 (T ) • e -T •γ . Dene T * := argmax{f (T )} and

Figure 5 :

 5 Figure 5: Impractical academic alternative to numerical control: the algebraic derivation of controls for a CNOT in the slowlyrelaxing subspace proceeds from (a) the logical two-qubit system via (b) the encoded schematic physical four-qubit system to the physical realisations in the settings of (c) System I and (d) System II as dened in the text. Note that by the denition of the encoding Bell states (Eqn. (14)) all angles halve upon going from (a) to (b). Eective Hamiltonians τ H are represented by large frames; n1, n2 indicate repetitions for the respective Trotter expansions. Black bars are local π2 -pulses with phases given as subscript unless other rotation angles given on top; empty bars denote local z-rotations by ip angle π. These pulses can be seen as bang-bang type controls; due to the high repetition rates (n1 = 2 but n2 = 64) the Trotter expansions accumulate eld strengths of p 500 kHz, while the pulse shapes from optimal control require p ≤ 50 Hz. Expanding the interior Hamiltonian in (b) is even more complicated (not shown here but documented in[START_REF] Spörl | Numerische und analytische Lösungen für quanteninformatisch relevante Probleme[END_REF]).

  The situation becomes signicantly more rewarding when moving to the category (II) of optimisations restricted to a weakly relaxing (physical) subspace used to encode logical qubits. A focus of this work has been on showing that for Markovian systems encoding logical qubits, the knowledge of the relaxation parameters translates into signicant advantages of relaxation-optimised controls over time-optimised ones. This is due to a dual eect: opengrape readily decouples the encoding subsystem from fast relaxing modes while simultaneously generating a quantum map of (close to) best match to the target unitary. Clearly, the more the decay of the subspace diers from its embedding, the larger the advantage of relaxation-optimised control becomes. Moreover, as soon as the relaxation-rate constants of the protected subsystem also disperse among themselves, modes of dierent decay should again only be interchanged to the least amount necessarythus elucidating the very intricate interplay of simultaneous optimisation tasks that makes them prone for numerical strategies.

  refer to the fullrank linear image F , one may readily project onto the states of concern by the appropriate projector Π to obtain the respective dynamics and quality factors of the subsystem Π Ḟ = -Π (i ad H + Γ) • F Re tr{Π F † target Π F (T ) } . (27)

Table I :

 I Gain Potential for Relaxation-Optimised Controls versus Time-Optimised Controls

	Category	Markovian	non-Markovian
	encoding: protected subspace	big	(dicult a )
	no encoding:
	full Liouville space	smallmedium mediumbig [32]

Acknowledgments

This work was supported by the integrated eu projects qap and q-essence as well as by Deutsche Forschungsgemeinschaft, dfg, in sfb 631 and by the Bavarian excellence network qccc.


B Niermann

A Kanitz

M Böke

J Winter

Impurity intrusion in radio-frequency micro-plasma jets operated in ambient air

Keywords: numbers: 52.25.Vy, 52.38.Dx, 52.25.Ya

Space and time resolved concentrations of helium He ( 3 S 1 ) metastable atoms in an atmospheric pressure radio-frequency micro-plasma jet were measured using tunable diode laser absorption spectroscopy. Spatial profiles as well as lifetime measurements show significant influences of air entering the discharge from the front nozzle and of impurities originating from the gas supply system. Quenching of metastables was used to deduce quantitative concentrations of intruding impurities. The impurity profile along the jet axis was determined from optical emission spectroscopy as well as their dependance on the feed gas flow through the jet.

Introduction

In recent years the research on atmospheric pressure micro-plasmas became a strong focus in plasma sciences, mostly due to their high potential for new plasma applications without the need for expensive vacuum equipment [1]. Among the large variety of microplasma sources, that make use of DC, pulsed DC and AC ranging from mains frequency to RF, atmospheric pressure plasma jets provide a simple design by featuring an α-mode RF discharge between two bare metallic electrodes.

Jet discharges are particularly well suited for applications in ambient atmosphere, like biological, medical, or a variety of surface treatment applications. Operating the discharge in direct contact with air, however, makes the plasma vulnerable to intrusions of nitrogen and oxygen molecules from the surrounding environment. This may reduce the operational space, induce instabilities and create unwanted reactive species. To study these influences in detail, excited species in the discharge provide an indirect measure for the amount of impurities in the discharge, since they react very sensitively to collisions with these species. Especially metastable species carry large amounts of potential energy and are a source of ionization. In this context helium metastable species pose an appropriate indicator, since helium is a widely used carrier gas in atmospheric pressure micro-discharges. Compared to molecular gases, the excitation threshold of helium metastable atoms is high, and can be exceeded just by high energetic electrons in the tail of the EEDF. Thereby the metastable excitation process is highly sensitive to impurities in the discharge. Due to their long lifetime they collide more frequently with other particles (about 15000 times per µs at atmospheric pressure), which makes them important for plasma chemistry processes and very sensitive to even tiny amount of molecular gases in the discharge. Especially nitrogen and oxygen molecules have some of the largest quenching cross sections for metastable helium atoms [2]. The metastable density in micro-discharges is several orders of magnitude lower than the density of the ground-state atoms. However, compared to most other species the density is significant and the electron collision excitation cross sections of some helium levels out of the metastable states exhibit values which are several orders of magnitude larger and have much lower thresholds than those for the ground state [3,4]. Thus, metastables strongly contribute to the ionization.

A reliable technique for the systematic investigation of metastable species is tunable diode laser absorption spectroscopy (TDLAS), since it is not invasive and provides absolute population densities. We applied TDLAS to record the spectral profiles and the lifetime of the lowest helium metastable state, deducing absolute densities with high sensitivity and spatial resolution, revealing detailed information about the distribution of impurities in the discharge. In addition to TDLAS we applied optical emission spectroscopy (OES) to observe nitrogen emissions lines that provide direct information about the nitrogen content in the discharge.

The atmospheric pressure micro-plasma jet

The atmospheric pressure micro-plasma jet is a capacitively coupled, non-thermal glowdischarge at high pressures. The design concept of this discharge is based on the plasma jet introduced by Selwyn et al. in 1998 [5] and advanced by Schulz-von der Gathen et al. [6]. Feed gas flows between two closely spaced stainless steel electrodes driven at 13.56 MHz radio-frequency in a parallel plate configuration (Fig. 1). Electrodes and plasma volume are enclosed by quartz windows, giving direct optical access to the discharge. The jet is operated in helium at various flow rates from 200 sccm to 5 slm. The electric field between the electrodes causes a breakdown in the gas and produces a plasma with electron temperatures and densities of about 1 to 2 eV and 10 10 cm -3 , respectively [7,8]. Atoms and molecules in the feed gas become excited, dissociated or ionized by electron impacts. Since the electrons are not in thermal equilibrium with the ions and neutrals, the gas temperature remains a few tens of K above room temperature [9]. The distance between the two windows is fixed to 1 mm while the electrode gap size is variable between 0.2 and 3 mm. In all here presented configurations the discharge operates as a typical α-mode rf glow discharge.

The whole jet setup is located in an airtight stainless steel vessel, whose atmosphere can be defined accurately. Thus, the jet can either be operated in ambient air, or in a pure helium environment, and a precise retracement of leak channels is possible.

Spectroscopic setup

The small dimensions of micro-discharges and their operation at atmospheric pressures are a challenge for optical diagnostics, since high sensitivity and high spatial resolution are required. For the absorption spectroscopic measurements a standard TDLAS set-up was used. The absorption profile was recorded by scanning the laser frequency across the absorption line. Figure 2 shows a sketch of the experimental setup. The laser beam from the DL passed through two beam splitters. A part of the beam was guided to a Fabry-Perot interferometer (1 GHz free spectral range), a second part through a low pressure reference cell to perform the calibration of the laser frequency. The part of the beam transmitted through the first beam splitter was attenuated by neutral density filters with an optical density in the order of 3, and focused into the discharge with a beam power of less than 2 µW at 100 µm spot size, to avoid any saturation effects. After passing the discharge, the beam was guided through a set of apertures and filters to suppress the emission from the plasma by reducing the collection angle and blocking wavelengths different than the observed transitions. The transmitted beam intensity was measured either by a photodiode with on-chip transimpedance amplifier for accurate density measurements, or with a biased high speed photodiode for lifetime measurements. For an effective measurement of the absorption signal across the jet axes, the discharge casing was mounted on a small movable stage featuring three electronically controlled stepping motors to adjust, with high precision, the discharge cell position in all spatial dimensions. This set-up allows the positioning of the jet with an accuracy of about 5 µm and automated xz-mapping of the complete plasma volume. In y-axis all measurements presented in this paper are line averaged.

OES measurements were carried out using a lens system to image the plasma emissions into a fiber, guiding the signal to a spectrometer with a spectral resolution of 0.75 nm. Emission profiles are averaged in y and z dimensions. The spatial resolution in x is 1 mm.

Since the absorption rate of the laser light by metastable atoms is very low, in the order of 10 -3 after 1 mm absorption length, lock-in technique was used to measure the changes in signal intensity. Applying lock-in technique requires the pulsing of the signal to be measured. This was realized by pulsing the RF-power coupled into the system, which consequently leads to a pulsing of the metastable density in the discharge. The pulse frequency was chosen to 4 kHz and a duty cycle of 50 %.

Absolute metastable densities were derived from the transmittance I I 0 and the Beer-Lambert law. Therefore four signals have to be acquired: L(ν) -Plasma and Laser on, L 0 (ν) -Plasma off, Laser on, P (ν) -Plasma on, Laser off, B(ν) -Plasma and Laser off (background), to calculate the transmittance spectra and correlate them with the plasma properties by

I(ν) I 0 (ν) = L(ν) -P (ν) L 0 (ν) -B(ν) = e -k(ν)•l , (1)
where I(ν) and I 0 (ν) are the intensities of transmitted radiation with and without the presence of absorbing species, k(ν) is the absorption coefficient and l the path length through the absorbing medium [10]. The absorption coefficient is connected to the population density of metastable atoms by

k(ν) = 1 4π 0 • π(e 2 ) cm e • f ik • N i • F (ν), (2)
where f ik is the oscillator strength of the line, N i the density of the lower level, and

F (ν) a normalized function ( ∞ 0 F (ν) • dν = 1)
representing the absorption lineshape. All other terms have their usual definitions. The absolute metastable density can then be given by

∞ 0 ln I 0 (ν) I(ν) dν = S = e 2 f ik l 4 0 m e c • N i , (3)
where S is the area under the absorption curve that provides the line-averaged density of the absorbing species.

Results and discussion

Metastable quenching by air

The data shown in figure 3 (upper graph) represent the metastable density in dependence on the feed gas flow. They have been taken in the center of the discharge (2 cm from the exit nozzle) while the jet was running surrounded by ambient air. The density increases with the gas flow rate and shows a weak tendency to saturate at higher flows. This strong dependence can not be explained by a fundamental change in the excitation mechanism, since neither the pressure nor the outer power coupling to the discharge is changing. To determine the origin of this mechanism in detail the decay rate of the density in the afterglow was measured under different conditions. Figure 3 (bottom graph) indicates that the lifetime of the species increases also with the gas flow and, like the density, shows a tendency to saturate for higher flows. Assuming infinite purity of the helium gas and proposing that the metastable lifetime is mainly determined by the three-body collision process:

He * + 2He → He * 2 + He, (4)
the theoretical lifetime of He* ( 3 S 1 ) in pure helium is expected to be

τ = (K 3 • N 2 He ) -1 = 5.8 µs, (5)
with K 3 = 2.5 • 10 -34 cm 6 s -1 the rate coefficient for three-body collisions and N He = 2.68 • 10 19 cm -3 the helium ground state density [12,13]. The loss of metastable atoms by two-body collisions as well as by diffusion can be omitted since the diffusion loss frequency is 2 orders of magnitude lower then the loss frequency by three-body collisions. Three situations have to be distinguished in this context. In the first setup (red triangles in figure 3) the airtight vessel around the jet was evacuated and filled with pure helium. Furthermore the gas supply system, that is the stainless steel tube connection from the gas bottle to the micro-plasma jet, was kept under a vacuum of 1 Pa for some time, to clean it from traces of water and air. As a result for high flows of 5 slm the theoretical lifetime for a clean helium feed gas can be reached, but the lifetime is dropping significantly for lower flows. Taking into account that the leakage rate of the gas supply system can be neglected, the metastable quenching can be attributed to tiny amounts of water and air molecules desorbing from the surface and contaminating the feed gas in the order of a 2 to 40 ppm. This fraction of contamination decreases linearly with the feed gas flow. In a second setup (pink triangles in figure 3) the airtight vessel was filled with ambient air. The measured lifetimes in this case are significantly below the theoretical value for a clean feed gas. Again the discrepancy can be attributed to the loss of metastable atoms by molecular excitation processes. Due to the low ionization thresholds these excitation processes are mostly Penning ionizations, in this case, with air introduced through the front nozzle:

He * + N 2 → He + N + 2 + e, (6)
He * + O 2 → He + O + 2 + e. (7)
Assuming that the dominant impurity contribution is due to residual N 2 and O 2 molecules, the impurity level can be estimated by a simple model. Taking into account quenching collisions with nitrogen and oxygen molecules the lifetime is limited to

τ = [K 3 • N 2 He + K N 2 • (N N 2 ,Bottle + N N 2 (Γ)) (8) + K O 2 • (N O 2 ,Bottle + N O 2 (Γ))] -1 , (9)
with K N 2 = 7 • 10 -11 cm 3 s -1 and K O 2 = 2.5 • 10 -10 cm 3 s -1 the rate constants for Penning ionization processes with molecular nitrogen and oxygen and N N 2 (Γ) as well as N O 2 (Γ) the flow dependant nitrogen and oxygen densities [13,14]. N N 2 ,Bottle and N O 2 ,Bottle are the base contaminations of the feed gas, given by the purity of the helium used in the experiment. This base contamination sums up to about 2 ppm. Since we are assuming the intrusion of air into the system, the ratio between N 2 and O 2 is fixed to 3.7. Furthermore we presume the intrusion air to vary anti-proportionally with the flow rate. The fitting parameter is the amount of air in the system and the fitting was done until the modeled data optimally reproduced the measured data. The result of the model is shown in figure 3 (blue line). It agrees well with the measured values over the whole flow rate range. We determine the level of air introduced through the nozzle to be in the order of a few dozen ppm (120 ppm for 1 slm flow, 25 ppm for 5 slm flow).

In a third setup (green triangles in figure 3) the micro-plasma jet as well as the gas supply system was exposed to ambient air. Repeating the flow dependant lifetime measurements the values show again a significant decrease. Nevertheless, the amount of air introduced into the system through the front nozzle must be the same like in the previous case. Additional metastable quenching can only be caused by traces of air and water desorbing from the surfaces of the gas supply system as well as from the metallic electrodes. This assumption is supported by mass spectrometry and OES measurements that show significant amount of water ions as well as nitrogen and oxygen emissions in the plasma and the effluent region. Matching these measurements with the previously described model (orange line) the observed quenching corresponds to a few hundred ppm of air molecules in the feed gas stream (500 ppm for 1 slm flow, 100 ppm for 5 slm flow).

For metastable density and lifetime measurements no error bars are shown, since statistical as well as systematic errors are negligible. For the presented data they are in the order of 10 9 cm -3 Due to the use of lock-in technique a sensitive measurement of low densities was possible with a high signal to noise ratio. A minor systematic error is caused by the unknown metastable profile in y-axis, and the implications on the absorption length. But since the metastable production and loss processes take place on a very local scale at atmospheric pressure, the profile in y-axis is expected to be homogeneous with the exception of a tiny volume in front of the glass windows that confines the plasma channel and acts as a sink for metastable atoms. Lifetime measurements were performed with a time resolution of less then 30 ns and revealed a perfectly exponential decay of the species. Thus errors can be omitted in this case as well.

Longitudinal impurity profile

Figure 4 shows a 2D-map of the He ( 3 S 1 ) metastable density in the discharge volume. Both, horizontal and vertical axis show the exact area between the electrodes. The map covers 2.000 reading points (40 vertical x 50 horizontal) of the absorption signal in the plasma volume. A variety of effects determine the metastable distribution in the vertical and in the longitudinal axis. While the vertical profile reflects the RF-sheaths structure and is determined by the electron density and temperature distribution, the longitudinal profile is governed mainly by environmental factors [11].

As discussed before the impurities entering the jet through the front nozzle are quenching the metastable atoms. Since the jet is running permanently in contact with the ambient air the intrusion of nitrogen and oxygen into the plasma channel is significantly high, and decreases the metastable density especially in the first millimeters from the nozzle. This observation is supported by OES measurements showing the emission of various nitrogen lines. Figure 5 presents the longitudinal profile of the molecular nitrogen emissions at 357 nm, representing the second positive system. The measurements reveal strong nitrogen emissions in the first 10 mm from the front nozzle. The intrusion depth increases with decreasing gas velocities, ranging from about 10 to 100 ms -1 . Especially for small gas flows a back diffusion of nitrogen can be observed over the whole length of the jet. The back diffusion of air into the system is promoted by the geometry of the gas channel, whose shape is rectangular. Especially in the corners of the channel the gas flow is likely not to be laminar and turbulences may occur.

OES measurements have been performed with a very limited spectral and spatial resolution, as described in the previous section. Therefore the expected errors are much larger than for TDLAS measurements. The signal to noise ratio strongly depends on the amount of nitrogen in the discharge. In the present case the errors for lower flows (around 0.5 slm) and therefore high nitrogen content are expected to be in the order of 10 %, while for higher flows (around 5 slm) the errors are much larger, in the order of 50 %.

Summary

It was shown that in atmospheric pressure micro-plasma jets operated in contact with ambient air the back diffusion of nitrogen and oxygen into the gas channel is a significant effect with strong impact on the energy transfer processes of the discharge. Furthermore the inner surfaces of the gas supply system and of the jet discharge itself pose a source of impurities contaminating the the feed gas. These influences and the likely consequences on the plasma chemistry have to be taken into account for any application that envisage the use of such discharges in ambient atmosphere.

Figure 1 .

 1 Figure1. Sketch of the micro-plasma jet discharge. Shown is a 2-dimensional crosssection through plane that is spanned by the two electrodes.

Figure 2 .

 2 Figure 2. The TDLAS setup.

Figure 3 .

 3 Figure 3. He ( 3 S 1 ) metastable density (top) and lifetime (bottom) for different setups in dependence on the feed gas flow. Also shown is the modeled lifetime taking into account significant contribution of nitrogen and oxygen quenching.

Figure 4 .

 4 Figure 4. 2-dimensional map of the measured metastable density in the discharge volume for the He 2 3 S 1 → 2 3 P 0 1,2 transition. Powered electrode at the top, grounded electrode at the bottom. Electrode gap size was 1.8 mm. Densities are given in cm -3 . Measurements were taken in under ambient atmosphere.

Figure 5 .

 5 Figure 5. Longitudinal emission profile of the second positive system of molecular nitrogen at 357 nm. Measurements were taken under ambient atmosphere.

Acknowledgments

This project is supported by DFG (German Research Foundation) within the framework of the Research Group FOR1123 and the Research Department 'Plasmas with Complex Interactions' at Ruhr-University Bochum.


M Lapert

Y Zhang

S J Glaser

D Sugny

Towards the time-optimal control of dissipative spin 1/2 particules in Nuclear Magnetic Resonance

Keywords: numbers: 32.80.Qk, 03.65.Yz, 78.20.Bh

We consider the time-optimal control of a spin 1/2 particle whose dynamics is governed by the Bloch equations with both longitudinal and transverse relaxation terms. We solve this control problem by using geometric optimal control techniques. We show the crucial role of singular extremals in the time-optimal synthesis. This role can mainly be attributed to the presence of dissipation. We also analyze the robustness of the optimal control sequence when both the maximum amplitude of the control field and the dissipative parameters are varied. Finally, we present an experimental implementation of the different solutions using techniques of Nuclear Magnetic Resonance.

Introduction

Nuclear Magnetic Resonance (NMR) is one of the most promising fields of applications of quantum control [START_REF] Rice | Optimal control of quantum dynamics[END_REF][START_REF] Shapiro | Principles of quantum control of molecular processes[END_REF][START_REF] Tannor | Introduction to quantum mechanics: A time-dependent perspective[END_REF]. In this domain, optimal control techniques can be used to design magnetic fields to control the dynamics of spin systems with applications extending from quantum computing to spectroscopy [START_REF] Ernst | Principles of Nuclear Magnetic Resonance in one and two dimensions[END_REF][START_REF] Levitt | Spin dynamics: basics of nuclear magnetic resonance[END_REF]. In this framework, numerical optimization procedures such as the GRAPE algorithm (Gradient Ascent Pulse Engineering Algorithm) have been developed [START_REF] Khaneja | [END_REF]7,8]. Recently, methods of geometric optimal control theory have also been applied with success [9,10,11,12,13]. This approach which is based on strong mathematical tools coming from differential geometry and Hamiltonian dynamics has a rapid development permitting to attack problems of increasing difficulty [START_REF] Jurdjevic | Geometric control theory[END_REF][START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Boscain | Optimal syntheses for control systems on 2-D manifolds[END_REF]. In this context, we have derived in Refs. [17,18] the complete solution of the optimal control of two-level dissipative quantum systems whose dynamics is governed by the Lindblad equation [19,20]. We have applied this analysis in Ref. [10] to the time-optimal control of a spin 1/2 particle in a dissipative environment. We have recently extended this work to the energy minimization control problem of a spin [22] and to the control of a spin in presence of both relaxation and radiation damping effects [23]. It is these different results that we propose to review in this paper. In particular, we will show the essential role of singular extremals in the optimal control law. Note that the existence of such singular solutions has been essentially ignored so far in the quantum control literature, and only few results exist [12,21]. Some parts of the material of this work overlap with Refs. [10,22,23], although the presentation is different. We will extend some of these results by constructing the optimal synthesis (i.e. the set of all the optimal solutions starting from the same initial point and reaching any point of the accessible set) and by analyzing the robustness properties of the optimal sequence with respect to variations of the amplitude of the control and of the dissipative parameters. Finally, we refer the reader to the works [17] for a more detailed background on geometric control and proofs of the results.

In this article, we consider the control of a spin 1/2 particle whose dynamics is governed by the Bloch equations with both longitudinal and transverse relaxation terms. As an example, we analyze the saturation problem in minimum time which consists in bringing the magnetization vector from the north pole of the Bloch sphere to the center of the Bloch ball. By its simplicity, this problem can be viewed as a standard textbook problem in quantum control. It is also one of the simplest but non trivial system where all the powerful machinery of geometric control techniques can be applied. The resolution of this control problem will be used to explain some basic tools of geometric optimal control theory. This will help the reader to enter into more specialized mathematical books of optimal control [START_REF] Jurdjevic | Geometric control theory[END_REF][START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Boscain | Optimal syntheses for control systems on 2-D manifolds[END_REF].

The paper is organized as follows. We introduce in Sec. 2 the model system used in the computation. Section 3 focuses on the Pontryagin Maximum Principle (PMP) [START_REF] Jurdjevic | Geometric control theory[END_REF][START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Boscain | Optimal syntheses for control systems on 2-D manifolds[END_REF] which is the main tool used to solve the control problem. We present some geometric techniques to construct the optimal solution. To simplify the discussion, we only consider the control of a spin 1/2 particle by one control field. The geometric optimal solutions for the saturation problem are presented in Sec. 4. Both the bounded and the unbounded cases are analyzed. An example of experimental implementation using techniques of Nuclear Magnetic Resonance (NMR) is proposed in Sec. 5. A conclusion and prospective views are given in a final section.

The model system

We consider a homogeneous ensemble of uncoupled spins 1/2 which are irradiated on resonance by a magnetic field. Homogeneous means here that the resonance offset (or detuning) of each spin is the same. In the rotating frame, the equation of motion can be written as follows [START_REF] Levitt | Spin dynamics: basics of nuclear magnetic resonance[END_REF]:

    Ṁx Ṁy Ṁz     =     -M x /T 2 -M y /T 2 (M 0 -M z )/T 1     +     ω y M z -ω x M z ω x M y -ω y M x     , (1)
where the vector of components (M x , M y , M z ) is the magnetization vector and M 0 the equilibrium point of the dynamics along the zaxis. The dynamics depends on two different terms: A drift term due to the longitudinal and transversal relaxation rates

T -1 1 and T -1 2
and a second term representing the effect of the control amplitudes ω x and ω y . We assume that the control field ω = (ω x , ω y , 0) satisfies the constraint | ω| ≤ ω max . We introduce the normalized coordinates x = (x, y, z) = M/M 0 , which implies that at thermal equilibrium the z component of the scaled vector x is by definition +1. The normalized control field which satisfies |u| ≤ 2π is defined as u = (u x , u y , 0) = 2π ω/ω max , while the normalized time τ is given by τ = (ω max /2π)t. Dividing the previous system by ω max M 0 /(2π), one deduces that the dynamics of the normalized coordinates is ruled by the following system of differential equations:

    ẋ ẏ ż     =     -Γx -Γy γ -γz     +     u y z -u x z u x y -u y x     ,
where Γ = 2π/(ω max T 2 ) and γ = 2π/(ω max T 1 ). We analyze a control problem where the initial point of the dynamics is the equilibrium point, i.e. the north pole of the Bloch sphere. The goal of the control will be to reach any point of the corresponding accessible set and in particular the center of the Bloch ball. In the setting of NMR spectroscopy and imaging, this latter question corresponds to saturating the signal, e.g., for solvent suppression or contrast enhancement, respectively [24,25]. Since the initial point is on the zaxis, the control problem admits a symmetry of revolution around this axis. Without loss of generality, it can be shown that one of the components of the control field can be taken to zero, e.g. here ω y [17]. In the reduced coordinates, this leads to u y = 0. We are thus considering a single input problem in a plane of the form:

ẏ ż = -Γy γ -γz + u -z y ,
where the subscript x has been omitted for the control parameter. This differential system can be written in a more compact form as follows: ˙

X = F ( X) + u G( X), (2)
where the coordinates of the state X of the system are (y, z) and F and G are two vector fields of components (-Γy, γ -γz) and (-z, y).

The Pontryagin Maximum Principle

We use the PMP to solve the optimal control problem with the constraint of minimizing the control duration [START_REF] Boscain | Optimal syntheses for control systems on 2-D manifolds[END_REF]. In addition, the field has to satisfy the relation |u| ≤ 2π. The PMP is formulated from the pseudo-Hamiltonian H which is given in our case by:

H = P • ( F ( X) + u G( X)), (3)
where we have introduced the adjoint state P of components (p y , p z ). The PMP tells us that the optimal solutions are a subset of the set of extremals which are Hamiltonian trajectories of H with

˙ X = ∂H ∂ P , ˙ P = - ∂H ∂ X . (4)
The optimal control law is given by the maximization condition

H( X, P , v) = max |u|≤2π H( X, P , u), ( 5
)
where v is the control field maximizing H and P a nonzero vector. In addition, since the Hamiltonian H does not depend on time, one deduces that it is a constant of the motion which is positive by construction of a maximum principle, H being negative in the time-maximum case and zero in the abnormal one [START_REF] Boscain | Optimal syntheses for control systems on 2-D manifolds[END_REF]. The Hamiltonian equations can be explicitly written in the case of Eq. ( 2) as ˙ X = F ( X) + u G( X) and

˙ P = -P • ( ∂ F ∂ X + u ∂ G ∂ X ), (6)
where ∂/∂ X denotes the gradient. ∂ F /∂ X is a 2 × 2-matrix whose element of the iline and jcolumn is equal to ∂F i /∂X j . The construction of the optimal solution requires the introduction of two set of points, the singular and the collinear sets which are respectively denoted by S and C. The set C is defined by the points X where the two vectors F and G are collinear. This condition is satisfied if the determinant denoted by det(•, •) of the two vectors is zero. We recall that

det( F , G) = F 1 G 2 -G 1 F 2 , (7)
where the F i and the G i are the components of F and G. Here this leads to -Γy 2 + γ(z -z 2 ) = 0. The set C is therefore the union of two parabolas if γ, Γ = 0. The definition of the set S is more involved since it is given by the relation det( G, [ F , G]) = 0, [ F , G] being the commutator of two vector fields which can be computed as follows in the coordinates x i :

[ F , G] j = i ( ∂G j ∂x i F i - ∂F j ∂x i G i ).
In our example, one arrives at 2Γyz + γ(y -2yz) = 0.

One deduces that S corresponds to the union of the vertical line y = 0 and of the horizontal one z = z 0 given by

z 0 = - γ 2(Γ -γ) = - T 2 2(T 1 -T 2 ) if Γ = γ (or equivalently if T 1 = T 2 ).
In the case Γ ≤ 3 2 γ then only the vertical line exists. The two sets S and C are displayed in Fig. 1 for a given set of parameters.

Returning back to the PMP, one sees that the maximization condition can be solved by introducing the switching function Φ = P • G = -p y z + p z y [START_REF] Boscain | Optimal syntheses for control systems on 2-D manifolds[END_REF]. Using Eq. ( 5) and the fact that the drift term P • F of H cannot be controlled, one immediately deduces that the control field u is given by

u = 2π × sign[Φ(t)]
if Φ(t) = 0. Such a point is said to be regular. The corresponding control of constant amplitude is called bang. A time such that the control changes sign is a switching time. Such a case occurs when Φ vanishes in an isolated point t = t 0 . This is associated to a bang-bang control sequence. The singular situation corresponds to the case where Φ is equal to zero on an interval [t 0 , t 1 ]. At this point, a crucial remark is the fact that the singular trajectories lie in the set S. This can be shown by using the relation Φ(t) = Φ(t) = 0 on the interval [t 0 , t 1 ]. The first derivative of Φ can de determined as follows. We have:

Φ = ˙ P • G + P • ˙ G.
Using Eq. ( 4) and

˙ G = ∂ G ∂ X • ˙ X, one gets that: Φ = P • [ F , G].
The system Φ(t) = Φ(t) = 0 admits a non zero P solution if the vectors [ F , G] and G are parallel, i.e. if det( G, [ F , G]) = 0. This first computation does not give the explicit form of the singular control u s , which is obtained from the second derivative of Φ. Indeed, a similar computation to the one for Φ leads to:

Φ = P • [[ G, F ], F ] + u P • [[ G, F ], G],
Towards the time-optimal control of dissipative spin 1/2 particules in Nuclear Magnetic Resonance6

and we get that

u s = - P • [[ G, F ], F ] P • [[ G, F ], G] if P • [[ G, F ], G] = 0.
If this latter relation is not satisfied, derivatives of higher orders have to be computed to determine u s .

In our case, one arrives at:

u s (y, z) = -yγ(Γ -2γ) -2yz(γ 2 -Γ 2 ) 2(Γ -γ)(y 2 -z 2 ) -γz , (8)
which simplifies into u s = 0 for the vertical singular line and into

u s = γ(γ -2Γ) 2(Γ -γ)y
for the horizontal one where z = z 0 . It can also be checked that the denominator is different from zero on the set S, except at the intersection of the two singular lines. Another important point is the admissibility of the singular control since by definition of the control problem we must have |u s | ≤ 2π. Using Eq. ( 8), one deduces that the singular control is admissible on the horizontal line if

|y| ≥ |γ(γ -2Γ)|/[2π(2Γ -2γ)].
For smaller values of y, the system cannot follow the horizontal singular arc and must leave this line by using a bang pulse. A consequence of this phenomenon, which is called a saturation of the control field, is the creation of a switching curve starting from the point where the field losses its admissibility [START_REF] Boscain | Optimal syntheses for control systems on 2-D manifolds[END_REF]. A switching line is a line such that the optimal control changes sign when crossing it. This curve can be constructed numerically as explained in Sec. [START_REF] Ernst | Principles of Nuclear Magnetic Resonance in one and two dimensions[END_REF].

From this first analysis, we know that the extremal trajectories will be the concatenation of bang arcs of maximum or minimum amplitude and of singular arcs on the set S. However, this study tells us nothing about the optimality status of the extremals. This optimality can be in some cases determined from the clock form.

The clock form

We derive in this paragraph the expression of the clock form denoted α [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]. By definition, the clock form is a 1-form, i.e. a linear functional on the space of vector fields, which fulfills the following conditions

α( F ) = 1 α( G) = 0 .
A solution of this system exists except on the set C where F and G are collinear. Indeed, assuming that F = β G where β is a scalar function, one then deduces that

1 = α( F ) = βα( G) = 0
which gives a contradiction. We write α as α = α y dy + α z dz where dy and dz are the differentials of the coordinates y and z. We recall that dy is a mapping which associates to any vector field its ycoordinate. Simple algebra shows that α y and α z are solutions of the system

α y (-Γy) + α z (γ -γz) = 1 α y z = α z y .
We obtain that

   α y = -y Γy 2 -γz+γz 2 α z = -z Γy 2 -γz+γz 2 .
From the 1-form α, we can define the 2-form dα which is given by

dα = ( ∂α z ∂y - ∂α y ∂z )dy ∧ dz,
and reads after some calculations

dα = 2Γyz + γy -2γyz [Γy 2 -γz + γz 2 ] 2 dy ∧ dz.
If we write dα as dα = g(y, z)dy ∧ dz then one sees that g(y, z) = 0 on S and that g is infinite on C. It is then straightforward to deduce that the function g has a constant sign in the regions delimited by the lines of S. This point is represented on Fig. 1 for a given set of parameters. Four quadrants with different signs of g are displayed. As suggested by its name, the clock form is a form which allows one to determine the time taken to travel a path and to compare the extremals. Let γ be a path in the plane (y, z) which does not cross C and T the time of travel along γ. We have

γ α = T 0 α( ˙ X)dt = T 0 α( F )dt = T.
We consider now two paths γ 1 and γ 2 starting and ending at the same points and respectively associated to the durations T 1 and T 2 (see Fig. 1 for an illustrative example). One shows by using the Stokes theorem that

T 1 -T 2 = γ 1 α - γ 2 α = D dα,
where D is the surface delimited by γ 1 ∪ -γ 2 . For paths γ 1 and γ 2 which lie in one of the four quadrants defined by S, we can straightforwardly determine the time-optimal trajectory. For instance on the Fig. 1, the time to travel the path γ 1 is lower than the one to travel γ 2 since the two paths belong to a region with g < 0.

Using such arguments, we can also analyze the optimality of singular lines. The proof uses a singular path and a path starting and ending on the singular line but associated to a bang-bang control. Note that this result is only a local optimality result in the sense that the singular trajectory is optimal for any neighboring initial and final points on this axis. The difference between local and global optimality will be detailed in Sec. 4. In the example of Sec. 4, it can be shown that the horizontal singular line is optimal and that the vertical one is optimal only if z > z 0 .

Time-optimal control of spin 1/2 particles

We apply in this section the tools of Sec. 3 to the control of spin 1/2 particles. We begin our study by considering the control problems defined by the relaxation parameters γ -1 and Γ -1 (expressed in the normalized time unit defined above) of 23.9 and 1.94, respectively and M 0 ≈ 2.15 × 10 -5 . Such values of the parameters correspond to a realistic experimental situation in NMR where T 1 = 740 ms, T 2 = 60 ms and ω max /(2π) = 32.3 Hz [10].

Optimal synthesis

This paragraph is rather technical and describes the way to construct the optimal synthesis, which is the set of all the optimal sequences starting from the north pole of the dynamics and reaching any point of the corresponding accessible set A. Since the control field cannot fully compensate the effect of dissipation [26,27], the system is not controllable and all the points of the Bloch ball cannot be reached. In this case, the reachable set can be determined by the explicit construction of all the trajectories satisfying the constraint |u| ≤ 2π (see Fig. 2 for an example). More precisely, we first search for the boundary of this set, i.e. the trajectories with u = ±2π originating from the north pole. Then we show that all the points inside the boundary are attainable by the explicit construction of the optimal control laws as displayed in Fig. 2.

The application of the PMP gives only the local behavior of the extremals. In particular, we know that the optimal solutions are the concatenations of bang arcs of amplitude ±2π and of singular optimal arcs. Different extremal solutions which are locally optimal can be used to reach the same target state. The last step of the optimal synthesis consists of choosing among these extremals the time-minimizing one, which corresponds to the optimal one. This choice can be done by using the clock form and numerical comparisons of extremals when the paths cross the singular or the collinear loci. Note also that neighboring optimal trajectories have the same structure, except in presence of particular lines, i.e. the switching curve and lines belonging to the sets S and C.

The corresponding optimal synthesis is represented in Fig. 2. In the numerical examples, the switching curve has been determined numerically by considering a series of trajectories with u = +2π originating from the horizontal singular set where Φ = 0. The points of the switching curve correspond to the first point of each trajectory where the switching function vanishes. In Fig. 2, we observe that the optimal trajectory can have different structures extending from bang-bang arcs to bang-singular-bang-bang and bang-singular-bang-singular arcs. A bang-singular-bang-bang control is, e.g., the concatenation of a bang-bang arc followed by a singular and a bang arc. (Color Online) Schematic representation of the optimal synthesis when the initial point of the dynamics is the north pole. An arbitrary zoom has been used to construct the figure. Regular curves are plotted in blue (dark) and red (dark gray) for control fields equal to +2π and -2π, respectively. The optimal singular trajectories are displayed in green (light gray). The dotted line is the switching curve, while the dashed one is the non-admissible part of the horizontal singular line. The small arrows indicate the way the curves are traveled.

The saturation control problem

In this paragraph, we analyze the saturation control problem which consists in bringing the initial state vector to the center of the Bloch ball in minimum time. We compare the optimal control law with an intuitive one used in NMR, the inversion recovery sequence [24,25]. The intuitive solution is composed of a bang pulse to reach the opposite point of the initial state along the zaxis followed by a zero control where we let the dissipation act up to the center of the Bloch ball. The optimal and the intuitive solutions are plotted in Fig. 3. Using results of the previous section, one deduces that the optimal control is the concatenation of a bang pulse, followed successively by a singular control along the horizontal singular line, another bang pulse and a zero singular control along the vertical singular line. The dynamics leaves the horizontal line before the limit of admissibility at a point of coordinates (y 0 , z 0 ). The coordinate y 0 can de determined numerically as follows. We consider the bang arc starting from (y 0 , z 0 ) and we compute Φ and Φ along this trajectory. y 0 is then characterized by the fact that Φ = Φ = 0 at the intersection point of the extremal with the vertical singular line. Figure 3 displays also the switching curve. We have checked that the second bang pulse of the optimal sequence does not cross the switching curve up to the vertical singular axis. In this example, the center of the Bloch ball is reached with an accuracy better than 10 -15 . The durations of the two sequences are 202 ms for the optimal solution and 478 ms for the IR one. One therefore sees that a gain of 58% is obtained for the optimal solution over the intuitive one.

This result also highlights the role of singular extremals. The two solutions considered use singular controls, along the horizontal and the vertical lines for the optimal sequence and along the vertical axis for the IR solution. The difference between the two trajectories comes from the fact that the singular controls are optimal for the optimal law, while it is not the case for the IR one.

The ratio between the two control durations can be analytically computed in the unbounded case where there is no bound on the control. More precisely, this means that we consider the limit ω max → +∞. In this case, the optimal solution is only composed of a bang arc followed by singular controls along the horizontal and the vertical lines since there is no admissibility condition on the control (see also the discussion of Sec. 4.3). We first compute the duration of the IR sequence. In the limit of an infinite bound on the control, the first bang allows to reach instantaneously the south pole of the Bloch sphere. The second part of the control uses the longitudinal dissipation with a zero control field to go to the center of the sphere. The dynamics is governed by:

ż = 1 T 1 (1 -z). (9)
In this computation, note that the parameters Γ and γ and the reduced time τ cannot be used since they depend on the bound ω max . The solution of Eq. ( 9) is

z(t) = (z(0) -1)e -t T 1 + 1 with z(0) = -1.
We next determine the time t 1 such that z(t 1 ) = 0. Inverting this equation, we find the time

t 1 = T 1 ln(1 -z(0)) = T 1 ln 2.
For the optimal solution, the first step is to use a bang pulse to reach instantaneously the horizontal singular line in a point of coordinates (y, z 0 ) (with z 0 = T 2 2(T 2 -T 1 ) and y 2 + z 2 0 = 1). The second step consists of following this line. We know that the dynamics and the singular control satisfy:

   ẏ = -1 T 2 y -u s z 0 u s = -1 2T 1 2T 1 -T 2 T 1 -T 2 1 y . (10)
The dynamics reduces to:

ẏ = - 1 T 2 y - 2T 1 -T 2 4(T 1 -T 2 ) 2 T 2 T 1 1 y , (11)
and if we consider the change of variable Y = y 2 , we get:

Ẏ = -2Y /T 2 -α, (12)
where

α = 2T 1 -T 2 4(T 1 -T 2 ) 2 T 2
T 1 . The solution of this equation is

Y (t) = ( αT 2 2 + Y (0))e -2 t T 2 - αT 2 2 .
Starting from the coordinate Y (0) = 1-z 2 0 and following the singular line up to Y (t) = 0, we find the duration t 2 = T 2 2 ln 1 + 2Y (0) αT 2 . The last part of the optimal solution uses the vertical singular line with u = 0. The same expression as for the I.R. solution leads with z(0) = z 0 to t 3 = T 1 ln(1 -z 0 ). Finally, one arrives at:

   T opt = T 2 2 ln 1 + 2(1-z 2 0 ) αT 2 + T 1 ln(1 -z 0 ) T IR = T 1 ln 2 (13)
where T opt = t 2 + t 3 and T IR = t 1 . As can be expected, the two times only depend on the two relaxation parameters T 1 and T 2 . The optimal solution uses actively the transverse relaxation to reach in minimum time the target state. In our example, the ratio of the two durations is equal to 0.383. A complete numerical study shows that this ratio decreases as ω max increases and tends asymptotically towards the value of the unbounded case [10].

Robustness of the optimal solution

We analyze in this section the robustness of the structure of the optimal solution with respect to variations of dissipative parameters and of the bound on the control. We first consider that the bound of the control field is varied. In this paragraph, we assume that the control field satisfies |u| ≤ 2π × u 0 where u 0 can be changed at will and we fix the dissipative parameters. In this case, three different structures can be found as shown in Fig. 5. When u 0 is sufficiently small, the first bang arc does not intersect the horizontal singular line and the optimal solution is the concatenation of a bang arc followed by a vertical singular arc. When u 0 1, we recover the standard situation described in Sec. 4.2 where the system leaves the horizontal singular line just before the limit of admissibility. It can be shown that this situation is generic in the sense that the intersection point of the initial bang arc with the horizontal singular line belongs to the admissible part of this line. The proof goes as follows. We know that the limit of admissibility is defined by:

y 0 = γ(γ -2Γ) 4πu 0 (Γ -γ) , z 0 = - γ 2(Γ -γ) .
We then determine the limiting case where the bang arc intersects the singular line. This corresponds to the tangency of the two curves, which is schematically represented in Fig. 4. The tangency occurs when ż = 0 and z = z 0 , which leads to the relations γ(1 -z 0 ) + 2πu 0 y = 0 and y = -γ(1 -z 0 ) 2πu 0 = y 0 .

We can conclude that the bang arc is tangent to the singular line when their intersection point is the point of coordinates (y 0 , z 0 ).

In the limit u 0 → +∞, the first bang arc belongs to the Bloch sphere and the limit point of admissibility tends to the vertical axis. This explains the change of structure in the unbounded case where the control is the succession of a bang, a horizontal singular arc and a vertical one.

Figure 6 displays the evolution of the optimal trajectories when the dissipative parameter T 2 varies, T 1 being fixed. Numerical values are taken to be 444, 386, 295, 138 and 60 ms from bottom to top, which correspond in the reduced coordinates to a parameter Γ equal to 0.0697, 0.0803, 0.1048, 0.2241 and 0.5160. If we denote by z = a the position of the horizontal singular extremal then we have T 2 = 2aT 1 /(2a -1). Figure 6 shows that the different extremals present the same qualitative structure with a pulse sequence composed of a bang, a horizontal singular, a bang and a vertical singular extremal to reach the origin. The optimal trajectories only differ in the durations of the bang and singular arcs. Two limiting cases can be considered. The first one corresponds to a = 0, which can be associated in a practical situation to the relation T 2

T 1 or γ Γ, i.e. the longitudinal relaxation is negligible with respect to the transverse one. The optimal solution is then composed of a bang arc followed by a horizontal singular arc with z = 0 up to the center of the Bloch ball. Along this arc, we have u s = 0 since γ = 0. The other limiting case is a = -1, i.e. T 2 = 2T 1 /3, for which the IR sequence is the optimal solution. This conclusion is also true for a < -1.

Experimental implementation

In this work, we consider NMR control experiments in liquid state. Such experiments present some key advantages with respect to other domains of quantum control such as, e.g., the control of atomic and molecular dynamics by laser fields. They can be summarized as follows. The models used in NMR reproduce with great accuracy the experimental conditions even for very complex molecular structures such as proteins in presence of an interaction with an environment. The typical error between theory and experiment is of the order of few percents. In addition, very complex control fields can be shaped experimentally, which reinforces the interest of theoretical studies. To illustrate these properties, we detail below the case of the saturation control problem. Figure 7 presents some experimental data for the time-optimal saturation control problem. The parameters of Sec. 4 have been used in the experiments which were performed on the 7 the very good agreement between theoretical curves and experimental points. This agreement confirms that such kind of optimal trajectories can be implemented with modern NMR spectrometers.

Conclusion and prospective views

Our goal in the research reported in this paper has been to explore the time-optimal control of a spin 1/2 particle in presence of relaxation effects. Completing previous research

Towards the time-optimal control of dissipative spin 1/2 particules in Nuclear Magnetic Resonance15 in this field, this work gives a general overview of the construction and of the determination of the optimal solution, which combines a geometric analysis and analytical and numerical computations. However, it is clear that much remains to be accomplished in the application of geometric optimal techniques to quantum control. Up to now, only very simple quantum systems have been considered. A challenging task in a near future will be to extend such works to more complicated systems. In this perspective, we have applied this method to the optimal control of two uncoupled spin systems in Ref. [11] and we are actually working on the control of three and fours spins. Such studies will be interesting from both a fundamental perspective and for possible applications in NMR spectroscopy and imaging in order to take into account experimental constraints such as resonance offsets and control field inhomogeneities. However, some theoretical and numerical limitations may exist. According to the Pontryagin Maximum Principle, the optimal solutions are computed from a shooting equation, by determining the initial adjoint vector, which selects the extremal trajectory reaching the target state. The resolution of this equation is more and more difficult when the complexity of the system increases or when the chaotic behavior of the associated classical dynamic becomes predominant. In this case, a possibility is to use a continuation method when the control problem depends on a parameter (e.g. the relaxation rates). The continuation approach is used to solve step by step the shooting equation, the solution of an initial problem being known. In addition, even in a very complicated situation where purely numerical algorithms such as the GRAPE one has to be used, the geometric techniques can be interesting. One can find the geometric solution of a simpler problem and then use this control law as a trial solution of the numerical algorithm. This helps the convergence of the algorithm and guides it towards a solution close to the geometric one. We have shown the efficiency of this coupling between the two methods in [23] in order to solve the saturation control problem of a large number of spins with different resonance offsets due to experimental imperfections.

Aknowledgment S.J.G. acknowledges support from the DFG (GI 203/6-1), SFB 631, the EU program Q-ESSENCE, and the Fonds der Chemischen Industrie. Experiments were performed at the Bavarian NMR center at TU Mnchen.

Figure 1 .

 1 Figure 1. (Color Online) Sign of the function g in the plane (y, z). The dashed and dotted lines indicate respectively the positions of the collinear locus and the horizontal singular line (see the text). Two paths γ 1 and γ 2 are plotted in blue (dark) and red (dark gray) to illustrate the use of the clock form.

Figure 2 .

 2 Figure 2. (Color Online) Schematic representation of the optimal synthesis when the initial point of the dynamics is the north pole. An arbitrary zoom has been used to construct the figure. Regular curves are plotted in blue (dark) and red (dark gray) for control fields equal to +2π and -2π, respectively. The optimal singular trajectories are displayed in green (light gray). The dotted line is the switching curve, while the dashed one is the non-admissible part of the horizontal singular line. The small arrows indicate the way the curves are traveled.

Figure 3 .

 3 Figure3. (Color online) Plot of the optimal trajectories (blue or dark gray) and of the Inversion Recovery sequence (green of light gray) in the plane (y, z). The corresponding control laws are represented in the lower panel. In the upper panel, the small insert represents a zoom of the optimal trajectory near the origin. The black line is the switching curve originating from the horizontal singular line. The cross indicates the position of the admissibility point. The vertical line corresponds to the intuitive solution. The blue (or dark gray) curve is the optimal trajectory near the origin.

Figure 4 .

 4 Figure 4. (Color Online) Schematic evolution of a bang arc in green (light gray)starting from the equilibrium point of the dynamics when the maximum bound on the control is decreased (from left to right). The horizontal dashed line is the horizontal singular line, while the red (black) cross represents the limit of admissibility.

Figure 5 .

 5 Figure 5. (Color Online) Plot of the optimal trajectories for different bounds on the control field. The green (light gray) extremal curve does not intersect the horizontal singular line. The blue (dark) trajectory has the generic structure described in Sec. 3, while the red (dark gray) one corresponds to the unbounded case.

Figure 6 .

 6 Figure 6. (Color Online) (top) Evolution of the magnetization vector along the optimal trajectory for different dissipative parameters (bottom) The corresponding control fields u.

Figure 7 .

 7 Figure 7. (Color Online) Plot of the optimal trajectories (blue or dark gray curve) and of the inversion recovery sequence (green or gray curve) in the plane (y, z) for T 1 = 740 ms, T 2 = 60 ms and ω max /(2π) = 32.3 Hz. Filled squares and open diamonds correspond to the experimentally measured points. The lower panel displays the control laws.

  1 H proton spins of H 2 O at room temperature (298 K) in a sample constituted of 10% H 2 O, 45% D 2 O, and 45% deuterated glycerol saturated with CuSO 4 . The experiments were done on a Bruker Avance 600 MHz spectrometer with linearized amplifiers. As already mentioned below, we observe in Fig.