Quantum brachistochrone for multiple qubits

Efficient control of qubits plays a key role in quantum information processing. In the current work, an alternative set of differential equations are derived for an optimal quantum control of single or multiple qubits with or without interaction. The new formulation enables a great reduction of the computation load by eliminating possible redundant complexity involved in previous algorithms. A relaxation technique is designed for numerically detecting optimal paths involving entanglement. Interesting continuous symmetries are identified in the Lagrangian, which indicates the existence of physically equivalent classes of paths and may be utilized to remove neutral directions in the Jacobian of the evolution. In the ‘ground state’ solution among the set of optimal paths, the time-reversal symmetry of the system shows up, which is expected to be universal for the symmetry-related initial and final state.


Introduction
Ever since people realized that quantum computers are much more powerful than their classical counterparts [1,2], researchers have been working hard to bring it to life for years, which is making considerable progress and leads to a dramatic increase in the complexity of controllable quantum systems, from prototyped small quantum devices to machines possibly with thousands of qubits [3]. Undoubtedly, in order to reduce the impact of dissipation or decoherence that is always present in such complex systems, people continue to work on time-optimal control strategy, which minimizes the necessary traversing time required to reach a quantum target (namely, a quantum state or a unitary operation) and proves to be important for building efficient gates in quantum computer architectures. Therefore, a lot of research about the quantum speed limit (QSL) has been reported, which gives a lower bound of the time needed for a quantum system to travel from the initial to the final state and is estimated with the average energy or its variance [4][5][6][7]. In addition, time optimal control provides a physical framework [8] for defining complexity of quantum algorithms, probably superior to traditional ones [9]. For example Nielsen et al proposed a new quantum computing criterion by using the Riemannian geometry and transforming the problem of finding an optimal quantum trajectory into a geometric one in the space of Hamiltonians [10], and revealed that the complexity of a quantum gate is related to the problem of optimal control [11].
In a series of works on the QSL [4][5][6][7], the estimation of minimum evolution time is quite general and applicable to many different situations, in which, however, the optimal evolution path is not given explicitly. In the literature, the Grape algorithm and the Krotov algorithm are very popular for locating optimal evolution path. In the Grape algorithm [12,13], the initial and the target state evolve forward and backward respectively for the same amount of time to obtain two state sequences and the control variables are updated by iteratively minimizing the difference between the two sequences. The Grape algorithm relieves the limitation in the number of control variables and thus enjoys high flexibility. However, the convergence is slowing down when the deviation tend to a local. Also, an improper selection of step size often leads to a

The variational principle
Without loss of generality, the Hamiltonian of a quantum system that is of interest could be written as: where ξ j (t)'s are real numbers and A j 's are Hermitian operators that satisfy Tr A j = 0 and Tr(A j A k ) = δ jk (j, k m). The traceless Hamiltonian could be divided into two parts: the controllable part A = n j=1 ξ j (t)A j with each ξ j (j = 1, . . . , n) subject to external control and the uncontrollable part B = m j=n+1 ξ j (t)A j with ξ j 's being constant or changing in a pre-determined way, often used to model interaction between qubits. For later convenience we use u and v to label coefficients of the controllable part u = (ξ 1 , ξ 2 , . . . , ξ n , 0 , . . . , 0) t and the uncontrollable part v = (0 , 0 , . . . , ξ n+1 , ξ n+2 , . . . , ξ m ) t , and hence ξ = u + v. In this work, we build a formulation applicable to the general case which includes both components.

