Optimal Pulse Design for Dissipative-Stimulated Raman Exact Passage

Quantum control of lossy systems is known to be achieved by adiabatic passage via an approximate dark state relatively immune to loss, such as the emblematic example of stimulated Raman adiabatic passage (STIRAP) featuring a lossy excited state. By systematic optimal control study, via the Pontryagin maximum principle, we design alternative more efficient routes that, for a given admissible loss, feature an optimal transfer with respect to the cost defined as (i) the pulse energy (energy minimization) or (ii) the pulse duration (time minimization). The optimal controls feature remarkably simple sequences in the respective cases: (i) operating far from a dark state, of π-pulse type in the limit of low admissible loss, or (ii) close to the dark state with a counterintuitive pulse configuration sandwiched by sharp intuitive sequences, referred to as the intuitive/counterintuitive/intuitive (ICI) sequence. In the case of time optimization, the resulting stimulated Raman exact passage (STIREP) outperforms STIRAP in term of speed, accuracy, and robustness for low admissible loss.


Introduction
Quantum control often faces losses or noise that have to be circumvented in order to make operational the modern quantum technologies, such as in quantum information processing [1,2]. Two intuitive opposite strategies have been generally developed: (i) operating sufficiently fast [3], i.e., in a timescale in which the loss does not act, or (ii) operating adiabatically along a dark state, which is by construction immune to loss, as in the emblematic example of stimulated Raman adiabatic passage (STIRAP) [4][5][6][7]. The first strategy is an intensive area of research, including optimal control [2,8], shortcut to adiabaticity [3], single-shot shaped pulse [9][10][11], and robust optimization [12][13][14][15]. Optimal control theory (OCT) is a powerful tool to mitigate intensities of the control pulses and speed up the evolution allowing one in principle to reach the ultimate bounds in the system (often referred to as the quantum speed limit when optimal time is considered) [16]. Besides numerical implementation of OCT, such as gradient ascent algorithms (GRAPE) [17], the Pontryagin maximum principle (PMP) [8,[18][19][20] allows analytic or semi-analytic derivation of the optimal controls (typically with respect to time or energy) in transforming the initial infinite-dimension control problem into a finite-dimension problem. In the second strategy, the process is, however, approximative due to adiabatic principles, but it is also affected by losses, since the dark state is only an approximation of the dynamics for a realistic finite time, i.e., not strictly adiabatic, process.
In this paper, we determine optimal routes in a Λ-system, where the upper excited state, coupled to both ground states, is lossy, yielding necessarily to a lossy dynamics when the two ground states are not directly coupled. We follow an inverse-engineering strategy: we start by assuming a total loss (at the end of the dynamics) that is defined a priori as an admissible loss. We then determine using the PMP the optimal dynamics with respect to a given cost (the pulse duration or the energy) that will reach the target state and realize the admissible loss. As a consequence, the transfer to the target state is incomplete due to the loss, but the fidelity of the transfer is a priori known and controlled. We then refer this process to an optimal dissipative stimulated Raman exact passage (optimal dissipative STIREP), similarly to its definition in the non-linear case [21].
When the pulse energy is considered as cost (for a given pulse duration), the optimal dynamics operates with a relatively strong field far from a dark state. In the limit of low admissible loss, one can reformulate the problem as the control of a planar pendulum and an analytic expression of the coinciding pulses, as a sine Jacobi elliptic function, is provided.
On the other hand, when the pulse duration (with constrained pulse amplitudes) is taken as cost to minimize, we explicitely determine the optimal dynamics. A key result of this work is established in the limit of low admissible loss: we derive a faster, more accurate and more robust dynamics than the traditional STIRAP. The resulting optimal pulse sequence is remarkably simple, featuring a slow counterintuitive order, reminiscent of the STIRAP sequence, but sandwiched by two sharp intuitive sequences.
The paper is organized as follows: in Section 2, we define the system. In Sections 3 and 4, we determine the energy-optimal and time-optimal dissipative STIREP processes, respectively. We conclude in Section 5, where a comparative analysis and in particular the robustness of the derived optimal processes are provided. Appendices with the details of the calculations complete the paper.

Definition of the Lossy-Driven Raman System
We consider a three-level Λ-system, in which the Hamiltonian in rotating wave approximation (in units such thath = 1) is given in the basis of the states {|1 , |2 , |3 } by: where u p and u s are the pump and Stokes controls, respectively (corresponding to half of the traditional Rabi frequencies), with the lossy upper state |2 , via the dissipation rate Γ. Instead of analyzing such a complicated lossy system requiring specific adaptation of the cost [22], large dissipation [23], or assumption on the controls [24], we consider an alternative procedure to treat approximately but accurately the problem in the situation of interest having a relatively low dissipation rate without adapting the cost, nor restricting the controls. From the unlossy system H ≡ H Γ=0 , the effects of the loss are taken into account at the second order perturbation theory from the knowledge of the state amplitude of state |2 (in absence of dissipation): where t i (t f ) is the initial (final) time, the state amplitude of state |j is denoted c Γ,j in presence of the dissipation, and c j when Γ = 0, respectively. This simplification can be shown to be numerically accurate for a low enough ratio Γ/u max 0.1, where u max is the peak value of the pulses. We consider that the pulses are exactly resonant because they produce the most efficient two-photon coupling [25]. We denote the state solution |ψ = [c 1 , c 2 , c 3 ] T . We use the alternative notation x 1 ≡ c 1 , x 2 ≡ ic 2 , and x 3 ≡ −c 3 so that all the elements of the solution |x ≡ [x 1 (t), x 2 (t), x 3 (t)] T are real when we consider the initial (real) state c 1 (t i ) = 1, c 2 (t i ) = c 3 (t i ) = 0, and the time-dependent Schrödinger Equation (TDSE) becomes: Throughout the paper, we consider an admissible loss, which is, according to (2), characterized by: where A, the time area of the population in state |2 , is a given constant, referred to as the (total) normalized loss (with respect to Γ). The dynamics is determined with the non-lossy Hamiltonian H Γ=0 and the expected loss is, thus, given by P loss ≈ AΓ if one assumes Γ/u max 0.1.

Construction of the Pseudo-Hamiltonian and Derivation of the Equations of Motion from PMP
The goal of the control is to steer the system from x 1 (t i ) = 1 to x 3 (t f ) = 1 (chosen as the target; similarly, we could have chosen x 3 (t f ) = −1) in a fixed time T = t f − t i while minimizing the energy of the controls: under the constraint (4).
To take into account constraint (4), we augment the dimension of the system with a new coordinate y(t) such that: of initial y(t i ) = 0 and final value: The constraint (4) reduces, thus, to a boundary problem on y.
Similarly to the unconstraint optimization case [8], it is convenient to use angle coordinates, which simplify the representation of the dynamics from three components to two angles: with the initial conditions ϕ(t i ) = 0, θ(t i ) = 0. The equations of the dynamics (3), complemented by (6), can be simplified as: after a rotation on the control fields: leading to an invariant cost on the new field variables since u 2 p + u 2 s = v 2 p + v 2 s . According to the PMP [18,26,27], the minimization of the energy (5): leads to the control Pontryagin Hamiltonian (see Appendix A, where we have considered the standard choice p 0 = 1/2): where Λ = [λ ϕ , λ θ , µ] T is the co-state gathering the conjugate momenta of ϕ, θ, and y, respectively. The Hamilton equations lead to the equations of motion (9) and (10) and to: This implies that λ θ and µ are constants of motion. The maximization condition of the PMP gives: which yields: leading to: which features an effective autonomous system (i.e., explicitly time-independent, since H c depends only on the dynamical variables and their conjugate momenta). The equations of motion finally read:φ with the boundary conditions (for a complete population transfer from state |1 to state |3 ):
As detailed in Appendix A, we first integrate Equation (26) and obtain sin ϕ (A24) as a function of an elliptic sine of time 0 ≤ t ≤ T/2: where the set of the three parameters {µ, λ θ , h} is replaced by the set {m, ϕ 0 , T} with the following correspondence: These are determined in order to satisfy the boundary conditions for a given normalized loss A.
In order to take into account the boundary condition y(t f ) = A of (21), sin ϕ has to satisfy: From the integral of the Jacobi elliptic function in (28), we obtain: where E(m) is the complete elliptic integral of the second kind. This gives for the optimal unconstrained system m = 0 [8]: We next integrate Equation (27), as detailed in Appendix B. Imposing by symmetry that θ(ϕ 0 ) = π/4 leads to the condition: which gives an implicit relation between sin ϕ 0 and m. For a given value of loss A [below 3/8 (32), the latter corresponding to the optimal unconstrained system], we aim at finding the optimal couple of parameters m and ϕ 0 , simultaneously satisfying Equations (31) and (33). The obtained data are shown in Table 1. A 3.12 3.14 3.14 3.14 3.14 One can conclude that for a decreasing admissible loss A, m decreases and ϕ 0 decreases to π/4, both monotonically since the right-hand side goes to ϕ 0 when m → −∞. The minimum peak value of ϕ(t) is, thus, π/4, asymptotically reached in the limit of no admissible loss, A → 0.

Derivation of the Pulses and the Dynamics
The original controls u p and v p are obtained by reversing Equation (12), where θ is given by Equations (A31) and (A32), and v p , v s by Equation (17) with ϕ obtained from Equation (28) and the correspondence (29).
From the values m and ϕ 0 shown in Table 1, one can derive the controls and determine numerically the corresponding population dynamics. Figure 1 shows such dynamics for three typical values of A.
We observe that the pump and Stokes controls operate in the so-called intuitive order (first pump and next Stokes), and that they get closer and coincide more and more with a larger peak for decreasing A. One can notice in Figure 1 that the pump and Stokes pulses are almost already undistinguishable at the scale of the figure already for A = 0.08T.
We determine that, in the limit of a small normalized loss A, i.e., A 0.08T, corresponding to a large negative m and ϕ 0 = π/4, the pulses appear as fully overlapping of sine Jacobi elliptic function form: with m solution of Equation (31): which is well approximated by: In this limit, the peak of the pulses is well approximated by: the area of each pulse is: and the generalized Rabi frequency pulse area: which has to be compared to the optimal Rabi frequency pulse area, which is √ 3π (the counterpart of the π-pulse for Λ systems), obtained in absence of constraint on loss [8].
We conclude that the present energy-optimal pulse operates far from a dark state, of πpulse type, in the sense that it is relatively close to the optimal Rabi frequency pulse area.

Comparison with Standard STIRAP and Parallel STIRAP
The energy of the derived optimal STIREP pulses with a low normalized loss A 0.036T (featuring almost overlapping pulses) is compared to that of the standard STIRAP with Gaussian pulses in a situation giving the same normalized loss A. We obtain the energy E 90¯h T for STIRAP, almost twice larger than energy-optimal STIREP E 56¯h T . The optimal pulses feature a shorter duration, compared to the long adiabatic process of STIRAP, and are more intense with the peak u max 19.6/T three times larger than those of the standard STIRAP u max 6/T. We next compare the energy-optimal dissipative STIREP with respect to Stimulated Raman parallel adiabatic passage (parallel STIRAP) with coincident pulses [28], where non-adiabatic transfer is ideally cancelled (in the adiabatic limit) [29], see Appendix D. Numerical implementations are shown in Figure 2 for Ω 0 = 14.453/T, in a situation giving the same energy of the pulses as that for energy-optimal STIREP E 56¯h T . We obtain the time area of the transient population in the excited state A 0.535T, much larger than for energy-optimal STIREP A 0.036T. Comparing the pulses, we first notice that the peak in parallel STIRAP u max 4.2/T is roughly five times smaller than the one in energy-optimal STIREP u max 19.6/T. The duration of the control process is greater than for energy-optimal STIREP.

Analogy with a Planar Pendulum
We can reformulate the optimization problem in the limit of fully overlapping equal control fields: and show that it becomes analogous to the planar pendulum tipping over from the initial unstable vertical position. We denote by a(t) the area at time t of u(t): In this case, the dynamics can be exactly integrated: showing that the target state is reached if: Minimizing the energy of the controls (5) under the constraint of a given loss (7) is equivalent to finding the optimal solutions which minimize the energy T 0 u(t) 2 dt = T 0ȧ 2 dt while satisfying the constraint condition of an admissible given loss For that purpose, we introduce the Lagrangian equation: where λ is a Lagrange multiplier, aiming at minimizing T 0 L λ (t)dt. A necessary condition is given by the Euler-Lagrange principle: giving the equation of a planar pendulum: where λ is directly connected to the frequency of the pendulum, a (up to a factor) is the angle (where 0 corresponds to the unstable vertical positon), and u(t) is the angular velocity. This equation can be integrated by using Jacobi functions with the boundary conditions a(0) = 0 and a(T) = π/ √ 2: with m a constant related to the initial velocity u(0) and λ: We derive, using b = √ 2a: and π/2 giving the relation between λ, m and T, i.e., by inverting (49): and One can notice that the initial velocity u(0) = √ 2K(m)/T is strictly positive, which induces the pendulum to tip over. In the limit of large negative m, this initial velocity goes to zero and we recover the pulse shape (34):

Construction of the Pseudo-Hamiltonian
In the time-optimal control problem, we minimize the following equation: where T is the control duration of the process to be determined optimally, under the constraint on the total peak amplitude of the fields, bounded by a constant u 0 : The pesudo-Hamiltonian, which includes the subtraction of the integral kernel of the above cost, i.e., the constant −p 0 , to which we add p 0 , reads: where the costate Λ has four components The adjoint equations of the costate are as follows: The pseudo-Hamiltonian is of the form (A45) H c = H 0 + u p H p + u s H s (see Appendix E), with the control variables u p and u s , and H p = λ 2 x 1 − λ 1 x 2 , H s = λ 3 x 2 − λ 2 x 3 , we can, thus, apply the results of Appendix E: This leads to: and the controls attain the maximum of the constraint at each time: We can make a change of variables for the time, renormalizing the field amplitude as follows:ũ i.e.,ũ 2 p +ũ 2 and the equation becomes: This means that we can always renormalize the field amplitudes by modifying the optimal time accordingly. We will consider below the tilde variables (corresponding to finding the optimal time for u 0 = 1).
Introducing the following angle coordinates: The equations of the dynamics (69) can be simplified as follows: after a rotation on the control fields: leading to an invariant cost on the new field variables sinceũ 2 p +ũ 2 s =ṽ 2 p +ṽ 2 s = 1. We arrive at: with Λ = [λ ϕ , λ θ , µ] T being the costate gathering the conjugate momenta of ϕ, θ, and y, respectively. The following Hamilton Equations: lead to the equations of motion (71) and to: This implies that λ θ and µ are constants of motion. H c (73) is again of the form (A45), implying: with The equations of motion read then: with the boundary conditions:

Construction of the Optimal Trajectories from the PMP
Since H c = h is a constant, we obtain: Following the same lines as in the energy minimum case, we have: with the normalized constants:μ We know that ϕ has by symmetry to reach its maximum value at t = T/2:φ(T/2) = 0, i.e., (1 −μ sin 2 ϕ 0 ) 2 =λ 2 θ tan 2 ϕ 0 (91) with the notation ϕ 0 ≡ ϕ(T/2). This equation shows a dependence uponλ 2 θ . We can, thus, limit our study to positiveλ θ . It is solved in Appendix F.
We next solve the differential Equation (89) numerically and determine which values ofμ andλ θ satisfy the following: where the left-hand side is the numerical solution of (89) and the right-hand side the possible solution(s) ϕ 0 , as determined in Appendix F. We note that forμ ≤ −8, (92) is satisfied only for the smallest root ϕ 0 [i.e., (A77) for k = 1], since ϕ (assumed positive) grows from 0 to ϕ 0 , where dϕ/dθ = 0. In addition, we emphasize that the results become strongly sensitive toλ θ for large negativeμ, which necessitates a high precision on the estimation of the parameter. Forμ ∈] − 8, 1], there is one root of ϕ 0 (A74), and (92) can be always satisfied. Forμ > 1, no solution satisfying (92) exists. The couplesλ θ ,μ and the corresponding ϕ 0 are shown in Table 2.
Using these values of the couple (μ,λ θ ), we can solve the differential Equation (88), to determine the value of t 0 when ϕ(t 0 ) = ϕ 0 , and thus to obtain the optimal time T = 2t 0 , and the resulting normalized loss: We show the corresponding data in Table 2, from which one can conclude that the optimal control time T decreases for increasing values of A with diminishing rates and gradually flattens out for large A. This is shown in Figure 3. We observe from Table 2 that, for the unconstrained case, i.e.,μ = 0, the determined time area is A 1.0203/u 0 , while the optimal time is T = 2.7207/u 0 , which recovers the unconstrained result obtained in the energy-optimal dissipative STIREP (see Table 1) [i.e., Equation (32)]: A = 3T 8 1.0203/u 0 , and ϕ 0 = π/3 ≈ 1.0472. Table 2. Time-optimal ϕ 0 as a function of the coupleλ θ ,μ satisfying (92), and the corresponding normalized loss A, minimal time T, ratio A/T, pulse area A, and energy E . The valuesμ are given exactly and all the others quantities are approximate.

