Maximizing entanglement in bosonic Josephson junctions using shortcuts to adiabaticity and optimal control

In this article we consider a bosonic Josephson junction, a model system composed by two coupled nonlinear quantum oscillators which can be implemented in various physical contexts, initially prepared in a product of weakly populated coherent states. We quantify the maximum achievable entanglement between the modes of the junction and then use shortcuts to adiabaticity, a method developed to speed up adiabatic quantum dynamics, as well as numerical optimization, to find time-dependent controls (the nonlinearity and the coupling of the junction) which bring the system to a maximally entangled state.


Introduction
Quantum entanglement [1,2] is the nonclassical correlation between particles whose individual quantum states cannot be described independently even if they are separated by large (even intercontinental) distances [3]; a single quantum state is rather necessary to describe the system as a whole. It is considered to be a unique physical resource, playing a central role in most of the technologies associated with the second quantum revolution [4], including quantum computation and communication [5].
A bosonic Josephson junction (BJJ) is a system of two boson ensembles, with each of them occupying a single quantum state, which interact through a tunnel barrier. Mathematically, this system is described as two interacting nonlinear quantum oscillators. It has been implemented experimentally in various physical settings, with Bose-Einstein condensates (BEC) confined in optical traps [6], atom chip [7], semiconductor microcavities (exciton-polariton systems) [8,9,10], superconducting circuits [11] and photonic systems [12]. It provides an ideal model to study correlations and entanglement in quantum systems.
Shortcuts to adiabaticity (STA) [13,14] is a method to speed up quantum adiabatic dynamics. The idea behind the method is to arrive at the same final state as with a slow adiabatic process, but without necessarily following the instantaneous eigenstates and eigenvalues. Closely related to this is the counterdiabatic-transitionless driving approach [15,16], where an extra Hamiltonian term is added such that the system can be driven along adiabatic paths of the original Hamiltonian. For both cases, the desired transfer can be theoretically completed in arbitrarily short times. In practise, there are always experimental restrictions which limit the STA duration. Optimal control theory (OCT) [17], originally developed during the cold war to answer questions related to the space race, for example the design of minimum-time or minimum-fuel trajectories to the moon, has been well integrated in the STA framework [18] and quantum control in general [19], to evaluate the limits of quantum performance in the presence of realistic constraints. Since the use of adiabatic processes is ubiquitous in quantum dynamics and generally in physics, it is no surprise that STAs have found a wide spectrum of applications. These include the fast cooling and transport of atoms [20,21], BECs [22] and trapped ions [23], the efficient manipulation of two-and three-level quantum systems [24,25], the design of waveguides and photonic lattices [26,27], the optimization of quantum heat engines [28,29,30,31,32,33], suppressing non-adiabatic excitations across a quantum phase transition [34,35], the fast optomechanical cooling [36] and quantum computation [37,38,39], and even the control of mechanical systems [40]. In the context of BJJs, STAs have been exploited for the fast generation of spin-squeezed states [41,42,43] and to expedite the superfluid to Mott-insulator transition [44,45].
In the present work we consider the situation where the two modes of a BJJ are initially loaded with two weakly populated coherent states. Starting from this classical separable state, we use the method of STA to find the time-dependent controls (the nonlinearity and the tunneling rate of the junction) which drive the system to a maximally entangled state. We also express the desired transfer as an optimal control problem and use numerical optimization to obtain the controls which achieve it in minimum time. Note that entanglement generation in BJJs is an active field of research for BECs trapped in optical lattices, with the emphasis given in the semiclassical limit of large occupation numbers [41,42,43,46,47,48,49,50]. Here on the contrary we consider the case of weak pumping, where the occupation numbers remain small. This limit has been investigated in the context of semiconductor microcavities but with constant controls [51], mainly for their potential use as single-photon sources [52,53,54,55].
The current article is structured as follows. In the next section we describe the model of a BJJ initially prepared in a product state of two weakly populated coherent states, while in section 3 we quantify the entanglement between the two modes of the junction. In section 4 we use STAs and numerical optimization to find the timedependent nonlinearity and tunneling rate of the junction which can drive the system to a maximally entangled state. In section 5 we consider the effect of dissipation to the desired transfer, and section 6 concludes this work.