A variational scheme for the QBE
The action integral is written as [21]: (2) where P(t) = |ψ(t) ψ(t)| is the projection to the state |ψ(t) which is a normalized wavefunction defined in the configuration space, φ is an auxiliary wavefunction and c.c. denotes complex conjugates.
(ΔE) 2 ≡ ψ|H 2 |ψ − ψ|H|ψ 2 is the energy variance. Parameters λ and λ j are Lagrange multipliers and all {λ j } j=1,2,...,n are set to zero, ω is a given constant, the Q j 's denote the coupling constants in the Hamiltonian and the Planck's constant is chosen to be 1 for simplicity. Now we take the variation of the action equation (2) to obtain the equation of motion and all the constraints imposed to the variables.
(a). The variation with respect to φ gives which is the usual Schrödinger equation. In the projective space CP n−1 , the Fubini-Study line element is written as ds 2 = dψ|(1 − P)|dψ , which measures the displacement in the angular direction. Here, it is easy to check that ds 2 = (ΔE) 2 dt 2 . Hence, the first part of equation (2) refers to the total time of the process and the other parts embody differential or other constraints. (b). Here, we take the variation with respect to λ and λ j separately, resulting in 1 2 m j=1 ξ 2 j = ω 2 , which defines the size of the Hamiltonian as a finite energy constraint, and ξ j = Q j (n < j m), which depicts the influence of generic interactions.
(c). The variation with respect to ψ gives where · refers to the average with respect to |ψ . Equations (3) and (4) are seen in the literature [21] but the equations below appear new. (d). The variation with respect to H could now be taken as multivariate change with respect to {ξ j } j=1,2,...,m , resulting in where F jk = A j A k + A k A j − 2 A j A k , D j = φ|A j |ψ + c.c. and k stands for m k=1 which is also the case in all the following equations. The detailed derivation of equation (5) has been relegated to appendix A. Next, we derive the equations of motion for λ and {ξ j } j=1,...,n .
A multiplication of both sides of equation (5) with {ξ j } j=1,...,n , the controllable components, and then a summation over j results in where D · ξ = m j=1 D j ξ j and ω 2 = n j=1 ξ 2 j is constant (the constraint on the controllable part of the Hamiltonian). Noticing that {ξ j } j=n+1,...,m ≡ 0, we obtain: whereF = ξ · F andF j = k F jk ξ k . Here we used the fact that d dt D · ξ = 0 (details in appendix B), ξ · C = C · ξ = 0 (details in appendix C) and where C jk satisfies [H, A j ] = −i k C jk A k , with C jk being a real number. For a detailed derivation of equation (8), please check appendix D. Feeding equation (8) into equation (7) results iṅ Taking the time derivative of both sides of equation (5), and together with equation (8), we getu With equations (8), (9), (10) and the Schrödinger equation equation (3), the optimal path can be calculated now.
Here let's discuss a special case with H ∈ A, namely ξ = u, which means that all parts are controllable and thus v = 0. It is obvious that this case is simple since Hence equation (9) implies thatλ = 0. equation (10) could be reduced tȯ Considering u · C = 0 (appendix C) and feeding equation (11) into equation (13) leads to a trivial solutioṅ u = 0, which indicates that u is a constant vector and the evolution corresponds to a simple rotation with a fixed axis on the Bloch sphere.

Gauge symmetries in the Lagrangian
In this section, we discuss the gauge symmetries in the Lagrangian equation (2), which create neutral directions in the evolution and may bring complication in numerical calculation. First, consider a transform U 1 : φ → φ + iaψ, a ∈ R, with which the Lagrangian becomes The expression ψ |ψ + c.c is a total differentiation, which could be integrated out. With the Hermiticity of H, the extra integral in equation (14) vanishes and henceÛ 1 is a symmetry transform. Next, we consider another transformÛ 2 which only changes the values of the invisible Lagrange multipliers φ, λ, λ while the equation for the observables ψ, u remain intact. Actually, it is most obvious in the simplified formulation equation (18) below, where Ω, u,Ω,u keep invariant under the transformÛ 2 . HenceÛ 2 is a symmetry transform.