Derivation of the Pulses and the Dynamics
From a given pair (μ,λ θ ), ϕ(t) can be obtained numerically by solving the differential Equation (88), and then λ ϕ (t) is derived from (87). The original controlsũ p ,ũ s are obtained from Equation (72), where the angle θ is derived numerically from Equation (89): Figures 4-6 show the parameters ϕ(t) and θ(t) and the dynamics for three typical couples of (μ,λ θ ) with decreasing losses. Starting with the full intuitive dynamics in the unconstrained case (Figure 4), we observe that, as the optimal time increases (i.e.,μ becomes larger in absolute value and negative), the pump pulse decreases sharply at early times before slowly increasing and the Stokes pulse increases sharply at early times before slowly decreasing (and a symmetric situation at late times). This corresponds to a slow counterintuitive sequence, reminiscent of the STIRAP sequence, sandwiched by two fast intuitive sequences. This remarkable simple optimal sequence, referred to as an intuitive/counterintuitive/intuitive (ICI) sequence, represents an important finding of our paper. Figures 4-6 show that this behavior is more pronounced, i.e., with a sharper intuitive sequence with higher peak amplitude, for a decreasing admissible loss.
which we define with the actual θ, gets closer to one for a smaller admissible loss. We can notice that this projection is in fact 1 minus the population in state |2 .