Weakly populated bosonic Josephson junctions
We consider a BJJ described by a quantized two-mode model. The system Hamiltonian in the Bose-Hubbard approximation is [44,45,54,55] whereâ i ,â † i are the creation and annihilation operators at site i, ω is the common resonant frequency of both modes, and U(t), J(t) are the strengths of the nonlinearity and coherent coupling, respectively, which are assumed to be controllable functions of time. Note that, although the coupling rate is a well-known control parameter, the nonlinearity can also be varied in time experimentally [56,57] and it has been exploited in the design of STAs [45,58,59]. As we shall latter explain, the desired transfer cannot be achieved with constant controls.
We assume that the BJJ is initially prepared in a separable product of coherent states where with a small average number of quanta In this weak population limit (4), the system evolution is approximately restricted to the manifold of up to two field quanta and thus the state can be well described by the following truncated wavefunction where |ij is the state with i and j quanta in the two modes, respectively. Note that this approximation has been employed in Refs. [51,54,55]. From Schrödinger equation we find the differential equations governing the evolution of coefficients c ij (t) i d dt Observe that the systems describing the evolution of coefficients in each submanifold of states with the same total number of bosons are independent of each other. This is a characteristic of the evolution under Hamiltonian (1), where each term contains an equal number of creation and annihilation operators. As a consequence, the probability amplitude within each submanifold remains constant. For the first neglected submanifold of states with a total number of three bosons, this probability is of the order of α 6 , as found from the expansion of the initial coherent states, and this is why the truncation (5) is a valid approximation in the limit (4). From the initial state of (2) and the expansion of (3) we find the following initial values for the coefficients The entanglement can be obtained from the relation where λ i , i = 1, 2, 3 are the eigenvalues of ρ 1 . They satisfy the characteristic equation where C = 2 |c 00 c 11 − c 10 c 01 | 2 + |c 10 c 02 | 2 + |c 01 c 20 | 2 + |c 11 c 02 | 2 + |c 11 c 20 | 2 + |c 02 c 20 | 2 is a quantity called concurrence [60,61,62] and In order to find λ i from (11), we first estimate parameters C, D in the limit of (4). From (6)- (8) and the initial conditions of (9a)-(9f) we derive the following constants of the motion Eq. (13c) implies that c 11 , c 20 , c 02 ∼ α 2 , thus D ∼ α 12 → 0 in the limit (4). Additionally, from (13b) we have c 01 , c 10 ∼ α and, if we combine this with the above estimates for the rest of the coefficients, we find from (12) the estimate C ∼ α 2 . Observe now that if D was equal to zero then one of the eigenvalues, let's say λ 3 , would also be zero. Since D is actually a very small perturbation, we can assume that λ 3 remains close to zero and ignore the power terms λ 3 3 , λ 2 3 in (11). Solving the remaining equation for λ 3 we obtain , a very small value, indeed. The other two eigenvalues can be found by solving the characteristic equation with D = 0 and ignoring the zero solution. They are from which we find λ 2 ∼ C 2 ∼ α 4 . The contribution of λ 3 in the entanglement value is actually very small compared to that of the other two eigenvalues. Indeed, thus we can omit the third eigenvalue in (10). We end up with the following expression Observe that the entanglement is an increasing function of the concurrence C, thus it is maximized when C is maximum. From the previously derived estimates of the coefficients c ij , we find that the dominant term in expression (12) for the concurrence is C ≈ 2|c 11 − c 10 c 01 | ≤ 2(|c 11 | + |c 10 c 01 |).
Using the constants of the motion (13b), (13c) we obtain the following bounds where the equalities hold when From (14) and (15a), (15b) we find the maximum concurrence corresponding to the maximum value of entanglement. In the following sections we will derive controls U(t), J(t) which achieve this value.

Maximization of entanglement using shortcuts to adiabaticity and optimal control
In this section we will use the methodology of shortcuts to adiabaticity in order to find controls U(t), J(t) which obtain the maximum concurrence value of (17). Before that we make a couple of useful observations.