The evolution of Hamiltonian as a rotation
The equation ξ · C = 0 tell us that v · C = −u · C. Substituting the equation (5) into equation (10) results inu where I is the identity matrix and P u is the projection operator |u u| |u| 2 (ω 2 = n j=1 ξ 2 j = n j=1 u 2 j = |u| 2 ). Here we use the Dirac notation to emphasize the projection, whereas u is actually real. With the fact that v = 0, it is easy to check that ξ ·ξ = u ·u = 0. And the uncontrollable part of the time derivative of equation (5) similarly gives the evolution of λ where P v = |v v| |v| 2 and |v| 2 |u| 2 is a real constant. After, defining Ω = − 2 √ Tr I λ λ , the angular velocity of the rotation of u can be described with Ω. Furthermore, the QBE problem can be characterized with a simpler set of equationsu along with the Schrödinger equationψ = −iHψ.
According to appendix C, C can be conveniently determined, which in specific circumstances may provide possible routes to analytic solutions. As a specifical case, in the case without entanglement, the phase space can be divided into separate Bloch spheres of qubits. Then, according to equation (18), the equation of a single qubit in three-dimensional space is appreciably simplified tȯ u = Ω × ξ −û · (Ω × ξ)û = Ω × u. Meanwhile, the uncontrollable part is generally limited to no more than one dimension per qubit, or there will be no solution for most boundary conditions. As a result, Ω is in the direction of λ which is along v in the case of three dimensions. Here the cross product holds for each SO(3) Bloch sphere, which indicates that u is rotating with the angular velocity Ω.
As Ω, v is in the same direction for a single qubit. After multiplying Ω with equation (18) we see that |Ω|, the rotation speed of u, is a constant vector. In fact, in this case λ = 0 and thus Ω = 0,u = 0 as given in the previous section. It is easy to see that the spin vector σ (σ j={x,y,z} are the Pauli matrices) is rotating around u, and the evolution of wave function is analytically computed with rotation matrices.

A reduced scheme to the QBE
In practice, there are usually a large number of uncontrollable terms in the Hamiltonian. However, for the convenience of computation and design, most interaction terms are fixed. In this case the number of control parameters is appreciably low, while, the computation complexity of the conventional numerical methods remains as high. To address and eliminate the redundant complexity, here we propose an alternative approach, reintroducing the Lagrange multiplier φ previously replaced by D.
In general situations we see that the tremendous complexity is mainly brought in by the factor C · D, since the commutation relation (D6) will incur a large number of or sometimes even infinitely many A k 's and thus D k 's in the QBE equations when there are multiple interacting qubits. It is seen that most of the previous searching schemes [21,31] are confronted with this problem but no special attention has been paid. Here, we introduce an efficient new scheme based on the fact that C · D is immediately derivable from φ. Technically, consider then equation (5) turns to Thus C · D j = C · φ |A j |ψ + c.c.. Furthermore, it is much easier to evolve φ than φ since equation (4) yields which gives In analogy to previous derivations in section 3, it is easy to check that equation (18) now holds in an alternative formu where from equation (20), c., which could be evaluated directly with the computed states φ and ψ and the known H and A i 's. In conclusion, the whole QBE problem is fully described by equation (23). Terms like C · F · ξ/2ΔE 2 in C · D are eliminated now because of the introduction of the state variable φ , and hence the previous large number of equations needed to evaluate C · D are simply replaced by the evolution of φ . Since the number of degrees of freedom of vector φ is much smaller than that of matrix C. Now a two-point boundary value problem with equation (23) is left for us to solve. To begin with, we count the number of independent variables. The number of different variables in equation (23) is given by equation (20) (N u + N λ constraints), the energy constraint n j=1 ξ 2 j = ω 2 , the prescribed initial and final states ψ(0) = ψ 0 and ψ(T) = ψ target (2N ψ constraints). Thus we have N = N (where N λ = 1 and N φ = N ψ ). Here we evaluate the number of differential equations that we have to solve. In the work of reference [31], the number of equations is N u (accessible space) +N λ (forbidden space) +N ψ . In our method, this number is N u + N ψ + N φ + N λ . Very often, especially for multiple qubits N φ + N λ N λ . The difference increases with the number of qubits, as we will see in the following examples. Obviously, our method needs to solve fewer equations in complex situations, which is due to the introduction of the variable φ .
In order to steer the final state ψ(T) to the target state ψ target by adjusting the initial conditions and the flight time T, the famous shooting method (explained in appendix F) and the idea of 'relaxation' will be carefully used in finding the optimal path when Hamiltonian H contains interaction terms. We give more details about relaxation below.
When it comes to large-scale qubit systems, the shooting method per se can be too costly and need a good trial solution that is sufficiently close to the exact solution. In order to overcome those shortcomings we suggest a 'relaxation' method to meet the boundary conditions step by step. First, one stars from the case without interaction, solving the problem with separate qubits. With interaction present, secondly, we invite part of the interactions back and use the no-interaction solution as the initial condition (i.e. taking the first-step coupling constant J = 0.01 in our examples), so a new solution in case of the weak interaction is obtained, which can be used as the initial condition for the case with a little bit stronger interaction. Repeat this process until the full interaction is restored and an optimal path is then obtained. This is a relaxation process since the effective Hamiltonian in the computation is gradually relaxed to the desired one.