Conclusions and Discussion
In this paper, we have derived energy-and time-minimum optimizations under the constraint of a given admissible loss leading to exact state-to-state transfer of three-level Λ-type quantum systems. The resulting control fields have very different shapes depending on the considered optimization.
In the case of energy-optimal dissipative STIREP, the pulses feature an intuitive sequence similar to the unconstrained situation, of amplitudes increasing with the decrease of the given loss, tending to coincident pulses below the loss A 0.08T. We obtain in this case the analytic expressions of the pulses, as a sine Jacobi elliptic function. The energyoptimal dynamics operates with a relatively strong field in the limit of a low loss, and far from a dark state.
In the case of time-optimal dissipative STIREP, the control time increases with the decrease of the given loss A. The process features a remarkable pulse sequence: a relatively slow counterintuitive sequence sandwiched by sharp intuitive sequences (see Figure 6) and the ICI sequence, sharing, thus, some similarities with STIRAP except the very beginning and the very end of the process, which are strongly non-adiabatic (see below for a detailed comparison). The time-optimal strategy operates relatively close to the dark state in the limit of a low admissible loss.
The process duration, peak pulse amplitude, energy, and area for several given admissible losses are reported in Tables 1 and 2. One can compare a few values in order to highlight the main features and differences of the two optimization strategies:

•
The energy optimization with the (low) given admissible loss A = 0.05T with T ≡ T EO the time of the process, referred to as EO, yields A ≈ 0.7/u max,EO with u max,EO the peak of the pulse, T EO ≈ 14/u max,EO , the energy E EO = 40h/T EO ≈ 2.8u max,EO , and the pulse area A EO = π; • The time-optimization with the same loss A = 0.05T and the same peak amplitude u 0,TO1 ≡ u max,EO , referred to as TO1, is roughly obtained forμ = −5: A ≈ 0.7/u 0,TO1 , T TO1 ≈ 3.4/u 0,TO1 , E TO1 ≈ 11.5h/T TO1 ≈ 3.4u 0 and A TO1 = 3.4; it shows a much smaller (roughly four times smaller) time of processing, but slightly larger pulse area and energy; • The time-optimization with the same loss A = 0.05T and the same duration as the energy optimization: T TO2 ≡ T ≡ T EO , referred to as TO2, is roughly obtained for µ = −20.5: this leads to a significant larger energy E TO2 = 49h/T and a twice larger pulse area A TO2 = 7, but to a (twice) smaller peak pulse amplitude u 0,TO2 ≈ 7/T. We show the dependence of the pulse amplitude on the duration corresponding to these three examples in Figure 7. The energy-minimization strategy can, thus, achieve the transfer for a given loss in a relatively small pulse area, but with a relatively large pulse peak amplitude due to its sharp shape. On the other hand, the time-minimization strategy can achieve it with a weaker pulse amplitude, but for a larger pulse area (and energy).
In Figure 8, we study the robustness as a function of a relative deviation α of the amplitudes, i.e., with the amplitudes u j (1 + α), j = p, s, taking into account explicitly the loss, i.e., using the Hamiltonian Equation (1). We consider the three cases: energy optimization with A = 0.05T (EO), time optimization (TO2), and the associated STIRAP. The latter is defined as the traditional counterintuitive configuration of the pumps and Stokes pulses with sine and cosine shapes, respectively, that fit well the actual TO2 pulses except the initial and final sharp intuitive sequences (see the upper frame of Figure 8). We observe that the time-optimal ICI sequence TO2 features a flat asymmetric profile on a relatively large zone, and that it is much more robust than its associated STIRAP. The lack of robustness of the latter is expected, as the total pulse area is too low to achieve efficient adiabaticity.  Comparing with the energy-optimal dissipative STIREP, the time-optimal ICI sequence is much more robust. The weak robustness of the energy-minimization pulse sequence can be understood by the fact that this process can be interpreted as the counterpart of the two-state π-pulse transfer to the Λ system with a large amplitude and a short time in order to respect the admissible loss.
We conclude that the time-optimal ICI sequence achieves precise and fast population transfer, with a chosen loss, and in a robust way (in the limit of low loss).
We notice that the time-optimal ICI sequence is very similar to the pulse sequence derived in [15] (compare Figure 6 with Figure 3 of Ref. [15]),where the optimization was determined with the explicit constraint of robustness, but without consideration of loss.
The implementation of the ICI sequence with a high fidelity of error P loss 10 −3 (P loss ≈ 10 −3 is the situation considered in Figure 8) corresponds to T 2 × 10 −2 /Γ with the typical value A = 0.05T and u 0 = 7/T, i.e., u 0 350Γ.
For a practical implementation, we have to consider the additional decay within the Λ system, including decay channels from state |e to the two ground states |1 and |3 and the associated decoherence, which has to be analyzed with the density matrix formulation and the Lindblad equation, as detailed in Appendix G. Numerical analysis has been conducted for the time-optimal situation, which is the situation of interest since it features robustness and low loss for large enoughμ. We obtain that the two additional channels |e → |1 and |e → |3 (associated with the respective rates γ 1 and γ 3 ) add each a contribution of half the external loss given by the Γ-channel (for the same rate). We obtain more specifically for the total loss: P loss ≈ P (Γ) with the external loss given by the Γ-channel defined by (2): P In practice, we have this to consider the total loss P loss (95).
Concerning a possible experimental implementation, one can mention the excited state 1 D 2 of a praseodymium ion in a Pr 3+ :Y 2 SiO 5 crystal with Γ −1 ≈ 164 µs [30,31], which requires the Rabi frequency u 0 2π × 340 kHz and the duration T 3.3 µs, well in the achievable experimental range.
Author Contributions: X.C. and S.G. proposed the initial idea about the investigation of optimal dissipative STIRAP by PMP. S.G. proposed its determination using perturbation theory to relate the loss to the area of the population in the upper state of a lossless dynamics. D.S. supervised the implementation of the PMP. K.L. and S.G. carried out the numerical simulations. All authors participated in the analysis of the results and in the paper redaction. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Pontryagin's Maximum Principle
To determine optimal control fields u(t) of a dynamical systemẋ = f x(t); u(t) (of dimension N) with respect to the minimization of a given cost: Pontryagin's maximum principle transforms the initial infinite-dimension control problem into a finite-dimension problem and allows discontinuous controls [18,20]. It provides necessary conditions for optimality. It states that the trajectories of the extremal vector x(t) and of the corresponding adjoint state p(t) formed by the Lagrange multipliers, p(t) ≡ [p 1 (t), · · · , p N (t)], fulfill Hamilton's equations: associated with the control pseudo-Hamiltonian equation: where the constant p 0 > 0 can be chosen for convenience since it amounts to multiply the cost function by a constant. For almost all t i ≤ t ≤ t f , the function H c p(t), x(t); u(t) is maximum at certain controls v(t) = u(t), for which one can write H c p(t), x(t); v(t) = const., i.e., The costate λ defined via the conjugate moments: λ = p T , is solution of the second Hamilton equation: Appendix B. Integration of Equation (26) In this Appendix, we solve the differential Equation (26): being one of the equations of motion derived from PMP in the energy-optimal case (see Section 3). From the fact that ϕ(T/2) ≡ ϕ 0 is maximum at t = T/2, we plug the identitẏ ϕ(t = T/2) = 0 into (26), which gives the relation: with the parameters: as a function of which one can express two possible values of sin 2 ϕ 0 : In general, this solution exists if β 2 ≥ 4α and (β ± β 2 − 4α)/α > 0. If α > 0 and β 2 > 4α, there are two roots ϕ 0,± . If α < 0 and β < 0, only ϕ 0,+ exists. If α < 0 and β > 0, only ϕ 0,− exists.

Appendix C. Integration of Equation
being the second equation of motion derived from PMP in the energy-optimal case (see Section 3). The sign ± is the same as the one ofφ, which, taking into account that θ = 0 when ϕ = 0, leads to: From the result: where Π(n; Φ|m) is an incomplete elliptic integral of the third kind: we identify: and Equation (A26) becomes for the positive branch (for which ϕ increases from 0 to ϕ 0 ): Imposing the symmetrical negative branch: as in the unconstrained case µ = 0, leads to: which gives an implicit relation between sin ϕ 0 and m. We can recover for µ = 0 (i.e., α = 0 and m = 0) ϕ 0 = π/3 (solution of (33)), giving , i.e., a = 0 and b = 4/3, and the positive branch: Using: we obtain for µ = 0:
Contour plots as a function ofμ andλ θ in Figure A1 show the solution(s) of ϕ 0 satisfying (91). Considering that the maximum number of solutions is three, we make three plots.