Preliminaries
First, note that the diagonal terms proportional to the resonant frequency ω in (7) and (8) simply add a phase factor e −iωt to the coefficients c 10 , c 01 , and another one e −2iωt to c 20 , c 11 , c 02 . This results in an overall phase factor e −2iωt for the complex number 2(c 11 − c 10 c 01 ), corresponding to the concurrence C, which is eliminated by the absolute value operation in (14). Thus we can proceed the analysis as if ω = 0.
The second observation will lead us to the appropriate values of α 1 , α 2 , characterizing the initial coherent states, under the restriction |α 1 | 2 + |α 2 | 2 = α 2 , where α is real and constant. First note the symmetry in system (8) between the first and third variables, in the sense that the system remains invariant if we interchange them. Additionally, for the entanglement to be maximized at the final time t = T , it is necessary that |c 11 (T )| attains its maximum value given in (15b), while c 20 (T ) = c 02 (T ) = 0 from (16b). Thus, the final conditions for the first and third variables of system (8) are also symmetric and, if we propagate this symmetric system backwards from t = T to t = 0 we obtain c 20 (0) = c 02 (0). But from initial conditions (9e), (9f) we find α 2 1 = α 2 2 , thus where, without loss of generality, we have assumed in-phase and real α 1 , α 2 . The antiphase choice α 1 = −α 2 = α/2 leads to a similar shortcut which requires a negative U(t) in order to be implemented, thus we consider only the in-phase case (18). Using this relation and (9b), (9c) we find the following initial conditions for system (7) If we define the vector |ψ = c 10 c 01 , then system (7) can be written in compact form as (recall that we can set ω = 0) where σ x is the Pauli spin matrix [63]. Since the initial state is an eigenstate of σ x , with eigenvalue 1, the above equation can be easily integrated as from which we find Thus, the choice of (18) is not only necessary in order to achieve the maximum |c 11 (T )|, but it also assures that the product |c 10 (t)c 01 (t)| has its maximal value, as determined in (15a). We next concentrate on system (8). Using (18) and (9d)-(9f) we end up with the following initial conditions for this system, which lead to the maximum final value of Starting from (20) and following the procedure described in Ref. [45], one can find controls which drive system (8) to the desired final state of maximum The phase factor in the right hand side of the above equation is acquired during the evolution and will be later chosen such that c 11 (T ) is in phase with the product −c 10 (T )c 01 (T ), whose absolute value is already maximal as explained above, so the concurrence (14) and thus the entanglement are maximized. Following Refs. [44,45], the transfer from (20) to (21) corresponds to the transition from a superfluid state, where each quantum is distributed with equal probability in both modes, to a Mott insulator state, where the two quanta are isolated in separate modes [64]. Here, we will derive exactly the same controls as in [45] without using the Lie algebra of U3S3, but the more familiar Lie algebra of SU(2).