One-qubit Bloch sphere
Here we first try to search for local optimal paths of single-qubit evolution on the Bloch sphere. The goal is to flip the spin state of a single Fermion from | + x to | − x .
Assuming = 1, the Hamiltonian with the full control is where γ = g q 2m is the gyromagnetic ratio and q is the charge carried by the particle. Since γ −2, then H B · σ.
Based on previous discussions section 3, the evolution path of ψ can be analytically addressed, as a rotation with the angular velocity ω eff around an axis u rotating itself at Ω. The effective angular velocity ω eff of spin vector is determined as where B is an effective magnetic field and the Hamiltonian H = σ · B = j A j u j . This is derived from the fact that the SU(2) element of a rotation with Δφ is e iσ Δφ 2 . In the frame of reference rotating with Ω, the motion of ψ is viewed as a fixed-axis rotation with the angular velocity Ω = ω eff − Ω.
For instance, consider the numerically simulated case [21] where . Now the evolution turns out to be simply a rotation around Ω followed by a reversed rotation around Ω back to the static frame. The rotation matrix Then it is easy to check that where the values of Ω and Ω are twice as much as those in reference [21]. Equation (26) gives a family of analytic solutions if we choose (ΩT, Ω T, Ω ω eff ) = (kπ, lπ, k √ l 2 −k 2 ), where l, k are non-negative integers and l + k is odd. The case k = 0 corresponds to the geodesic solution with no anisotropic constraints and cases of k = 0 correspond to the oscillating solutions with l − 1 nodes. We will denote L = l − 1 for later convenience. The oscillating curve in figure 1 plots the solution with k = 1, L = 0 and 1.
The simulated paths obtained by our numerical method are consistent with the analytical solution and displaced as well in figure 1. From previous discussions, we know that in this caseu = 0, i.e.Ḃ = 0. So the optimal path is naturally a geodesic curve (see figure 1 (left), which shows one of the many geodesic solutions). The red line on the equatorial plane (marked by the green grid inside the sphere) indicates the evolution of the external magnetic field B, which does not change throughout the evolution. In principle, the constraints Tr (Hσ j ) j=x,y,z = 0 allow oscillating solutions withḂ = 0. Let's take the case of H = B x σ x + B y σ y , B z ≡ 0 as an example. In numerical computation, we use the shooting method to identify the initial value φ 0 for the prescribed boundary conditions ψ(0) = | + x , ψ(T) = | − x . Figure 1 (right) shows the resulting optimal path which oscillates with one node on the equator, where Ω = (0, 0, − 2 √ 3 ). The evolution of B (the area swept by the red line in figure 1 (right)) corroborates the theory of rotation in section 3. Of course, it is not hard to get a longer locally optimal path with more than one node, which we will not show here.

Two-qubit Bloch sphere
In the case of two qubits without interaction, two Bloch spheres may be used to represent a non-entangled state and the geodesic solution is simple according to the discussion in section 5.1 above. Here, we focus on the oscillating solutions under the anisotropic constraint featuring B z = 0. We give an example with Figure 2. An optimal path with L = 1, ω (1) eff = 2, ω (2) eff = 2 connecting the initial state | + x (1) , | + y (2) to the final state | − x (1) , | − y (2) of qubits 1 and 2 plotted on two Bloch spheres with the Hamiltonian given by equation (27) and T = 2.7202.
Similarly, our purpose is to flip both qubits (where ψ(0) (1) = | + x and ψ (1) target = | − x for qubit one, ψ(0) (2) = | + y and ψ (2) target = | − y for qubit two). Obviously, these two qubits are both reversing their states, and because x-and y-directions are totally equivalent these qubits should share the same energy at the ground state (namely ω (1) eff = ω (2) eff ). Figure 2 shows a local optimal solution of the constrained two-qubit system with L = 1. Ω (1) . It is apparent that the evolution of B turns out to be a rotation as expected. (Since the numerical solution is consistent with the analytical one, it will not be distinguished in the following. Also, we can obtain longer paths with larger L, not shown here.) Actually, the local solutions can be derived analytically as long as the interaction terms are missing in the Hamiltonian. According to the rotation picture, the analytical approach can be applied to no-interacting many-qubit systems with an arbitrarily large number of qubits. Here in figure 2 the solution corresponds to (k (1) , L (1) ) = (k (2) , L (2) ) = (1, 1) as the oscillating solution in figure 1 (right), so they naturally share the same geometric property. Later on, these solutions would be good initial guesses in the presence of interaction in the searching algorithm based on the idea of relaxation (see 4).

Two or three qubits with interaction
In the case of a single qubit or double qubits without interaction, it is easy to obtain the solutions analytically, which are very regular. With interaction, however, it is very difficult to obtain such an analytical solution, and thus we resort to the numerical method discussed in appendix F. In this section, three different cases are investigated together with a general description of the consequences brought by possible entanglement in the wave functions.
Case 1, the initial and final spin states are the same as in section 5.

2, but there is an extra interaction term in the Hamiltonian
where J yy = 1 is a fixed interaction constant. In the numerical computation, in order to utilize a reasonable initial condition, we gradually increase the value of J yy from 0.01 to 1 with a step-size 0.01. After the interaction term J yy is introduced, the wave function ψ(t) is no longer separable, but becomes entangled. Following the work in the literature [32,33], we may define an entanglement index E p for the two subsystems A and B, where ρ A = Tr B (|ψ AB ψ|) is the reduced density matrix of the pure state |ψ AB over subsystem A, and ρ B has a similar definition. Figures 3(a)-(d) depict the entanglement of ψ(t) with L = 1, 3, 2, 4, which have a common feature: the entanglement of the wave function connecting the two separable states (the initial and the target state) increases first, reaches a maximum and then decreases back to zero. The entanglement evolution profile for L = 1, 3 (figures 3(a) and (b)) is symmetric and the maximum is exactly at the symmetry point, while for L = 2, 4 (figures 3(c) and (d)), the profile is skewed. Given a suitable initial value of φ 0 , the computation converges exponentially when approaching an optimal path (See figure 5(c)). Next, we keep the initial state unchanged and set the target state as an entangled one with E p = 0.6061. The wavefunction ψ(t) connecting the two states is computed with J = 1, 0.9, 0.65, the entanglement evolution  of which is depicted in figure 4. We observe that with the decrease of the coupling strength J, the change of entanglement E p (t) tends to be a linear curve. In addition, we observe that there is an approximate center symmetry about t = 0.5T, indicating a fast initial rise and final relaxation but a steady transition between. Case 2, we consider a two-qubit model with the following Hamiltonian where σ (l) m with (m = x, y, z, l = 1, 2) are the Pauli matrices for the lth qubit. The initial and the final spin states are the same as in section 5.2 and the coupling constant J = 1. The results obtained with our method Figure 5. The evolution and the convergence with the Hamiltonian given by equation (30). The evolution of the measurement σ(t) (a) and the control magnetic field shows that the control magnetic field B(t). In this case, N u = 6, N λ = 9, N φ = N ψ = 8, N λ = 1, the number of QBEs that we solve is 23 both in the reference [31] and in our method.
Next, we make a comparison of the efficiency of finding such a connecting wavefunction among different methods with the same initial conditions. The convergence of different methods is compared in figure 5(c). Our method (scheme I, red points), which though oscillates at the initial stage (in the first dozen iteration), converges much more quickly (within less than 30 iterations) and displays good stability as well. Scheme II and III denote the Krotov algorithm [13][14][15][16] and the Grape algorithm [12,13], respectively. The error of these two algorithms are always monotonously decreasing. For the Krotov algorithm, the error drops very fast at the beginning (the Krotov scheme shows better convergence for short times), but the decrease slows down after 10 iterations and hence takes long to reach a small value. For the Grape algorithm, we clearly see that the convergence is quite slow (with more than 500 iterations), and the rapid convergence is not achieved until the iteration proceeds many steps. From the perspective of energy, as shown in figure 5(d), our method ensure that the strength of the magnetic field remains constant (|B(t)| = 2.23 607, red line) during the evolution, while in the Krotov algorithm (black line) and the Grape algorithm (blue line) the magnitude of the magnetic field changes over time. The optimal time is computed as T =