4.2.
Derivation of the shortcut using the Lie algebra of SU (2) It is not hard to verify using (8) that, if we define the vector |ψ as then it obeys the Schrödinger equation with the two-level Hamiltonian where and σ i , i = x, y, z are the Pauli matrices [63]. The initial and final conditions for |ψ are determined from (20) and (21), respectively, and they are Observe that the initial state is an eigenstate of S x , located in the equator of the Bloch sphere, while the desired final state is the spin-down eigenstate of S z , located in the south pole. The phase factor multiplying the final state will be clarified later. At this point, it is instructive to exploit the Bloch sphere picture and explain why constant controls are not suitable for the desired transfer. Observe that, moving the system from the equator to the south pole under Hamiltonian H 0 with constant controls would require U = 2J. In this case, the state is rotated with angular frequency (2U) 2 + (4J) 2 = 2 √ 2J, and this rotation contributes to the phase acquired by c 11 . On the other hand, the product c 10 c 01 acquires a phase determined by the angular frequency 2J, see (19). Since the two frequencies are not commensurate, the two terms c 11 , c 10 c 01 cannot acquire the necessary phase difference which maximizes concurrence, when the controls are restricted to be constant.
In order to find the time-dependent controls which drive system (23) between states (25) along an adiabatic shortcut, we need to diagonalize Hamiltonian (24). If we parametrize U, J as in Refs. [44,45] with time dependent E 0 (t), ϕ(t), then with instantaneous eigenvalues and normalized eigenvectors The time-dependent reference Hamiltonian H 0 (t) can be expressed as with approximate time-dependent adiabatic solutions where the phases are since the inner product term in (32) is zero. According to the transitionless drivingcounterdiabatic approach [15,16], in order to drive the system along the adiabatic path of the reference Hamiltonian H 0 (t), it is necessary to use a modified Hamiltonian where the extra term is given by since the term in the second line of (34) is zero. If the state |ψ satisfies the Schrödinger equation with the counterdiabatic Hamiltonian H(t) then the system evolves exactly along the adiabatic solutions of (31) of the reference Hamiltonian H 0 (t), no matter how short is the duration T .
The construction of the extra term H cd = −φS y through for example a fast switching between S x , S z resulting in their commutator S y , is not a very practical approach [44]. In order to implement the shortcut with a Hamiltonian of the same form as H 0 , we follow an alternative method suggested in Ref. [45]. Consider the wavefunction |ψ I (t) , connected to the state |ψ through the unitary operator B(t) It obeys the alternative dynamics where then |ψ I (t) is an alternative shortcut (we note here that some of the above boundary conditions may be relaxed in certain cases). If we specifically choose where b(t) is a real function of time to be determined, then from (38a), (38b) we find The choice eliminates the undesirable extra term in the second line of (44), and we finally get where Observe that the actual Hamiltonian H I which is used to implement the shortcut has the same form as the reference Hamiltonian H 0 , where the functions U(t), J(t) of the latter have been replaced by the actual controls U I (t), J I (t) in the former.
We summarize the procedure that should be followed in order to obtain correctly the shortcut. We start from system equation (8) with ω = 0 and the actual controls U I , J I instead of the reference functions U, J. Then we define |ψ I as an expression similar to (22) but with U replaced by U I . It obeys Schrödinger equation (37) with H I given in (46). We next move to find the appropriate functions of time ϕ(t), E 0 (t) which determine the reference adiabatic paths. We first discuss the choice presented in Ref. [45] and later provide a new pair of functions which leads to a shorter shortcut for the transfer that we study. Observe from the boundary conditions (25) and the eigenvectors (29) which determine the adiabatic solutions (31) that, for the in-phase (α 1 = +α 2 ) choice of the initial coherent states, the evolution should takes place along |ψ − (t) . The boundary conditions (25) are correctly reproduced when The conditions assures thatḂ(0) = 0. Using polynomials to interpolate the functions E 0 (t), ϕ(t) at intermediate times and imposing on them the above boundary conditions, we find where s = t/T and E max 0 is the maximum value of E 0 (t). Using (53), (54) and (45) at the boundaries s = 0, 1, we have From (43) and (55) we finally obtain Observe that only three out of the four boundary conditions (39)- (42) for B(t) are satisfied, but this does not cause any problem as we explain below.

Calculation of the shortcut duration T
We find the values c 11 (T ), c 10 (T )c 01 (T ) obtained with the shortcut. From (36) we have where as derived from (31), (29), (32), and the final condition (49). Note that expression (58) has exactly the form (25) where ψ I,2 (T ) is the second component of the vector |ψ I (T ) . Using the above equation along with (56), (57) and (58) we obtain where note that we have used (47a), (28) from the first line to the second, and (55), from the second line to the third. Observe that c 11 (T ) has exactly the anticipated form (21), thus the fact that B(T ) = 1 does not affect our analysis. On the other hand, from (19) we have where since we use J I (t) in the original equation (7), not J(t). with k integer. The choice k = −1 provides the shortest feasible phase difference (corresponding also to the minimum time T ) Using (60), (62) and (47b), the above equation becomes and, if we use the substitution s = t/T , we finally obtain where ϕ ′ = dϕ/ds = T dϕ/dt = Tφ. Observe that equation (65), after the integration of its left hand side (LHS) with respect to s, becomes an algebraic equation for the duration T of the shortcut which is necessary to build the desired phase difference. In Fig. 1 we plot the LHS of (65) as a function of T (blue dashed line) and find that the necessary duration to build the desired phase difference π is T = 77.  Fig. 2(e) we display in the complex plane the normalized quantity [c 11 (t) − c 10 (t)c 01 (t)]/α 2 , where note that the surrounding black circle has radius equal to 1 + √ 2, which is the maximum value of the concurrence when normalized with respect to α 2 , see (17). In Fig. 2(f) we plot the time evolution of the normalized concurrence C/α 2 , until it reaches the maximum value 1 + √ 2 (horizontal black line).
We close this subsection by pointing out that one may would like to follow an alternative approach and use other available shortcuts for two-level systems, like for example those in Ref. [25]. The problem in this case is that in the relations corresponding to Eq. (65), which determine the phases for maximum concurrence, the integrals for these shortcuts are independent of the duration T . They only depend on the shape of the shortcut, i.e. its functional form with respect to s, thus it is necessary to introduce extra design variables, something which may complicate the procedure.