√
Nπ ω for a quantum state evolving between two end states of a diameter on the Bloch sphere, which is consistent with the lower bound computed by Margolus and Levitin [34,35]. The Fubini-Study distance to the target state seems to decrease linearly with time as shown in figure 6, and we note that with the decrease of the energy ω, the traveling time T increases. Therefore, because the initial state ψ 0 and target state ψ target are given, the total Fubini-Study distance of the two states is consistent. Our method ensures that the energy of the system is constant in the process of evolution. Larger energy corresponds to faster movement, so the larger the energy, the faster the F(t) decreases. The linear dependence of F(t) on time indicates a uniform approach to the target state on the optimal path. Figure 6. The Fubini-Study distance F(t) between ψ(t) and ψ target with ω = 2.8284, 2.7803, 2.6926, 2.5554, 2.3601 and the Hamiltonian given by equation (30). The distance between neighboring states along a trajectory is defined by using the Fubini-Study line element ds = dψ|(1 − P)|dψ , where dψ is the change of the state along the path and P = |ψ ψ| is a projection operator. Re(ψ(t)) and Im(ψ(t)) represent real part (a) and imaginary part (b) of wave function ψ(t) respectively. The notation |000 ∼ |111 in the legend represents the eight components of ψ(t).

Case 3, we study a model with three qubits for which the Hamiltonian is
where σ (l) m with (m = x, y, z, l = 1, 2, 3) are the Pauli matrices for the lth qubit. We choose the initial state (3) . The obtained optimal wavefunction ψ(t) connecting them with the interaction strength J = 0.2 and L = 1 is depicted in figure 7. Figure 7(a) depicts the real and (b) the imaginary part of ψ(t). In this case, N u = 9, N λ = 54, N φ = N ψ = 16, N λ = 1, the number of quantum Brachistochrone equations that have to be solved is 79 with the method in reference [31], however, the number is only 42 for our method. At this time, it is obvious that the number of equations we need to solve is much less, especially with the increase of the qubit number.

A symmetry operation
In this section, we check a symmetry operation of the Hamiltonian equation (28) and ferret out possible mechanism underlying the phenomenon discussed in section 5.3 above. It is interesting to observe a discrete symmetry operation in Hamiltonian (28) associated with x and B (2) y have a reflection symmetry with respect to t = 0.5 and B (1) y and B (2) x have a center symmetry with respect to (0.5, 0), or vice versa. All the above indicates that the maximum of E p is associated with a fixed point of the operation F for the case L = 3.
In the cases L = 2, 4, the symmetry in ψ(t) and B(t) discussed previously disappear. Taking L = 2 for an example, figures 9(a) and (b) do not coincide with figures 9(c) and (d) and the maximum value of E p is not at t = 0.5. On the other hand, under the discrete symmetry operation F, we have B (1) y . It's still true for L = 2, 4 (see figures 9(a)-(d)). When we overlap graphs figures 9(a) and (c), we find that they are axisymmetric with respect to t = 0.5. As the symmetry partner ψ t of the wave function ψ t is obtained with the symmetry operation F, even if the entanglement evolution of ψ (t) or ψ(t) (see figure 10(a)) itself is asymmetric, they are axially symmetric with respect to each other about t = 0.5, which also indirectly reflects the symmetry of the corresponding Hamiltonian equation (28).
In the case that ψ 0 and ψ f do not flip to each other, e.g. the angle between ψ 0 and ψ f is less than π, we then have Fψ 0 = −ψ f (take the case L = 1 as an example), which indicates that the entanglement evolution cannot be symmetric (shown in figure 10(b)). However, the maximum entanglement still exists which is  unique but not at t = 0.5, and also from the figure we see that this maximum increases with the angle between the initial and the final state.