A faster shortcut
In the previous subsections we showed that the maximum normalized concurrence 1+ √ 2 belongs to the reachable set of our system, but the necessary duration T to reach this value with the above presented shortcut is quite large. In the case where the undesirable effect of relaxation is present, this long duration may lead to a severe degradation of the performance. For this reason, in the present subsection we derive an alternative, faster shortcut.
We start by finding an estimate of the minimum necessary time to build the π phase difference in (65). The procedure will lead us to some useful observations for the construction of the faster shortcut. First of all we set E 0 (t) = E max 0 , in order to maximize the integral. Next, we ignore for simplicity the second term under the square root in (65), the one which is multiplied by the relatively small quantity 1/T 2 . The LHS of (65) becomes approximately (T E max 0 /2) 1 0 ds(cos ϕ + sin ϕ − 1). Since cos ϕ + sin ϕ = √ 2 sin(ϕ + π/4), obviously the choice ϕ = π/4 maximizes the integrand. Using this optimal constant value for the whole interval 0 ≤ s ≤ 1, in order to find an approximate expression for the LHS of (65), and then solving for T we obtain In order to achieve a duration T close to the above estimate, we construct a shortcut that mimics the desirable characteristic identified above, i.e. E 0 (t) stays close to E max 0 and ϕ(t) close to π/4. The following function remains equal to E max 0 until s = s 0 , where s 0 is a design parameter. For s ≥ s 0 we choose a polynomial form to satisfy the boundary conditions at s = 1(t = T ) as well as the smoothness conditions at the junction point Note that at the junction point we require the continuity ofË 0 such that the control U I , which depends onĖ 0 , see (47a), has a continuous derivative there. The coefficients for a specified value of the design parameter s 0 . For the function ϕ(s) we choose the following form thus it remains equal to the constant value π/4 in the interval s 1 ≤ s < s 2 , where s 1 , s 2 are design parameters. For s < s 1 and s ≥ s 2 we choose polynomial forms to satisfy the boundary conditions at the initial and final points, as well as the smoothness conditions at the junction points Note again that at the junction points we require the continuity of the third order derivative so the control U I , which depends onφ, has a continuous derivative there. From (72) we obtain while the rest coefficients a j , j = 3 . . . 6 are chosen to satisfy (74) and are found by solving numerically the linear system for a specific value of s 2 . For a concrete example, we pick s 0 = 9/10, s 1 = 2/10 and s 2 = 1 − s 1 = 8/10. The durations of the transient intervals are chosen short enough but not too short, to avoid negative values in U I . In Fig. 1 we plot the LHS of (65) as a function of T (red solid line) and find that the necessary duration to build the desired phase difference π with this shortcut is only T = 15.665 units of time (E max 0 ) −1 . This value is much shorter than the previous one and close to the estimate of (66). In Figs. 2(a) and 2(b) we plot ϕ(s) and E 0 (s) from (71) and (67)

Maximization of concurrence using optimal control
Having shown that a maximally entangled state is reachable and in order to evaluate how well performs the faster shortcut introduced above, we apply an optimal control approach and find the minimum necessary time to reach the maximum value of the normalized concurrence with bounded controls where the upper bounds are chosen close to the maximum values of the shortcut controls, see Fig. 2. This complementary procedure, see for example our work [65] on the expansion of Bose-Einstein condensates and Refs. [66,67] in the context of quantum statistical mechanics, is important since the corresponding controls may be useful under different experimental constraints, while note that optimal control has been exploited for entanglement maximization between two qubits [68,69]. We use the freely available optimal control solver BOCOP [70] to numerically solve a series of optimal control problems with increasing duration T and objective the maximization of the final normalized concurrence C(T )/α 2 . Note that in the BOCOP software package, the continuous-time optimal control problem is approximated by a finite-dimensional optimization problem, using time discretization. The resultant nonlinear programming problem is subsequently solved using the nonlinear solver Ipopt. For the current problem we use a time discretization of 1000 points. With the controls restricted as in (79) we find that the minimum necessary time to achieve the maximum value 1 + √ 2 is T = 6.71 units of time. In Fig. 2 we plot (magenta dashed-dotted line) the corresponding controls, as well as [c 11 (t) − c 10 (t)c 01 (t)]/α 2 and C(t)/α 2 . Observe that the shorter time obtained with the optimal control approach is achieved with non-smooth controls (a typical behavior for minimum-time problems), in contrast to the smooth controls corresponding to the adiabatic shortcut, thus optimal controls might be more difficult to implement experimentally. Another characteristic of the minimum-time controls is that, the larger is the maximum allowed amplitude, the shorter is the necessary time to reach the target.