Summary and discussion
In this work, we derive a general set of differential equations (equation (23)) for the optimal quantum control of single or multiple qubits with or without interaction. In the derivation and application of equation (23), we discussed the following aspects: first, previous numerical methods conventionally involve a large number of physically irrelevant variables derived from the commutation relation with system Hamiltonians. We thus introduce a new set of observables to eliminate the physically irrelevant part, whereby an accelerated computation becomes possible. As a result, the efficiency of computation are boosted significantly. With the help of a 'relaxation' idea, the drawback of the shooting method, of requiring a good initial guess, is surpassed by slowly pushing the trial solution to the correct one with a gradual restoration of the interaction. Second, we discussed symmetries of the Lagrangian formalism of the QBE. In the numerical calculation, these entail unwanted degenerate directions in the Jacobian matrix, which may cause serious trouble in the solution of the resulting boundary value problem. Finally, groups of analytic solutions are given as rotations of qubits on Bloch spheres in the absence of interaction, which may be used as the starting point of the above relaxation scheme.
In an application of the new scheme to the case of two qubits with interaction, the solution to the QBE displays a unimodal evolution profile of the entanglement. If the initial and final state transform to each other under a symmetry operation of the Hamiltonian. The evolution of the optimal path either bears a related symmetry or has a symmetry partner. More explicitly, the solutions with even number of oscillations in the profile are asymmetric and the symmetry operation gives another optimal path-the symmetry partner of the original solution.
With the current formalism, we are able to design optimal control strategy for multi-qubits with possibly complex interaction. However, the evolution is not unique and many local optimal paths exist which could be indexed by the number of oscillations of the paths in simple cases. It could be interesting to classify these paths and check their bifurcation routes when system parameters change. On the other hand, this work may provide experimental researchers a practical computation tool for quantum information regulation, since an optimal control strategy of the magnetic field is readily designed based on our scheme.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix A. The demonstration of equation (5)
The variation of the Hamiltonian H could now be taken as multivariate change with respect to {ξ j } j=1,2,...,m . Let's respectively check different parts of the action integral equation (2). First, check one term of the integrand where D j = φ|A j |ψ + c.c. are real numbers, thus the Hamiltonian variation of the term I is Next, let's consider the second part Let G j = δg(ξ,λ) δξ j and we obtain the variation of II Third, before treating first term of the equation (2), we need to do some preparation. Here, we define: which is a symmetric tensor, and the energy variance could be written as: leading to which gives So, the variation with respect to ξ j of the first term of equation (2) gives where δ ψ |(1 − P)|ψ = 0, because here we only take the variation with respect to ξ j , and where (A5) has been used. With equation (A2), (A3) and (A9), the variation with respect to ξ j gives where G j = λξ j , ( j n).
Appendix C. The demonstration of C · ξ = ξ · C = 0 Here we demonstrate that C · ξ = ξ · C = 0 is universally valid with respect to spin algebra. In the case of H ∈ A i.e. ξ = u, it yields C · u = u · C = 0. and I is 2 q dimensional unit matrix (q is the number of qubits) and A k = 1 √ 2 σ k for q = 1. Hence Tr I njk ξ n A k which yields C jk = n − 1 √ Tr I njk ξ n . Thus ξ · C = −C · ξ = 0.

C.2. Many qubits with and without entanglement
We will prove that ξ · C = C · ξ = 0 remains valid. Below, we denote the unit matrix I 0 with σ 0 and extend our index to include 0.
Lemma. C jk = i Γ ijk ξ i , where Γ is an antisymmetric tensor.
Proof. First let us check some properties of Pauli matrices. With the commutation and anti-commutation relations {σ i , σ j } = 2δ ij I(i, j = 1, 2, 3), it is easy to check even if i, j = 0, 1, 2, 3 where σ 0 = I 0 which is two-dimensional unit matrix and α is the order index of qubits,