The effect of dissipation
We can incorporate dissipation in our system's evolution using the following master equation for the density matrix ∂ρ ∂t are Lindblad terms expressing losses to the environment at rate κ. This equation is actually derived from a stochastic Schrödinger equation which includes random quantum jumps. These random jumps become rare for vanishing occupation numbers of the modes, which is the case in the weak pumping limit and under the presence of dissipation. As a consequence, the non-diagonal Lindblad terms 2â j ρâ † j , j = 1, 2 can be neglected and the density matrix equation becomes [51,54,55,71] ∂ρ ∂t where the effective non-Hermitian Hamiltonian is Under this evolution the density matrix can be factorized as where state |ψ(t) satisfies the Shrödinger equation The evolution described by the above equation can be correctly accounted for with the simple substitution ω → ω − iκ/2 in system equations (7), (8). Consequently, the effect of dissipation is to multiply the dissipationless values of c 10 , c 01 and c 11 with e −κt/2 and e −κt , respectively, while the concurrence (14) is reduced by a factor of e −κt .
In Fig. 3(a) we plot the effect of dissipation in the evolution of normalized concurrence for the faster shortcut and for various values of the dissipation rate (from top to bottom κ = 0, 0.01, 0.05, 0.1, units of E max 0 ). The top curve (κ = 0) is actually the same with the red solid line shown in Fig. 2(f), while the rest of the plots are obtained by multiplying this curve with the corresponding dissipation factor e −κt . Observe that there is an overall degradation of the performance, while the maximum of each curve is shifted towards earlier times, since the dissipation factor is a decreasing function of time, and this shift is larger for larger dissipation rates.
In Fig. 3(b) we plot the maximum normalized concurrence which can be obtained under constraints (79) with an optimal process of duration T , from T = 1 to T = 7 with an increment ∆T = 0.1, for various values of the dissipation rate (from top to bottom κ = 0, 0.01, 0.05, 0.1, units of E max 0 ). For the top curve, which corresponds to the absence of dissipation (κ = 0), the performance is actually a non-decreasing function of the duration T . This can be easily explained since, for a larger duration T ′ > T , the same performance can be obtained in the interval [0 T ] and then set the controls to zero and do nothing in the remaining interval (T T ′ ]. Observe also that there is a discontinuity in the slope of this curve, which is due to the fact that the optimal pulse sequences change shape at this point from less to more switchings, as shown in Figs. 2(c), 2(d). This is a kind of behavior that we have encountered several times in our previous work on optimal control of quantum systems, see, for example, Ref. [72]. The three lower curves are obtained from the upper curve by multiplying it with the corresponding dissipation factor e −κt . The overall performance is decreased as before but now the degradation is milder, since here we deal with the optimal processes. The duration T corresponding to the maximum normalized concurrence is shifted towards earlier times for larger dissipation rates since, in the presence of dissipation, the waiting with zero controls is not free but comes with an exponential cost.

Conclusion
In this article, we used the methods of shortcuts to adiabaticity and optimal control to obtain time-dependent controls which can drive a bosonic Josephson junction, initially prepared in a product of weakly populated coherent states, to a state of maximum entanglement between the two junction modes. As controllable variables, we considered the nonlinearity and the tunneling rate of the junction. The present work may find application in the variety of physical contexts where a bosonic Josephson junction can be implemented.  Figure 2. In all the above figures, blue dashed lines correspond to the shortcut described in [45] and in Ref. 4.2, red solid lines to the shortcut introduced in subsection 4.4, and magenta dashed-dotted lines to the time-optimal process obtained in subsection 4.5 (a,b) Functions ϕ(s), E 0 (s) for the two shortcuts, s = t/T and T = 77.724, T = 15.665 units of time, respectively (b,c) Controls U I , J I for the two shortcuts and the optimal process as functions of s = t/T , where T as before for the shortcuts and T = 6.71 for the optimal case (e) Evolution of the quantity [c 11 (t) − c 10 (t)c 01 (t))]/α 2 on the complex plane for the three cases. The black circle of radius 1+ √ 2 corresponds to the maximum normalized concurrence (f) Evolution of the normalized concurrence C(t)/α 2 for the three cases. The black solid line corresponds to the maximum value 1 + √ 2, which is obtained for the three cases at the different durations mentioned above.