# Benchmarking LLM powered Chatbots: Methods and Metrics

Debarag Banerjee, Pooja Singh, Arjun Avadhanam, Saksham Srivastava

###### Abstract

Abstract

Autonomous conversational agents, i.e. chatbots, are becoming an increasingly
common mechanism for enterprises to provide support to customers and partners.
In order to rate chatbots, especially ones powered by Generative AI tools like
Large Language Models (LLMs) we need to be able to accurately assess their
performance. This is where chatbot benchmarking becomes important. In this
paper, we propose the use of a novel benchmark that we call the E2E (End to
End) benchmark, and show how the E2E benchmark can be used to evaluate
accuracy and usefulness of the answers provided by chatbots, especially ones
powered by LLMs. We evaluate an example chatbot at different levels of
sophistication based on both our E2E benchmark, as well as other available
metrics commonly used in the state of art, and observe that the proposed
benchmark show better results compared to others. In addition, while some
metrics proved to be unpredictable, the metric associated with the E2E
benchmark, which uses cosine similarity performed well in evaluating chatbots.
The performance of our best models shows that there are several benefits of
using the cosine similarity score as a metric in the E2E benchmark.

##  I Introduction

“A chatbot is a computer program that uses artificial intelligence (AI) and
natural language processing (NLP) to understand customer questions and
automate responses to them, simulating human conversation. – (”What is a
chatbot?”) As the field of AI continues to grow rapidly, chatbots have become
more and more refined, now capable of learning from large volumes of data and
improving their conversational abilities over time. This has given rise to a
new section of chatbots known as ML (Machine Learning) based chatbots that are
not pre-programmed with rules, but rather learn from large amounts of data and
respond to user input by analysing said data. A good guide on ML based
chatbots can be found at (IBM Watson Advertising,“The ultimate guide to
machine-learning chatbots and conversational AI”,2022).Advertising  Today
chatbots are being used in various domains such as customer service,
healthcare, and even education brandtzaeg2017people (Brandtzaeg, 2017). A
specific common chatbot application is one where they are used to provide
product support of an enterprise’ products, thus requiring less human
involvement in customer support. The foundation of many modern chatbots is the
Large Language Model (LLM), a model that processes natural language input and
generates words based on the data seen. Their impact and capabilities are
discussed in the paper tamkin2021understanding (Alex Tamkin,“Understanding the
Capabilities, Limitations, and Societal Impact of Large Language
Models”,2021). While LLMs have made chatbots very human-like in their customer
interactions, they also gave rise to a new problem known as Hallucination
(Nabeel Ahmed, Feb 2022). Hallucination is a phenomenon where an AI model
creates results that are not real, do not match any data the algorithm has
been trained on, or do not follow any other discernible pattern. It is thius
important to employ a combination of evaluation techniques, to assess the
presence of hallucination in large language models, especially when used in
the context of chatbots. With the increasing popularity of chatbots, it is
important to develop methods to evaluate their performance. A chatbot
benchmark is a standardized evaluation framework used to assess the
performance and capabilities of chatbot systems. It involves defining a set of
tasks or criteria that the chatbot must fulfil, and then measuring its
performance against those tasks or criteria. There are two important factors
to consider while benchmarking a chatbot: accuracy and usefulness. Chatbot
accuracy is a measure of how factually correct it is, while usefulness
measures to what extent the chatbot meets the user’s needs. Not that while
earlier an apparently useful answer would almost always be useful, with the
deployment of hallucination-prone LLM-powered chatbots, that is no longer the
case -i.e. it is possible for a chatbot to hallucinate up an answer that is
not supported by any of the available documents, and can be completely wrong.
Thus to successfully evaluate the chatbot’s performance, both these factors
need to be considered. There are various types of chatbot benchmarks, each
with their own sets of advantages and disadvantages. Some examples of such
benchmarks are discussed in wang2022modern (Wang, 2022). An example of a
category of chatbot benchmarks are those that compare the answer to a question
from a chatbot with that generated by a human. While there are several other
benchmarks such as Summarisation benchmarks, Information Retrieval benchmarks,
and more, such benchmarks provide the best option for a user-centric
evaluation, as it emphasizes the user’s perspective. Inspired by this
category, we propose the End-to-End (E2E) Chatbot Benchmark which compares the
answer provided by the chatbot to an expert human answer (also known as a
golden answer) using certain semantic similarity metrics that try to compare
the “meanings” of each answer, as opposed to just comparing their exact words
which is the current state-of-the-art (SOTA) in benchmarking. The E2E
benchmark also measures all the optimization factors previously mentioned for
a successful evaluation. Overall, the E2E benchmark provides a user-centric,
standardized approach and promotes real-world relevance. However, it does need
a set of Golden answers. In this paper, we propose the semantic similarity
aspect of the E2E benchmark to be measured using the cosine similarity measure
among their corresponding text embeddings. This paper discusses the E2E
Chatbot Benchmar, and how it improves the assessment and evaluation of
chatbots, compared to another well-used SOTA alternative, ROGUE scores.

##  II Related Work

Chatbot output testing involves giving inputs to the chatbot and analyzing the
output given by the chatbot. The analysis is done to gain insight about the
chatbot and to measure various relevant metrics like the relevance,
completeness, accuracy, and recall. Various metrics have been used by
researchers to evaluate the performance of chatbots which include but are not
limited to the use of benchmarks and metrics like GLUE, BLEU, perplexity,
ROUGE score and human evaluation. In the subsequent section we shall explore
the various metrics. Chatbots should be accurate, relevant, and useful for
them to satisfy user needs. Accuracy is indicative of how many correct answers
the chatbot is able to provide, relevance is usually measured with the help of
precision and recall. Perplexity is a commonly used intrinsic metric for
evaluation of any language model, it is used to evaluate how well a model
predicts text samples, that is- it is used to measure the degree of
uncertainty when a model generates a new token therefore perplexity scores
having lower values are indicative of a better performance by the model. Alex
Wang et al. (2019) in their conference paper published at ICLR presented
another benchmark- the General Language Understanding Evaluation (GLUE)
benchmark for language model evaluation which is used to evaluate the general
language understanding capabilities of a language model. It consists of a
collection of diverse tasks for NLU (Natural Language Understanding) along
with its own datasets and evaluation metric, it attempts to provide a single
framework for assessing model across various tasks however criticism about it
not encompassing the entire range of tasks, lack of a diverse dataset and lack
of multi-lingual assessments have been made. Bilingual Evaluation
Understudy(BLEU) [Kishore Papineni et al.(2002)] is another metric commonly
used for evaluation of machine translation of a text from one language to
another. BLEU determines the accuracy of individual n-grams by tallying the
count of matching n-grams between the candidate and reference translations.
Furthermore, it takes into account the conciseness of the candidate
translation relative to the references. These precision scores are
harmoniously combined through a weighted geometric mean, wherein the brevity
penalty plays a role, resulting in the ultimate BLEU score. BLEU was one of
the first metrics to claim a high correlation with human judgement, and still
is one of the most popular and inexpensive metrics for translation evaluation,
however it is not comparable across different datasets and does not take into
account intelligibility or grammatical correctness. ROUGE or Recall Oriented
Understudy for Gisting Evaluation is used for evaluating the text summarizing
capabilities of the model wherein the goal is to generate a concise summary of
a longer text. ROUGE score measures the similarity between the machine-
generated summary and the reference summaries using overlapping n-grams, word
sequences that appear in both the machine-generated summary and the reference
summaries. The most common n-grams used are unigrams, bigrams, and trigrams.
ROUGE score is important as it gives flexibility to use n-grams based on
requirements, however it may not fully capture semantic meaning or coherence
of summary, which is important for quality of the summary. Human evaluation is
also a very important criterion that any language model or chatbot should
excel in, however it can be very subjective and prone to biases, also large
scale human evaluation can be very expensive and time consuming.

##  III Benchmarks in Relation to Chatbots

This figure gives an insight into a typical chatbot designed for information
retrieval from a knowledge graph to user’s natural queries that are
(optionally) subsequently summarized by an LLM. It also provides different
benchmarks that can be constructed to tap into the different stages

Figure 1: A typical LLM-powered chatbot for answering questions based on a
document corpus and the various benchmarks that can be used to evaluate it.

The E2E benchmark uses a set of “Golden Answers” to accurately measure the
chatbot’s performance. Golden Answers are a set of pre-defined answers to
which the chatbot’s answers are compared. Through this process, we can
evaluate the chatbot’s performance. Golden answers also come with their own
set of advantages and disadvantages. These answers need to be defined before
runtime (manually) and thus can prove to be tedious, time-consuming, and do
not provide real-time feedback. However, several advantages make golden
answers an extremely useful evaluation method such as consistency,
standardized evaluation, and providing vital information for improving chatbot
performance. The E2E benchmark aims to simulate real-life scenarios and assess
the chatbot’s performance in handling diverse user queries. A combination of
objectivity and subjectivity is created by mixing metrics such as accuracy,
precision, etc., and human evaluators that provide subjective ratings and
feedback. Thus, the benchmark provides an experience based around the user, as
they are the ones that assess the chatbot’s performance based on criteria such
as relevance, correctness, etc. Through these methods the E2E benchmark
considers both factors involved in chatbot benchmarking, accuracy, and
usefulness and successfully measures all optimizations. Other benchmarking
methods have their advantages and disadvantages. The Relevance benchmark also
creates a user-centric experience and does not require the aid of golden
answers; however, it does not account for hallucination, which is the biggest
drawback of Large Language Models. The Summarization benchmark does not suffer
from hallucination issues and does not require golden answers, but it has a
narrow focus and is not accurate in measuring the important aspects of LLM
efficacy, such as language understanding, response generation, or context
retention. The previous two benchmarks are also unsupervised learning methods.
The Information Retrieval benchmark serves as a search benchmark rather than a
chatbot benchmark. Considering all these factors, the E2E benchmark proves to
be the best way to evaluate chatbot performance in this scenario. In this
paper we are introducing an additional technique to E2E benchmarking -beyond
word / n-gram based metrics like ROGUE and BLEU. Namely we are proposing that
if we measure the semantic similarity between golden answers and chatbot
generated answers, they constitute a more sensitive metric to measure any
improvements or degradations in the chatbot system. Specifically in this paper
we compare three different benchmarks, the first two -one using the Universal
Sentence Encoder (USE) (Daniel Cer, 2018)cer2018universal  as a semantic
embedding, another using Sentence Transformers (ST) (Nils Reimers,
2019)reimers2019sentence  as a semantic embedding, and we compare them with
the various components of ROGUE as a legacy word or n-gram based metric. To
conduct the analysis, three libraries have been used, namely:
1\. Universal Sentence Encoder
2\. Sentence Transformer
3\. ROUGE (Recall-Oriented Understudy for Gisting Evaluation)
As previously discussed, this paper will primarily be focusing on the E2E
benchmark as compared to other methods. In the following sections of the
paper, the results of this evaluation will be discussed.

##  IV Results

We propose here a semantic-similarity based End-to-End (E2E) Benchmark for
benchmarking performance of chatbots. This benchmarking compares and checks
the answers from the chatbot against a ‘golden answer’ provided by a human
expert using a semantic similarity measure based on text embeddings. We also
compare the result we get from the such-constructed E2E Benchmark vs a more
traditional approach of comparing the word, bi-gram, and n-gram based ROGUE
scores computed on the same pairs of ‘golden answer’ vs chatbot answers.
We used three libraries to conduct our analyses which included
1\. Universal Sentence Encoder
2\. Sentence Transformer
3\. ROUGE
Before continuing the analysis, we would like to give a brief idea about
Universal Sentence Encoder and Sentence Transformers.

Figure 2: E2E Benchmark results with use of Cosine Similarity over embeddings
obtained from the Golden Answer vs the Chatbot generated Answer -using both
Sentence Transformer (ST) and Universal Sentence Encoder (USE) embeddings,
respectively. Result obtained from a product support chatbot.

Universal Sentence Encoder (USE) is a pre-trained model developed by Google
that converts sentences into fixed-length vector representations, or
embeddings, capturing their semantic meanings.

The Sentence Transformer (ST) is a model architecture designed to produce
high-quality sentence embeddings, which represent sentences as vectors. Unlike
traditional word embeddings, these sentence embeddings consider the complete
context and meaning of the sentences. Drawing inspiration from transformer
architectures like BERT, the Sentence Transformer model utilizes pre-training
and fine-tuning techniques.

For our purposes, we evaluated the performance of a publicly available product
support chatbot made by a large network equipment company to support their
enterprise users get answers to questions based on their documentation.

For both the embeddings, we take a golden answer (G) and the chatbot generated
answer (A), the we apply the embedding $X=F(T)$ where F can be either a
Universal Sentence Encoder or a Sentence Transformer on our text T- which
gives us $X_{G}$ and $X_{A}$ both being n-dimensional vectors generated as a
result of applying the operations. Then we calculate the cosine similarities
between $X_{G}$ and $X_{A}$ . Cosine similarity is defined as dot product of
two vectors, divided by the product of their magnitudes.

$S_{(G,A)}=\frac{X_{G}\cdot X_{A}}{\left|X_{G}\right|\left|X_{A}\right|}$

In case of using USE as the embedding operation, there is usually a constant
bias of 0.5 remaining in our final answer which is due to the inherent
characteristics of USE vectors. We computed the cosine similarities for both
Universal Sentence Encoder and Sentence Transformer, and for ROUGE we use
precision and recall to check for words and n-grams included in our generated
answers while comparing against the expert human answers. As shown in the
Figure 2, the chatbot performs reasonably well for both Universal Sentence
Encoder and Sentence Transformer, with the USE mean of about 0.64 and the ST
mean of about 0.47, we also calculated the correlation between the outputs of
both the methods which is displayed in Figure 3.

Figure 3: Correlation between E2E Benchmark results with use of Cosine
Similarity over Sentence Transformer (ST) embeddings, and Universal Sentence
Encoder (USE) embeddings from the same product support chatbot

We calculate the $R^{2}$ scores between the cosine similarity of both the
Sentence Transformer and Universal Sentence Encoder to gauge the similarity
between the models for our chatbot which gives us an impressive 0.7 which
indicates that the models both work pretty similar and have very similar
answers to our provided expert answers. Next, we turn our attention to ROUGE
scores. ROUGE-1 or word-level ROUGE score which evaluates the quality of
answer based on the overlapping of unigrams between the answers provided by
the chatbot, ROUGE-1 is useful as it captures the keywords however it fails to
take into account the order or structure of the words. We can also use the
metric (1-precision) to get an estimate of the hallucination the chatbot has.
The chart below provides the ROUGE-1 scores when we tested our chatbot answers
against the expert answer:

Figure 4: Precision vs Recall for ROGUE-1 scores between Golden Answers vs
Chatbot Answers from the same product support chatbot

The results are very interesting. They don’t show any direct mapping between
recall vs precision at unigram level however we also can use other forms of
ROUGE to get a better understanding of the model- ROUGE-2 (bigrams) and ROUGE-
LCS (longest common sequence)- the recall vs precision of which are given in
Figure 5 and Figure 6, respectively.

Figure 5: Precision vs Recall for ROGUE-2 scores between Golden Answers vs
Chatbot Answers from the same product support chatbot Figure 6: Precision vs
Recall for ROGUE-1 scores between Golden Answers vs Chatbot Answers from the
same product support chatbot

For bigrams, we see most of recall improve compared to unigram level,
precision however was better when we used ROUGE-1 as it only captures one word
at a time.

##  V Golden Answers vs. Random Words : A Negative Example

We also calculated the cosine similarity as described in the above sections-
this time between the golden answers and random words in order to better judge
how well our experiments perform in relation to a random sample, we found out
that the metrics are completely uncorrelated, with the R2 score being 0,
however we also notice that Sentence Transformer outperforms USE here as USE
still gets a cosine similarity of about 0.5 whereas Sentence Transformers give
a cosine similarity of 0. Image below shows the cosine scores calculated
between golden answers and random words.

Figure 7: E2E Benchmark results with use of Cosine Similarity over embeddings
obtained from the Golden Answer vs random words -using both Sentence
Transformer (ST) and Universal Sentence Encoder (USE) embeddings,
respectively.

We even calculated the correlation between the metrics as mentioned above the
results of which are in Figure 7.

Figure 8: Correlation between E2E Benchmark results with use of Cosine
Similarity over Sentence Transformer (ST) embeddings, and Universal Sentence
Encoder (USE) embeddings comparing Golden Answers vs random words

As we expected, there is no correlation between the metrics in this case.

##  VI Measuring improvement in the chatbot: GPT-3 with prompt engineering

It has been observed in the literature that LLM-based chatbots can be improved
using enhanced prompts, what is often called prompt engineering The chatbot we
used for our experiments uses GPT-3, and as such is sensitive to prompt
engineering. So we compare the E2E Benchmark results obtained using standard
prompts vs the results obtained after using enhanced prompts- using the
methods described :
1.) Standard prompt : ‘from this text: , can you answer the question: ’
2.) Enhanced prompt: ”role”: ”system”, ”content”: ”You are a network
engineering expert helping summarize the following message to help a question
from the network operator which follows the message”
The E2E results from these experiments are in Figures 9 and 10 respectively.

Figure 9: E2E Benchmark results with use of Cosine Similarity over embeddings
obtained from the Golden Answer vs random words -using both Sentence
Transformer (ST) and Universal Sentence Encoder (USE) embeddings,
respectively, using the Standard Prompt Figure 10: E2E Benchmark results with
use of Cosine Similarity over embeddings obtained from the Golden Answer vs
random words -using both Sentence Transformer (ST) and Universal Sentence
Encoder (USE) embeddings, respectively, using the Enhanced Prompt.

We see a significant rise in Sentence Transformer outputs on engineered
prompts, suggesting it is much more sensitive to prompt engineering tweaks on
the query. We also checked the performance of the chatbot against ROUGE
metrics, with Standard prompts, vs Enhanced prompts. The ROGUE-1 score results
from these experiments are in Figures 11 and 12 respectively, while the
ROGUE-2 score results from these experiments are in Figures 13 and 14
respectively.

Figure 11: ROGUE-1 score comparisons between the Golden Answer vs Chatbot
Answers, using the Standard Prompt. Figure 12: ROGUE-1 score comparisons
between the Golden Answer vs Chatbot Answers, using the Standard Prompt.

We see that for most outputs, the improvement in precision scores do not
correlate with increase in recall however ROUGE-2 capture this trend better as
we see in Figures 13 and 14.

Figure 13: ROGUE-2 score comparisons between the Golden Answer vs Chatbot
Answers, using the Standard Prompt.

When given engineered prompts- ROUGE-2 shows rather interesting output as
shown below-

Figure 14: ROGUE-2 score comparisons between the Golden Answer vs Chatbot
Answers, using the Standard Prompt.

ROUGE-2 shows clustered outputs, however the maximum values and the number of
outputs close to them seems to dip for both the axes. ROUGE scores do not show
improvement on applying prompt engineering. However the E2E benchmark metrics,
i.e. cosine similarity of the USE and ST embeddings show considerable
improvement on using engineered prompts.

##  VII Conclusion

Over the past few years, there has been tremendous growth in the number and
quality of chatbots. To ensure that these chatbots are up to the mark we need
to be able to effectively assess their performance. This has created the need
for efficient chatbot benchmarking methods. Throughout this paper, many
methods of chatbot benchmarking have been discussed, highlighting their
various advantages and disadvantages. In addition we have proposed a new
benchmark that we cal the E2E (End-to-End) Benchmark. The E2E Benchmark is
constructed using cosine similarity between sentence embeddings of the
chatbot’s answers to questions vs the corresponding “golden answers” to the
same questions generated by human experts. For generating such embeddings we
considered two different embedding libraries, the Universal Sentence Encoder
(USE) and the Sentence Transformer (ST). While the E2E Benchmark using USE
provides a higher mean score than the same using ST, the chatbot performed
well for both libraries and through the plotting their scores against each
other, we observed a similarity in their scores. Replacing the golden answers
by a set of random words should have resulted in a mean score of near 0 which
was observed in the case of ST, but USE showed a mean score of 0.5 which was
likely due to its systemic bias. By providing enhanced prompts to the chatbot,
we observed an increase in mean score for both libraries, however the ST score
showed a sharper increase, indicating high sensitivity to the enhanced prompt.
The ROGUE scores were also computed using two different metrics; precision
which tells how concise the answer is and recall, which tells us how accurate
the answer is when compared to the Golden answer. Different ROGUE scores are
produced when run on different n-gram levels (e.g., unigram-ROGUE-1). The
ROGUE-1 score did not show any obvious correlation between precision and
recall. For ROGUE-2, we observed an improvement in recall however, the
precision score was not as good compared to ROGUE-1. Applying Prompt
Engineering on the LLM part of the chatbot, we did not significantly improve
the ROGUE scores -while the E2E Benchmark did show significant improvement.
Through our analysis, we concluded that the E2E Benchmark, when it is used
with the ST library, works best as it performed well in all scenarios and
showed significant improvements at each stage. E2E Benchmarks using USE
followed that as the second best. While the ROGUE scores were too
unpredictable, and as a result, we can’t recommend them as a reliable means to
benchmark chatbots’ performance. In summary, while chatbot benchmarking
remains an open area of research, after careful analysis and applying
different methods to evaluate performance we firmly believe that the use of
the E2E Benchmark as we proposed, is the best option.

## References

  * (1) I. W. Advertising, “The ultimate guide to machine-learning chatbots and conversational ai: Ibm watson advertising.”
  * (2) P. B. Brandtzaeg and A. Følstad, “Why people use chatbots,” in Internet Science: 4th International Conference, INSCI 2017, Thessaloniki, Greece, November 22-24, 2017, Proceedings 4, pp. 377–392, Springer, 2017.
  * (3) A. Tamkin, M. Brundage, J. Clark, and D. Ganguli, “Understanding the capabilities, limitations, and societal impact of large language models,” arXiv preprint arXiv:2102.02503, 2021.
  * (4) Z. Wang, “Modern question answering datasets and benchmarks: A survey,” arXiv preprint arXiv:2206.15030, 2022.
  * (5) D. Cer, Y. Yang, S.-y. Kong, N. Hua, N. Limtiaco, R. S. John, N. Constant, M. Guajardo-Cespedes, S. Yuan, C. Tar, et al., “Universal sentence encoder,” arXiv preprint arXiv:1803.11175, 2018.
  * (6) N. Reimers and I. Gurevych, “Sentence-bert: Sentence embeddings using siamese bert-networks,” arXiv preprint arXiv:1908.10084, 2019.


# Accelerating the computation of quantum brachistochrone

Ding Wang1 Haowei Shi1 Yueheng Lan1,2<EMAIL_ADDRESS>1School of Science,
Beijing University of Posts and Telecommunications, Beijing 100876, China.
2State Key Lab of Information Photonics and Optical Communications,Beijing
University of Posts and Telecommunications, Beijing 100876, China

###### Abstract

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 redundant complexity involved in previous
formulations. 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.

Quantum control; Quantum brachistchrone; Nonlinear dynamics; Spin qubits;
Entanglement

###### pacs:

03.67.Lx, 02.30.Xx, 02.30.Yy, 03.65.-w

††preprint: APS/123-QED

##  I Introduction

Ever since people realized that quantum computers are much more powerful than
their classical counterpartsFeynman (1982); SHOR (1994), researchers have been
working hard on the problem for years, which has been making considerable
progress and leads to a dramatic increase in the complexity of controllable
quantum systems, from prototyped small quantum devices to machines with
thousands of qubits (a typical example is the D-wave quantum annealing system
Boixo  _et al._ (2014)). 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 quantum control, which is trying to
minimize 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 computing architectures. Therefore, a
lot of research about the quantum speed limit (QSL) has been reported, which
gives a lower bound of the time needed by a quantum system to travel from the
initial to the final state and is estimated with the average energy or its
variance del Campo  _et al._ (2013); Deffner and Lutz (2013); Taddei  _et al._
(2013); Mirkin  _et al._ (2016). In addition, time optimal control provides a
physical frameworkSchulte-Herbruggen  _et al._ (2005) for defining complexity
of quantum algorithms, probably superior to traditional concepts Nielsen and
Chuang (2011). For example Nielsen et al. proposed a new quantum computing
criterion by using the Riemannian geometry, namely transforming the problem of
finding an optimal quantum trajectory into a geometric one in the space of
HamiltoniansNielsen  _et al._ (2006a), and revealed that the complexity of a
quantum gate is related to the problem of optimal control Nielsen  _et al._
(2006b).

In a series of work about QSL del Campo  _et al._ (2013); Deffner and Lutz
(2013); Taddei  _et al._ (2013); Mirkin  _et al._ (2016), 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
Khaneja  _et al._ (2005); Maximov  _et al._ (2008), 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 shows
high flexibility. However, the convergence is slowing down when the difference
is close to some extreme values. Also, improper selection of step size often
leads to a tortuous search route. In order to overcome this difficulty, the
Krotov method Maday and Turinici (2003); Maximov  _et al._ (2008); Zhu and
Rabitz (1998); Shi and Rabitz (1990) introduces a second term to penalize the
undesired effects. A key feature of the algorithm is that all information is
utilized at each time step which often consistently reduces the objective
functional and accelerates the iteration. However, when parameters are too
small, its convergence rate deteriorates, and sometimes it becomes even worse
that the algorithm turns unstable. Another way to find some time-optimal
solution and its precise characterization in case of different constraints on
the Hamiltonian is through invoking interesting theories such as Pontryagin
maximum principle and the geometry of the unitary group Khaneja  _et al._
(2001); Boozer (2012); Hegerfeldt (2013). A significant piece of work was
reported by Carlini Carlini  _et al._ (2006), who derived a quantum
brachistochrone equation (QBE) with variational method which is widely used
and plays a crucial role Meiss (1992); Bunimovich (1995); Lan and Cvitanovic
(2004); Ghoussoub and Moameni (2007); Dong and Lan (2014); Wang  _et al._
(2018) in classical systems. The scheme produces a Schrodinger equation
together with the optimal control strategy under an energy constraint. Except
for some special cases Carlini  _et al._ (2007, 2008, 2014), analytic
solutions are hard to obtain and hence numerical solution is the only resort.
A scheme based on a geometric point of view is designed recently Wang  _et
al._ (2015, 2017), which depicts the quantum brachistochrone paths as
geodesics on the constraining manifold and solves the problem by solving a
family of geodesic equations. However, there is an infinite number of geodesic
families present even for locally time-optimal solutions, which makes the
scheme less convenient to use. Another drawback is the complexity of the
commutators introduced by different terms in the Hamiltonian to the Euler-
Lagrange equation (see Eq.(2) in Ref.Wang  _et al._ (2015)). As a result, the
computation load for multiple qubits could get very high and hardly be applied
in practice. The problem of how to reach an effective dimension reduction and
achieve fast maneuver of quantum bits is worth further exploration.

In this paper, we uncover the source of computational complexity from the
perspective of a variational principle, and in the process propose an
alternative set of equations with much improved efficiency, especially for
multiple bit with limited control parameters. We first slightly extend the
formulation in Ref. Carlini  _et al._ (2006) and derive all the necessary
equations for later convenience. By introducing a set of new variables, the
redundancy originated from the commutation relation between operators gets
eliminated, which enables a very efficient computation of optimal paths in the
presence of multiple qubits. A relaxation scheme is the designed and
implemented to solve the problem involving quantum entanglement when
interaction between different bits comes into paly. With several continuous
symmetries identified in the Lagrangian formulation, it is possible to remove
the annoying neutral directions in the numerical computation. In all the
examples we tried, the simplest optimal paths seem to always bear the time-
reversal symmetry if the initial and target state are related by a symmetry of
the Hamiltonian.

The paper is organized as follows. In Sec. II, a Lagrangian formation of the
QBE is reproduced with an introduction of uncontrollable terms (the
coefficients of which are fixed values). All the necessary variational
equations are derived and written in convenient vector forms. The solution of
these equations is checked analytically and verified numerically in Sec. III,
in the absence of interaction. To deal with multiple qubits with interaction,
a new set of equations are written down in Sec. IV by change of variables from
the original set, which achieves a great reduction of complexity. The
evolution of the entanglement on an optimal path is always unimodal, which is
demonstrated in all the examples in Sec. V. The results are summarized in the
final section.

##  II the variational principle

Without loss of generality, the Hamiltonian of a quantum system that is of
interest could be written as:

$H=\sum_{j=1}^{m}\xi_{j}(t)A_{j}\,,$ (1)

where $\xi_{j}(t)$’s are real numbers and $A_{j}$’s are Hermitian operators
that satisfy $TrA_{j}=0$ and $Tr(A_{j}A_{k})=\delta_{jk}~{}(j,k\leq m)$. The
traceless Hamiltonian could be divided into two parts: the controllable part
$\bm{\mathcal{A}}=\sum\limits_{j=1}^{n}\xi_{j}(t)A_{j}$ with each
$\xi_{j}(j=1,...,n)$ subject to external control and the uncontrollable part
$\bm{\mathcal{B}}=\sum\limits_{j=n+1}^{m}\xi_{j}(t)A_{j}$ with $\xi_{j}$’s
being constant or changing in a pre-determined way, often used to model
interaction between qubits. For later convenience we use $\bm{u}$ and $\bm{v}$
to label coefficients of the controllable part
$\bm{u}=(\xi_{1}\,,\xi_{2}\,,\cdots\,,\xi_{n}\,,0\,,\cdots\,,0)^{t}$ and the
uncontrollable part
$\bm{v}=(0\,,0\,,\cdots\,,\xi_{n+1}\,,\xi_{n+2}\,,\cdots\,,\xi_{m})^{t}$, and
hence $\bm{\xi}=\bm{u}+\bm{v}$ . In this work, we build a formulation
applicable to the general case which includes both components.

###  II.1 A variational scheme for the QBE

The action integral is written as Carlini  _et al._ (2006):

$S(\psi\,,H\,,\phi\,,\lambda,,\lambda^{\prime})=\displaystyle{\int}dt\dfrac{\sqrt{\langle\dot{\psi}|(1-P)|\dot{\psi}\rangle}}{\Delta
E}+\big{(}i\langle\dot{\phi}|\psi\rangle+\langle\phi|H|\psi\rangle+c.c.\big{)}+\lambda(\sum\limits_{j=1}^{m}\frac{\xi_{j}^{2}}{2}-\omega^{2})+\sum\limits_{j=n+1}^{m}\lambda^{{}^{\prime}}_{j}(\xi_{j}-Q_{j})\
\,,$ (2)

where $P(t)=|\psi(t)\rangle\langle\psi(t)|$ is the projection to the state
$|\psi(t)\rangle$ which is a normalized wavefunction defined in the
configuration space, $\phi$ is an auxiliary wavefunction and $c.c.$ denotes
complex conjugates. $(\Delta
E)^{2}\equiv\langle\psi|H^{2}|\psi\rangle-\langle\psi|H|\psi\rangle^{2}$ is
the energy variance. Parameters $\lambda$ and $\lambda_{j}^{{}^{\prime}}$ are
Lagrange multipliers and all $\\{\lambda^{\prime}_{j}\\}_{j=1,2,...,n}$ are
set to zero, $\omega$ is a given constant, the $Q_{j}$’s denote the coupling
constants in the Hamiltonian and the Planck’s constant $\hbar$ is chosen to be
1 for simplicity.

Now we take the variation of the action Eq. (2) to obtain the equation of
motion and all the constraints imposed to the variables.

####

The variation with respect to $\phi$ gives

$i|\dot{\psi}\rangle=H|\psi\rangle\,,$ (3)

which is the usual Schrödinger equation. In the projective space
$\mathbb{C}P^{n-1}$, the Fubini-Study line element is written as
$ds^{2}=\langle d\psi|(1-P)|d\psi\rangle$, which measures the displacement in
the angular direction. Here, it is easy to check that $ds^{2}=(\Delta
E)^{2}dt^{2}$. Hence, the first part of Eq. (2) refers to the total time of
the process and the other parts embody differential or other constraints.

####

Here, we take the variation with respect to $\lambda$ and
$\lambda^{\prime}_{j}$ separately, resulting in
$\frac{1}{2}\sum\limits_{j=1}^{m}\xi_{j}^{2}=\omega^{2}$, which defines the
size of the Hamiltonian as a finite energy constraint, and $\xi_{j}=Q_{j}$
($n<j\leq m$), which depicts the influence of generic interactions.

####

The variation with respect to $\psi$ gives

$i\dfrac{d}{dt}\bigg{[}\dfrac{H-\langle H\rangle}{2(\Delta
E)^{2}}\bigg{]}|\psi\rangle-i|\dot{\phi}\rangle+H|\phi\rangle=0\,,$ (4)

where $\langle\cdot\rangle$ refers to the average with respect to
$|\psi\rangle$. Eqs. $(\ref{Eq1})$ and $(\ref{Eq4})$ are seen in the
literatureCarlini  _et al._ (2006) but the equations below appear new.

####

The variation with respect to $H$ could now be taken as multivariate change
with respect to $\\{\xi_{j}\\}_{j=1,2,...m}$, resulting in

$D_{j}=\frac{1}{2\Delta
E^{2}}\sum\limits_{k}F_{jk}\xi_{k}-\lambda\xi_{j}-\lambda^{{}^{\prime}}_{j}~{}~{}~{}~{}(\text{with}~{}\lambda^{{}^{\prime}}_{j}=0~{}\text{for}~{}j\leq
n),$ (5)

where $F_{jk}=\langle A_{j}A_{k}+A_{k}A_{j}\rangle-2\langle
A_{j}\rangle\langle A_{k}\rangle$, $D_{j}=\langle\phi|A_{j}|\psi\rangle+c.c.$
and $\sum\limits_{k}$ stands for $\sum_{k=1}^{m}$ which is also the case in
all the following equations. The detailed derivation of Eq (5) has been
relegated to Appendix A. Next, we derive the equations of motion for $\lambda$
and $\\{\xi_{j}\\}_{j=1,...,n}$.

A multiplication of both sides of Eq. (5) with $\\{\xi_{j}\\}_{j=1,...,n}$,
the controllable components, and then a summation over j results in

$\lambda=\dfrac{1-\bm{D}\cdot\bm{\xi}-\sum\limits_{j=n+1}^{m}\bigg{(}\sum\limits_{k}\dfrac{\xi_{j}F_{jk}\xi_{k}}{2\Delta
E^{2}}-D_{j}\xi_{j}\bigg{)}}{\omega^{\prime 2}}\,,$ (6)

where $\bm{D}\cdot\bm{\xi}=\sum\limits_{j=1}^{m}D_{j}\xi_{j}$ and
$\omega^{\prime 2}=\sum\limits_{j=1}^{n}\xi_{j}^{2}$ is constant (the
constraint on the controllable part of the Hamiltonian). Noticing that
$\\{\dot{\xi}_{j}\\}_{j=n+1,...m}\equiv 0$, we could obtain:

$\dot{\lambda}=\dfrac{\sum\limits_{j=n+1}^{m}\bigg{\\{}-\dfrac{\xi_{j}}{2\Delta
E^{2}}\bigg{[}\sum\limits_{k}(\dot{F_{jk}}\xi_{k}+F_{jk}\dot{\xi_{k}})-\dfrac{\tilde{\bm{F}}\cdot\dot{\bm{\xi}}}{\Delta
E^{2}}\tilde{F_{j}}\bigg{]}+\dot{D_{j}}\xi_{j}\bigg{\\}}}{\omega^{\prime
2}}\,,$ (7)

where $\tilde{\bm{F}}=\bm{\xi}\cdot F$ and
$\tilde{F_{j}}=\sum_{k}F_{jk}\xi_{k}$. Here we used the fact that
$\dfrac{d}{dt}\bm{D}\cdot\bm{\xi}=0$ (details in Appendix C),
$\bm{\xi}\cdot\bm{C}=\bm{C}\cdot\bm{\xi}=0$ (details in Appendix E) and

$\dfrac{d}{dt}D_{j}=\dfrac{1}{2\Delta
E^{2}}(\sum\limits_{k}F_{jk}\dot{\xi}_{k}-\dfrac{\tilde{\bm{F}}\cdot\dot{\bm{\xi}}}{\Delta
E^{2}}\tilde{F_{j}})+\sum\limits_{k}C_{jk}D_{k}\,,$ (8)

where $C_{jk}$ satisfies $[H,A_{j}]=-i\sum_{k}C_{jk}A_{k}$, with $C_{jk}$
being a real number. For a detailed derivation of Eq. (8), please check
Appendix B. Feeding Eq. (8) into Eq. (7) results in

$\displaystyle\dot{\lambda}=$
$\displaystyle\dfrac{\sum\limits_{j=n+1}^{m}\sum\limits_{k=1}^{m}\xi_{j}C_{jk}(D_{k}-\frac{\tilde{F}_{k}}{2\Delta
E^{2}})}{\omega^{\prime 2}}$ $\displaystyle=$ $\displaystyle-\dfrac{\bm{v\cdot
C\cdot\lambda}^{\prime}}{\omega^{\prime 2}}\,,$ (9)

with $\dot{\bm{F}}=\bm{C}\cdot\bm{F}+\bm{F}\cdot\bm{C}^{T}$ (Appendix D).
Taking the time derivative of both sides of Eq. (5), and together with Eq.
(8), we get

$\dot{\bm{u}}=-\dfrac{\bm{C}\cdot\bm{D}}{\lambda}+\dfrac{\bm{C}\cdot
F\cdot\bm{\xi}}{2\lambda\Delta E^{2}}+\dfrac{(\bm{v\cdot
C\cdot\lambda}^{\prime})}{\lambda\omega^{\prime 2}}\bm{u}\,.$ (10)

With Eq. (8), Eq. (II.1), Eq. (10) and the Schrödinger equation Eq. (3), the
optimal path can be calculated now.

Here let’s discuss a special case with $H\in\mathcal{A}$, namely
$\bm{\xi}=\bm{u}$, which means that all parts are controllable and thus
$\bm{v}=0$. It is obvious that this case is simple since

$\displaystyle D_{j}$ $\displaystyle=\frac{1}{2\Delta
E^{2}}\sum\limits_{k}F_{jk}u_{k}-\lambda u_{j},$ (11)
$\displaystyle\dfrac{d}{dt}D_{j}$ $\displaystyle=\dfrac{1}{2\Delta
E^{2}}(\sum\limits_{k}F_{jk}\dot{u}_{k}-\dfrac{\tilde{\bm{F}}\cdot\dot{\bm{u}}}{\Delta
E^{2}}\tilde{F_{j}})+\sum\limits_{k}C_{jk}D_{k}.$ (12)

Hence Eq.(II.1) implies that $\dot{\lambda}=0$. Eq. (10) could be reduced to

$\dot{\bm{u}}=-\dfrac{\bm{C}\cdot\bm{D}}{\lambda}+\dfrac{\bm{C}\cdot\bm{F}\cdot\bm{u}}{2\lambda\Delta
E^{2}}.$ (13)

Considering $\bm{u}\cdot\bm{C}=0$ (Appendix E) and feeding Eq. (11) into Eq.
(13) leads to a trivial solution $\dot{\bm{u}}=0$, which indicates that
$\bm{u}$ is a constant vector and the evolution corresponds to a simple
rotation with a fixed axis on the Bloch sphere.

###  II.2 Gauge symmetries in the Lagrangian

In this section, we discuss the gauge symmetries in the Lagrangian Eq. (2),
which create neutral directions in the evolution and may bring complication in
numerical calculation. First, consider a transform
$\hat{U}_{1}:\phi\to\phi+ia\psi$, $a\in\mathcal{R}$, with which the Lagrangian
becomes

$S^{\prime}=S+\int dt~{}(a\langle\dot{\psi}|\psi\rangle-
ia\langle\psi|H|\psi\rangle+c.c.)\,,$ (14)

The expression $\langle\dot{\psi}|\psi\rangle+c.c$ is a total differentiation,
which could be integrated out. With the Hermiticity of $H$, the extra integral
in Eq. (14) vanishes and hence $\hat{U}_{1}$ is a symmetry transform. Next, we
consider another transform
$\hat{U}_{2}:\phi\rightarrow\dfrac{\lambda+b}{\lambda}\phi,~{}\lambda\rightarrow\lambda+b,~{}\bm{\lambda}^{\prime}\rightarrow\dfrac{\lambda+b}{\lambda}\bm{\lambda}^{\prime},~{}b\in\mathcal{R}$
which sends Eq. (2) to

$\displaystyle S(\psi,H,\phi,\lambda,\lambda^{\prime})$ $\displaystyle=\int
dt\dfrac{\sqrt{\langle\dot{\psi}|(1-P)|\dot{\psi}\rangle}}{\Delta E}$
$\displaystyle+\dfrac{\lambda+b}{\lambda}\left[(i\langle\dot{\phi}|\psi\rangle+\langle\phi|H|\psi\rangle+c.c.)+\lambda(\sum_{j=1}^{m}\frac{\xi_{j}^{2}}{2}-\omega^{2})+\sum_{j=n+1}^{m}\lambda^{\prime}_{j}(\xi_{j}-Q_{j})\right]\,,$
(15)

which only changes the values of the invisible Lagrange multipliers
$\phi,~{}\lambda,~{}\bm{\lambda}^{\prime}$ while the equation for the
observables $\psi,~{}\bm{u}$ remain intact. Actually, it is most obvious in
the simplified formulation Eq. (18) below, where
$\bm{\Omega},~{}\bm{u},~{}\dot{\bm{\Omega}},~{}\dot{\bm{u}}$ keep invariant
under the transform $\hat{U}_{2}$. Hence $\hat{U}_{2}$ is a symmetry
transform.

##  III The evolution of Hamiltonian as a rotation

The equation $\bm{\xi}\cdot\bm{C}=0$ tell us that
$\bm{v}\cdot\bm{C}=-\bm{u}\cdot\bm{C}$. Substituting the Eq. (5) into Eq. (10)
results in

$\dot{\bm{u}}=(I-P_{u})\cdot\dfrac{\bm{C}\cdot\bm{\lambda}^{\prime}}{\lambda}\,,$
(16)

where $I$ is the identity matrix and $P_{u}$ is the projection operator
$\dfrac{|\bm{u}\rangle\langle\bm{u}|}{|\bm{u}|^{2}}$ ($\omega^{\prime
2}=\sum\limits_{j=1}^{n}\xi_{j}^{2}=\sum\limits_{j=1}^{n}u_{j}^{2}=|\bm{u}|^{2}$).
Here we use the Dirac notation to emphasize the projection, whereas $\bm{u}$
is actually real. With the fact that $\bm{\dot{v}}=0$, it is easy to check
that $\bm{\xi}\cdot\bm{\dot{\xi}}=\bm{u}\cdot\dot{\bm{u}}=0$. And the
uncontrollable part of the time derivative of Eq. (5) similarly gives the
evolution of $\bm{\lambda}^{\prime}$

$\dot{\bm{\lambda}}^{\prime}=\bm{C}\cdot\bm{\lambda}^{\prime}+\dfrac{\bm{v}\cdot\bm{C}\cdot\bm{\lambda}^{\prime}}{\omega^{\prime
2}}\bm{v}=(I+\dfrac{|\bm{v}|^{2}}{|\bm{u}|^{2}}P_{\bm{v}})\cdot\bm{C}\cdot\bm{\lambda}^{\prime}.$
(17)

where $P_{v}=\dfrac{|\bm{v}\rangle\langle\bm{v}|}{|\bm{v}|^{2}}$ and
$\dfrac{|\bm{v}|^{2}}{|\bm{u}|^{2}}$ is a real constant. After, defining
$\bm{\Omega}=-\dfrac{2}{\sqrt{TrI}}\dfrac{\bm{\lambda^{\prime}}}{\lambda}$,
the angular velocity of the rotation of $\bm{u}$ can be described with
$\bm{\Omega}$. Furthermore, the QBE problem can be characterized with a
simpler set of equations

$\displaystyle\dot{\bm{u}}$
$\displaystyle=-\dfrac{\sqrt{TrI}}{2}(I-P_{u})\cdot\bm{C}\cdot\bm{\Omega}$
(18) $\displaystyle\dot{\bm{\Omega}}$
$\displaystyle=(I+\dfrac{|\bm{v}|^{2}}{|\bm{u}|^{2}}P_{\bm{v}})\cdot\bm{C}\cdot\bm{\Omega}-\dfrac{\sqrt{TrI}}{2\omega^{\prime
2}}(\bm{v}\cdot\bm{C}\cdot\bm{\Omega})\bm{\Omega}$

along with the Schrödinger equation $\dot{\psi}=-iH\psi$.

According to Appendix E, $\bm{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 Eq. (18),
the equation of a single qubit in 3-dimensional space is appreciably
simplified to
$\dot{\bm{u}}=\bm{\Omega}\times\bm{\xi}-\hat{\bm{u}}\cdot(\bm{\Omega}\times\bm{\xi})~{}\hat{\bm{u}}=\bm{\Omega}\times\bm{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, $\bm{\Omega}$ is in the direction of
$\bm{\lambda}^{\prime}$ which is along $\bm{v}$ in the case of three
dimensions. Here the cross product holds for each $SO(3)$ Bloch sphere, which
indicates that $\bm{u}$ is rotating with the angular velocity $\bm{\Omega}$.

As $\bm{\Omega},\bm{v}$ is in the same direction for a single qubit. After
multiplying $\bm{\Omega}$ with Eq. (18) we see that $|\bm{\Omega}|$, the
rotation speed of $\bm{u}$, is a constant vector. In fact, in this case
$\bm{\lambda}^{\prime}=0$ and thus $\bm{\Omega}=0$, $\dot{\bm{u}}=0$ as given
in the previous section. It is easy to see that the spin vector
$\langle\bm{\sigma}\rangle$ ($\sigma_{j=\\{x,y,z\\}}$ are the Pauli matrices)
is rotating around $\bm{u}$, and the evolution of wave function is
analytically computed with rotation matrices.

###  III.1 Analytical solution of a single qubit

The evolution path of $\psi$ can be analytically addressed, as a rotation with
the angular velocity $\bm{\omega}_{\mathrm{eff}}$ around an axis $\bm{u}$
rotating itself at $\bm{\Omega}$. The effective angular velocity
$\bm{\omega}_{\mathrm{eff}}$ of spin vector is determined with

$\bm{\omega}_{\mathrm{eff}}=2\bm{B}=\dfrac{2}{\sqrt{TrI}}\bm{u}=\sqrt{2}\bm{u}\,.$
(19)

where $B$ is an effective magnetic field and the Hamiltonian
$H=\bm{\sigma}\cdot\bm{B}=\sum\limits_{j}A_{j}u_{j}$. This is derived from the
fact that the SU(2) generator of a rotation with $\Delta\phi$ is
$e^{i\bm{\sigma}\frac{\Delta\phi}{2}}$. In the frame of reference rotating
with $\bm{\Omega}$, the motion of $\psi$ is viewed as a fixed-axis rotation
with the angular velocity
$\bm{\Omega}^{\prime}\\!=\bm{\omega}_{\mathrm{eff}}\\!-\bm{\Omega}$.

For instance, consider the numerically simulated case Carlini  _et al._ (2006)
where $\psi(0)=|+x\rangle$, $\psi_{target}=|-x\rangle$, $\bm{B}(0)$ is
$\omega$ in the $-y$ direction and $\bm{\Omega}$ lies in the $z$ direction
with $Tr(H\sigma_{z})=0$. Explicitly,
$\langle\bm{\sigma}\rangle(0)=(1,0,0),\bm{B}(0)=(0,-\omega,0),\bm{\Omega}\equiv(0,0,\Omega)$.
Now the evolution turns out to be simply a rotation around
$\bm{\Omega}^{\prime}$ followed by a reversal rotation around $\bm{\Omega}$
back to the static frame. The rotation matrix
$R(\phi,\theta,\Psi)=R_{\bm{\Omega}}^{-1}(\phi)R_{\bm{\Omega}^{\prime}}(\theta)=R_{z}^{-1}(\phi)R_{x}^{-1}(\Psi)R_{y}(\theta)R_{x}(\Psi)$
(see Appendix G ). Here $(\phi,\theta,\Psi)$ is valued $(\Omega
t,\Omega^{\prime}t,\mathrm{tan^{-1}}\dfrac{\Omega}{\omega_{\mathrm{eff}}})$,
where $\Omega^{\prime}=\sqrt{\Omega^{2}+\omega_{\mathrm{eff}}^{2}}$. Then it
is easy to check that

$\langle\bm{\sigma}\rangle(t)=R(\phi,\theta,\Psi)\cdot\langle\bm{\sigma}\rangle(0)=\begin{pmatrix}&~{}\mathrm{cos}\\!~{}\Omega
t\\!~{}\mathrm{cos}\\!~{}\Omega^{\prime}t+\dfrac{\Omega}{\Omega^{\prime}}~{}\mathrm{sin}\\!~{}\Omega
t\\!~{}\mathrm{sin}\\!~{}\Omega^{\prime}t\\\ &-\mathrm{sin}\\!~{}\Omega
t\\!~{}\mathrm{cos}\\!~{}\Omega^{\prime}t+\dfrac{\Omega}{\Omega^{\prime}}~{}\mathrm{cos}\\!~{}\Omega
t\\!~{}\mathrm{sin}\\!~{}\Omega^{\prime}t\\\
&\dfrac{\omega_{\mathrm{eff}}}{\Omega^{\prime}}\mathrm{sin}\\!~{}\Omega^{\prime}t\end{pmatrix}\,,$
(20)

which gives a family of analytic solutions if we choose $(\Omega
T,\Omega^{\prime}T,\dfrac{\Omega}{\omega_{\mathbf{eff}}})=(k\pi,l\pi,\dfrac{k}{\sqrt{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 contraints and cases
of $k\not=0$ correspond to the oscillating solutions with $l-1$ nodes, and
setting L= $l-1$ for later convenience. The oscillating curve in FIG. 2 plots
the solution with $k=1,L=1$.

###  III.2 Analytical solution for multiple qubits

First we give an example of two qubits with
$H=B_{x}^{(1)}\sigma_{x}^{(1)}+B_{y}^{(1)}\sigma_{y}^{(1)}+B_{x}^{(2)}\sigma_{x}^{(2)}+B_{y}^{(2)}\sigma_{y}^{(2)}$.
Similarly, our purpose is to flip both qubits (where
$\psi(0)^{(1)}=|+x\rangle$ and $\psi_{target}^{(1)}=|-x\rangle$ for qubit one,
$\psi(0)^{(2)}=|+y\rangle$ and $\psi_{target}^{(2)}=|-y\rangle$ for qubit
two). Obviously, these two qubits are both reversing their states, and because
$x-$ and $y-$ directions are totally equivalent there qubits should share the
same energy at the ground state(namely
$\omega_{\mathbf{eff}}^{(1)}=\omega_{\mathbf{eff}}^{(2)}$). Fig. 3 shows a
local optimal solution of the constrained two-qubits system with $L=1$.
$\bm{\Omega}^{(1)}=(0,0,-\dfrac{2}{\sqrt{3}}\omega),~{}\bm{\Omega}^{(2)}=(0,-\dfrac{2}{\sqrt{3}}\omega,0)$.
It is appearent that the evolution of $\bm{B}$ turns out to be a rotation as
expected.

Actually, the local solutions can be derived analytically as long as the
entangled terms are missing in the Hamiltonian. According to the rotation
picture, the analytical approach can be applied to non-entangled many-qubit
systems with an arbitrarily large number of qubits. Here in Fig. 3 the
solution corresponds to $(k^{(1)},L^{(1)})=(k^{(2)},L^{(2)})=(1,1)$ as the
oscillating solution in Fig. 2, so they naturally share the same geometric
property. Later on, these solutions would be good initial guesses in the case
of weak entanglement (see IV.2).

##  IV An Accelerated Approach to the QBE

In practice, there are usually a large number of uncontrollable terms in the
Hamiltonian. Sometimes for the convenience of computation and design, almost
all the interaction terms are made to vanish. At this moment the dimension of
control parameters is appreciably low, however, 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 $\phi$ previously replaced by $\bm{D}$.

In general situations we see that the tremendous complexity is mainly brought
by the factor $\bm{C}\cdot\bm{D}$, since the commutation relation (52) will
incur a large number of or infinitely many more $A_{k}$’s and thus $D_{k}$’s
in the QBE equations when there are multiple interacting qubits. It is
frustrating that most of the previous researches Carlini  _et al._ (2006);
Wang  _et al._ (2015) are confronted with similar problems but no special
attention has been paid to this problem. Here, we introduce an effective
deduction scheme based on the fact that $\bm{C}\cdot\bm{D}$ is immediately
derivable from $\phi$. Technically, consider

$\phi^{\prime}=\phi-\dfrac{1}{2\Delta E^{2}}(H-\langle H\rangle)\psi,$ (21)

then Eq. (5) turns to

$D^{\prime}_{j}=\langle\phi^{\prime}|A_{j}|\psi\rangle+c.c.=-\lambda\xi_{j}-\lambda^{\prime}_{j},~{}~{}~{}~{}(\text{with}~{}\lambda^{\prime}_{j}=0~{}\text{for}~{}j\leq
n).$ (22)

Thus ${\bm{C}\cdot
D^{\prime}_{j}}={\bm{C}\cdot\langle\phi^{\prime}|A_{j}|\psi\rangle}+c.c.$.
Furthermore, it is much easier to evolve $\phi^{\prime}$ instead of $\phi$
since Eq. (4) yields

$i\dot{\phi}-H\phi=\frac{i}{2\Delta
E^{2}}[(\dot{H}-\langle\dot{H}\rangle)\psi-\dfrac{\tilde{\bm{F}}\cdot\dot{\bm{\xi}}}{\Delta
E^{2}}(H-\langle H\rangle)\psi]\\\ \,,$ (23)

which gives

$\displaystyle i\dot{\phi}$ $\displaystyle{}^{\prime}=H\phi^{\prime}\,.$ (24)

In analogy to previous derivations in Sec. III, it is easy to check that Eq.
(18) now holds in an alternative form

$\displaystyle\dot{\bm{u}}$
$\displaystyle=(I-P_{u})\cdot\dfrac{\bm{C}\cdot\bm{\lambda}^{\prime}}{\lambda},$
(25) $\displaystyle i\dot{\phi}$ $\displaystyle{}^{\prime}=H\phi^{\prime},$
$\displaystyle\dot{\lambda}$ $\displaystyle=-\dfrac{\bm{v\cdot
C\cdot\lambda}^{\prime}}{2\omega^{2}},$ $\displaystyle\dot{\psi}$
$\displaystyle=-iH\psi\,,$

where from Eq. (22),
$\bm{C}\cdot\bm{\lambda}^{\prime}=-\bm{C}\cdot\bm{D}^{\prime}$ with
$\sum\limits_{j}C_{ij}D^{\prime}_{j}={\langle\phi^{\prime}|[H,A_{i}]|\psi\rangle}+c.c.$,
which could be evaluated directly with the computed states $\phi^{\prime}$ and
$\psi$ and the known $H$ and $A_{i}$’s. In conclusion, the whole QBE problem
is fully described by Eq. (25). Terms like
$\bm{C}\cdot\bm{F}\cdot\bm{\xi}/2\Delta E^{2}$ in $\bm{C}\cdot\bm{D}$ are
eliminated now because of the reintroduction of the state variable
$\phi^{\prime}$, and hence the previous large number of equations needed to
evaluate $\bm{C}\cdot\bm{D}$ are simply replaced by the evolution of
$\phi^{\prime}$. Since the number of degrees of freedom (DOF) of vector
$\phi’$ is much smaller than that of matrix $\bm{C}$, the computation will be
significantly accelerated.

Based on Eq. (5), (22) and the fact that $TrH=0$ (i.e. $\xi_{0}=0$), we
immediately obtain

$\displaystyle 2Re(\langle\phi|\psi\rangle)$
$\displaystyle=D_{0}=\sum_{j}\dfrac{F_{0j}}{2\Delta
E^{2}}\xi_{j}-\lambda\xi_{0}\equiv 0\,,$ (26) $\displaystyle
Re(\langle\phi^{\prime}|\psi\rangle)$
$\displaystyle=D^{\prime}_{0}=-\lambda\xi_{0}\equiv 0\,,$

where $F_{0j}=0$ and $A_{0}=\sigma_{0}$, which is the unit matrix. Hence the
component of $\phi$ or $\phi^{\prime}$, being parallel to $\psi$, vanishes.
Notice that $\bm{u}$, in the absence of uncontrollable components, is parallel
to $Re(\langle\phi^{\prime}|\bm{\sigma}|\psi\rangle)$ according to Eq. (22).
So $\bm{u}\cdot\langle\bm{\sigma}\rangle\propto
Re(\langle\phi^{\prime}|\bm{\sigma}|\psi\rangle)\cdot\langle\bm{\sigma}\rangle=0$
(proved in Appendix F) for every qubit.

Now a two-point boundary value problem with Eq.(25) is left for us to solve.
To begin with, we count the number of independent variables. The number of
different variables in Eq. (25) is
$N=N_{u}+N_{\phi^{\prime}}+N_{\lambda}+N_{\psi}+N_{\lambda^{\prime}}$, while
the number of the constraints is
$N^{\prime}=N_{u}+N_{\lambda^{\prime}}+2N_{\psi}+1$ given by Eq.(22) (giving
$N_{u}+N_{\lambda^{\prime}}$), the energy constraint
$\sum\limits_{j=1}^{n}\xi_{j}^{2}=\omega^{\prime 2}$, the prescrible initial
and final states $\psi(0)=\psi_{0}$ and $\psi(T)=\psi_{target}$ (giving
$2N_{\psi}$), then the number of independent variables is $N-N^{\prime}=0$
(where $N_{\lambda}=1$ and $N_{\phi^{\prime}}=N_{\psi}$).

In the current paper the _shooting method_ (IV.1) and the idea of _relaxation_
play a major role in finding the optimal path when Hamiltonian $H$ contains
interaction terms. We give some details below.

###  IV.1 Shooting Method

The shooting method aims to implement the multidimensional Newton-Raphson
method to locate zeros of $N_{target}$ functions with $N_{target}$ varibales,
which are obtained by integrating $N$ differential equations in our examples.

For a general dynamical system, the ordinary differential equation (ODEs) is
written as

$\dot{\bm{x}}=\bm{v}(\bm{x}(t)),$ (27)

where $\bm{x}\in\mathbb{R}^{d}$, $t\in\mathbb{R}$, and the state of the system
obtained by integrating Eq.(27) at time $t$ is $\bm{x}(t)=f^{t}(\bm{x})$. The
corresponding Jacobian matrix is
$J(\bm{x},t)=\frac{\partial\bm{x}(t)}{\partial\bm{x}(0)}$, which could be
obtained by integrating

$\frac{dJ}{dt}=AJ,~{}~{}~{}A_{ij}=\frac{\partial{\bm{v}_{i}}}{\partial{\bm{x}_{j}}},~{}~{}~{}\text{with
}J(\bm{x},0)=\bm{1}.$ (28)

Thus, the initial displacement $\delta\bm{x}(0)$ becomes $\delta\bm{x}(T)$ at
time $T$, namely

$\delta\bm{x}(T)=J(\bm{x},T)\cdot\delta\bm{x}(0)\,.$ (29)

Define a discrepancy vector $\bm{X}=\bm{x}(T)-\bm{x}_{target}$. To procure a
vanishing discrepancy, the correction $\Delta x_{0}$ should satisfy

$J\cdot\Delta\bm{x}_{0}=-\bm{X},$

to reach a better approximation

$\bm{x}_{0,new}=\bm{x}_{0,old}+\Delta\bm{x}_{0}\,.$ (30)

Typically, the Jacobian matrix $J$ is numerically evaluated by

$J_{ij}=\dfrac{\bm{X}_{i}(\bm{x}_{0,1},\ldots,\bm{x}_{0,j}+\Delta\bm{x}_{0,j},\ldots,T)-\bm{X}_{i}(\bm{x}_{0,1},\ldots,\bm{x}_{0,j},\ldots,T)}{\Delta\bm{x}_{0,j}}\,.$
(31)

In the current case, the adjustable part is the free parameter $\phi$ and the
discrepancy is $\psi(T)-\psi_{target}$.

Let’s slow down here and have a discussion on possible additional constraints
and the treatment of the symmetries mentioned earlier. Firstly, based on
Eq.(26), $Re(\langle\delta\phi^{\prime}|\psi\rangle)=0$, which with the help
of the symmetry transformation $\hat{U}_{1}$ could always be satisfied by
taking the difference $\phi^{\prime}_{0,new}-ia\psi$, where
$ia\psi=|\psi\rangle\langle\psi|\phi^{\prime}_{0,new}\rangle$. The unitary
evolution of $\psi$ will produce a singularity for the Jacobian matrix $J$. To
cope with this problem, our strategy is to totally ignore this eigen-direction
since it corresponds to irrelevant symmetry directions. The constraint
$\sum_{j=1}^{n}\xi_{j}^{2}=\omega^{\prime 2}$ could be garanteed by the
symmetry $\hat{U}_{2}$ since $\lambda_{0}$ and thus $b$ could be obtained by
substituting $\phi^{\prime}_{0,new}$ into Eq. (22) with $j\leq n$ when
starting a new iteration.

It is likely that the discrepancy may significantly fluctuate in the first few
iterations. So it is important to properly choose the target hyperplane for
the shooting on which the target state lies. The final state $\psi(T)$ is
adjusted to ’hit’ the target state $\psi_{target}$ by adjusting the initial
position Eq.(30) and the flight time $T$.

###  IV.2 Multiple Qubits: ‘Relaxation’

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, it is possible to solve the problem progressively.
First, we invite part of the interactions back and use the non-interacting
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.

Figure 1:  Evolution of Qubit 1 on a Bloch sphere. The black line represents
the state $\psi(t)$ that evolves from the initial state (purple point on the
left) to the final state (blue points on the right) on the Bloch sphere. The
green grid inside the sphere marks the equatorial plane, and the red part is
the area swept by the magnetic field vector. Similar representations are used
in later diagrams.

Figure 2:  The evolution of $\psi(t)$ of a single qubit with
$\omega_{eff}=2|\bm{B}|=2$ and L=1 (left), 2 (right) (namely, the black path
intersects the equatorial plane L times apart from the initial and final
states).

##  V examples

###  V.1 One-Qubit Bloch Sphere

So far most of the practical quantum control systems are based on spin
magnetic moment. 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\rangle$ to $|-x\rangle$.

Assuming $\hbar=1$, the Hamiltonian in the isotropic condition is

$H=-\bm{\mu}_{S}\cdot\bm{B}=-\dfrac{\gamma}{2}\bm{\sigma}\cdot\bm{B}=-\dfrac{\gamma}{2}(B_{x}\sigma_{x}+B_{y}\sigma_{y}+B_{z}\sigma_{z})\,,$

where $\gamma=g\dfrac{q}{2m}$ is the gyromagnetic ratio. For electrons
$\gamma\approx-2$ then simply $H=\bm{B}\cdot\bm{\sigma}$. From previous
discussions, we know that $\dot{\bm{u}}=0$, i.e. $\dot{\bm{B}}=0$. So the
optimal path is naturally a geodesic curve (see Fig. 2, which shows one of the
many geodesic solutions). The simulated paths are displaced as well in Fig. 2.
The red line on the equatorial plane (green grid inside the sphere) indicates
the evolution of the external magnetic field $\bm{B}$, which does not change
throughout the evolution in the case.

In principle, the constraints $Tr(H\sigma_{j})_{j=x,y,z}=0$ allow oscillating
solutions with $\dot{\bm{B}}\not=0$. Let’s take the case of
$H=B_{x}\sigma_{x}+B_{y}\sigma_{y},~{}B_{z}\equiv 0$ as an example. In
numerical computation, we use the shooting method to identify the initial
value $\phi^{\prime}_{0}$ for the prescribed boundary condition
$\psi(0)=|+x\rangle,~{}\psi(T)=|-x\rangle$. Fig. 2 shows the resulting optimal
path which oscillates with one node on the equator, where
$\bm{\Omega}=(0,0,-\dfrac{2}{\sqrt{3}})$. The evolution of $\bm{B}$ (the area
swept by the red line in Fig. 2) corroborates the theory of rotation in Sec.
III. Of course, it is not hard to get a longer locally optimal path with more
than one node, which we will not show here.

###  V.2 Two-Qubit Bloch Sphere

Figure 3:  An optimal path with $L=1$, $\omega^{(1)}_{eff}=2$,
$\omega^{(2)}_{eff}=2$ connecting the initial state
$|+x\rangle^{(1)},|+y\rangle^{(2)}$ to the final state
$|-x\rangle^{(1)},|-y\rangle^{(2)}$ of Qubits 1 and 2 plotted on two bloch
spheres.

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 Secion V.1 above. Here, we focus on the
oscillating solutions under the anisotropic constraint featuring $B_{z}=0$. We
give an example with

$H=B_{x}^{(1)}\sigma_{x}^{(1)}+B_{y}^{(1)}\sigma_{y}^{(1)}+B_{x}^{(2)}\sigma_{x}^{(2)}+B_{y}^{(2)}\sigma_{y}^{(2)}.$
(32)

Fig. 3 depicts the local optimal path obtained by our method, which oscillates
with $L=1$. Also, we can obtain longer paths, not shown here.

Figure 4:  The time evolution of the entanglement Ep of $\psi(t)$ for
different $J_{yy}$ and L=1(a), 3(b), 2(c), 4(d). t/T represents the normalized
time where T is the total time. The initial state are
$|+x\rangle^{(1)},|+y\rangle^{(2)}$ and the final state are
$|-x\rangle^{(1)},|-y\rangle^{(2)}$ for Qubits 1 and 2 Figure 5:  The
evolution of the entanglement $E_{p}(t)$ of wavefunction $\psi(t)$ from the
initial state $\psi(0)=(0.5,0.5i,0.5,0.5i)^{T}$ to the final
state(-0.216558-0.450669i, 0.528077+0.256587i, 0.070204+0.391403i,
-0.294949+0.400223i) with J=0.65, 0.9, 1 and L=1.

Figure 6:  The evolution and the convergence. The evolution of the observed
measurement $\langle\boldsymbol{\sigma}(t)\rangle$ (a) and the control
magnetic field $\bm{B}$(t) (b), where
$\langle\sigma(t)_{m}^{(l)}\rangle=\langle\psi(t)|\sigma(t)_{m}^{(l)}|\psi(t)\rangle~{}(m=x,y,z;l=1,2)$
represents the observed value of the spin in the m-direction on the $l~{}th$
bit. The convergence of our method (c) and the GRAPE method (d), where Error
$d_{1}=|\psi(T)-\psi_{target}|$ and
$d_{2}=|1-\langle\psi_{target}|\psi(T)\rangle|$, N is the number of iteration
steps and J=1, L=1. Figure 7:  The Fubini-Study distance F(t) between
$\psi(t)$ and $\psi_{target}$ with
$\omega=2.8284,~{}2.7803,~{}2.6926,~{}2.5554,~{}2.3601$. The distance between
neighboring states along a trajectory is defined by using the Fubini-Study
line element $ds=\sqrt{\langle d\psi|(1-P)|d\psi\rangle}$, where $d\psi$ is
the change of the state and $P=|\psi><\psi|$ is a projection operator

Figure 8:  The evolution of the wave function $\psi(t)$ with $J=0.2$ and
$L=1$. $\text{Re}(\psi(t))$ and $\text{Im}(\psi(t))$ represent real part (a)
and imaginary part (b) of wave function $\psi(t)$ respectively. The notation
$|000\rangle\sim|111\rangle$ in the legend represents the eight components of
$\psi(t)$.

###  V.3 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. For two and
three qubits with interaction, however, it is very difficult to obtain such an
analytical solution, and thus we resort to the numerical method discussed in
IV.1.In this section, we give 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 V.2, but
there is an extra interaction term in the Hamiltonian

$H=B^{(1)}_{x}\sigma^{(1)}_{x}+B^{(1)}_{y}\sigma^{(1)}_{y}+B^{(2)}_{x}\sigma^{(2)}_{x}+B^{(2)}_{y}\sigma^{(2)}_{y}+J_{yy}\sigma^{(1)}_{y}\sigma^{(2)}_{y}$
(33)

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 $\psi(t)$ is no
longer separable, but becomes entangled. Following the work in the literature
Wootters (1998); Bennett  _et al._ (1996), we may define an entanglement index
$E_{p}$ for the two subsystems A and B, $E_{p}(|\psi\rangle_{AB})\equiv-
Tr_{A}(\rho_{A}log_{2}\rho_{A})=-Tr_{B}(\rho_{B}log_{2}\rho_{B})$, where
$\rho_{A}=Tr_{B}(|\psi\rangle_{AB}\langle\psi|)$ is the reduced density matrix
of the pure state $|\psi\rangle_{AB}$ over subsystem A, and $\rho_{B}$ has a
similar definition. Fig. 4(a-d) depict the entanglement of $\psi(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 target states) increases
first, reaches a maximum and then decreases back to zero. The entanglement
evolution profile for L=1, 3 (Fig. 4 a and b) is symmetric and the maximum is
exactly at the symmetry point, while for L =2, 4 (Fig. 4 c and d), the profile
is skewed. Given a suitable initial value of $\phi^{\prime}_{0}$, our
numerical method converges exponentially when approaching an optimal path
connecting the initial and target states (See Fig. 6(c)). Next, we keep the
initial state unchanged and set the target state as an entangled one with
$E_{p}=0.6061$. The wavefunction $\psi(t)$ connecting the two states is
computed with J=1, 0.9, 0.65, the entanglement evolution profile of which is
depicted in Fig. 5.

Case 2, we consider a two-qubit model, which is used in Wang  _et al._ (2015),
with the following Hamiltonian

$H=\sum\limits_{l,m}B_{m}^{(l)}\sigma_{m}^{(l)}+J\sum\limits_{m}\sigma_{m}^{(1)}\sigma_{m}^{(2)},$
(34)

where $\sigma_{m}^{(l)}~{}\text{with}~{}(m=x,y,z,l=1,2)$ are the Pauli
matrices for the $l^{th}$ qubit. the initial and final spin states are the
same as in section V.2 and the coupling constant J=1. The results obtained by
our method are portrayed in Fig. 6(a) and (b), where (a) depicts the
components of the observable and (b) shows that the control magnetic field
B(t).

Next, we compare the efficiency of finding such a connecting wavefunction
among different methods, such as the GRAPE and the “Geodesic-search” Wang  _et
al._ (2015). On the same machine and programming language, our method finds
such a connection wavefunction and the modulation magnetic field (shown in
Fig. 6(b)) in less than one minute. The GRAPE needs about 3.3 minutes and the
“Geodesic-search” needs about 18 minutes. Fig. 6(c) and (d) compare the
convergence upon iteration between our method and the GRAPE for the same
convergence criterion $10^{-5}$. Obviously, although the computation with our
method oscillates at the initial stage, it converges much more quickly and
displays good stability as well. However, the convergence of the GRAPE is
quite slow (with more than 500 iterations), which explains its relatively long
computation time. For more complex multi-qubit systems, we believe that these
advantages will become more obvious. In this case, the optimal time is
computed as $T=\frac{\sqrt{N}\pi\hbar}{\omega}$ 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 Margolus and Levitin
(1998); Levitin and Toffoli (2009). The Fubini-Study distance to the target
state seems to decrease linearly with time as shown in Fig. 7, and we note
that with the decrease of the energy $\omega$, the traveling time T increases.

Case 3, we study a model with three qubits for which the Hamiltonian is

$H=\sum\limits_{l,m}B_{m}^{(l)}\sigma_{m}^{(l)}+J(\sigma_{y}^{(1)}\sigma_{y}^{(2)}+\sigma_{y}^{(2)}\sigma_{y}^{(3)})$
(35)

where $\sigma_{m}^{(l)}~{}\text{with}~{}(m=x,y,z,l=1,2,3)$ are the Pauli
matrices for the $l^{th}$ qubit. We choose the initial state
$|+x\rangle^{(1)}\otimes|+y\rangle^{(2)}\otimes|+x\rangle^{(3)}$ and the final
state $|-x\rangle^{(1)}\otimes|-y\rangle^{(2)}\otimes|-x\rangle^{(3)}$ . The
obtained optimal wavefunction $\psi(t)$ connecting them with the interaction
strength J=0.2 and L=1 is depicted in Fig. 8.

Figure 9:  The wave function $\psi(t)$ for L=3 and $J_{yy}=1$. The red points
are invariant under the symmetry operation F. $|00>,|01>,|10>$ and $|11>$ are
the four components of $\psi(t)$, ‘Re’ and ‘Im’ represents the real and image
part.

###  V.4  The Symmetry Operation Discussion

Figure 10:  The function $\bm{B}(t)$ for L=2 ($a\sim d$) with $J_{yy}=0.05$
and L=3 ($e\sim h$)with $J_{yy}=0.2$. The black lines plot the controlled part
$\bm{B}(t)$ for qubit one and the red lines plot the controlled part
$\bm{B}(t)$ for the second qubit. The results of $c,~{}d,~{}g$ and $h$ were
obtained for the transformed Hamiltonian and states $H^{\prime}=F^{T}HF$,
$\psi_{0}^{\prime}=F\psi_{0}$, $\psi_{f}^{\prime}=F\psi_{f}$. Figure 11:  The
entanglement distribution of $\psi(t)$ (the black line) and $\psi^{\prime}(t)$
(the red line) with $J_{yy}=0.05$ and L=3, $\psi^{\prime}(t)$ is the symmetry
counterpart of $\psi(t)$ under the action $F$.

In this section, we check a symmetry operation of the Hamiltonian Eq. (33) and
ferret out possible mechanism underlying the phenomenon discussed in Part V.3
above. It is interesting to observe a discrete symmetry operation in
Hamiltonian (33) associated with

$F=\left[\begin{matrix}~{}~{}0&~{}~{}0&~{}~{}0&~{}~{}i\\\
~{}~{}0&~{}~{}0&~{}~{}i&~{}~{}0\\\ ~{}~{}0&-i&~{}~{}0&~{}~{}0\\\
-i&~{}~{}0&~{}~{}0&~{}~{}0\end{matrix}\right]$ (36)

which connects the initial state $\psi_{0}$ and the target state $\psi_{f}$,
i.e. $F\psi_{0}=-\psi_{f}$ in the current case. Obviously, the magnetic field
$\bm{B}(t)$ in the Hamiltonian is artificially controlled, which should be
changed accordingly (by $F^{T}HF$) to complete the symmetry operation, as
confirmed in Fig. 10. For the case of L=1,3, especially, the middle state
$\psi_{mid}$ (at t=0.5) remains invariant under this symmetry operation,
namely $F\psi_{mid}=\psi_{mid}$. In fact, to search the symmetry, a clue could
be directly taken from Fig. 9, where the red points represent the middle state
$\psi_{mid}=(0.09483-0.08982i,~{}-0.69486-0.02312i,~{}-0.02309+0.69443i,~{}-0.08993-0.09488i)^{t}$.
It is easy to check that the time course of the wave function $\psi_{t}$ is
also symmetric about the middle state $\psi_{mid}$. From the view of the
controlled part $\bm{B}(t)$ of $H$, take the cases of L=2 and 3 as an example.
Fig.10 $a,~{}b~{}(L=2)$ and $e,~{}f~{}(L=3)$ depict the time course of
$\bm{B}(t)$ from $\psi_{0}$ to $\psi_{f}$ while Fig.10 $c,~{}d~{}(L=2)$ and
$g,~{}h~{}(L=3)$ plot the time course of $\bm{B}(t)$ from
$\psi_{0}^{\prime}=F\psi_{0}$ to $\psi_{f}^{\prime}=F\psi_{f}$ with the
transformed Hamiltonian $H^{\prime}=F^{T}HF$. Apparently, Fig.10 $e$ and $f$
coincide with $g$ and $h$ respectively up to a minus sign in the relevant
component. In addition, the evolution of $B^{(1)}_{x}$ and $B^{(2)}_{y}$ is
axisymmteric with respect to $t=0.5$ and the evolution profiles of
$B^{(1)}_{y}$ and $B^{(2)}_{x}$ have a center symmetry with respect to (0.5,
0). All the above indicates that the maximum of $E_{p}$ is associated with a
fixed point of the operation $F$ for the case L=1, 3.

In the cases L=2, 4, the symmetry in $\psi(t)$ and $\bm{B}(t)$ discussed
previously disappear. Taking L=2 for an example, Fig.10 $a$ and $b$ do not
coincide with Fig.10 $c$ and $d$ and the maximum value of $E_{p}$ is not at
t=0.5. On the other hand, as the symmetry partner $\psi^{\prime}_{t}$ of the
wave function $\psi_{t}$ is obtained with the symmetry operation $F$, even if
the entanglement evolution of $\psi^{\prime}_{t}$ (see Fig. 11, the red line)
is asymmetric, those of $\psi^{\prime}(t)$ and $\psi(t)$ are axially symmetric
with respect to each other about t=0.5.

Figure 12:  The entanglement evolution of $\psi(t)$ with $L=1$ and
$J_{yy}=0.4$

in Eq. (33), the angle between $\psi_{0}$ and $\psi_{f}$ is from $\pi$ to
$0.7\pi$.

In the case where $\psi_{0}$ and $\psi_{f}$ do not flip to each other, e.g.
the angle between $\psi_{0}$ and $\psi_{f}$ is less than $\pi$, we then have
$F\psi_{0}\neq-\psi_{f}$ (take the case $L=1$ as an example), which indicates
that the entanglement evolution is asymmetric (shown in Fig. 12). However, the
maximum entanglement still exists which is unique but not at $t=0.5$.

##  VI SUMMARY AND DISCUSSION

In this work, we derive a general set of differential equations (Eq. (25)) for
an optimal quantum control, which is valid for single or multiple qubits with
or without interaction. In the derivation and application of Eq. (25), 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,
for the even oscillation number symmetry disappears in the profile but the
symmetry operation gives another optimal path – the symmetry partner of the
original solution.

With 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.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
under Grants No. 11775035, and also by the Fundamental Research Funds for the
Central Universities with Contract No.2019XD-A10.

## References

  * Feynman (1982) R. P. Feynman, Simulating physics with computers, Int. J. Theor. Phys. 21, 467 (1982).
  * SHOR (1994) P. SHOR, Algorithms for quantum computation : Discrete logarithms and factoring, in  _Proceedings of 35th Annual Symposium on Foundations of Computer Scienece_ (1994).
  * Boixo  _et al._ (2014) S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Evidence for quantum annealing with more than one hundred qubits, Nat. Phys. 10, 218 (2014).
  * del Campo  _et al._ (2013) A. del Campo, I. L. Egusquiza, M. B. Plenio, and S. F. Huelga, Quantum speed limits in open system dynamics, Phys. Rev. Lett. 110, 050403 (2013).
  * Deffner and Lutz (2013) S. Deffner and E. Lutz, Quantum speed limit for non-markovian dynamics, Phys. Rev. Lett. 111, 010402 (2013).
  * Taddei  _et al._ (2013) M. M. Taddei, B. M. Escher, L. Davidovich, and R. L. de Matos Filho, Quantum speed limit for physical processes, Phys. Rev. Lett. 110, 050402 (2013).
  * Mirkin  _et al._ (2016) N. Mirkin, F. Toscano, and D. A. Wisniacki, Quantum-speed-limit bounds in an open quantum evolution, Phys. Rev. A 94, 052125 (2016).
  * Schulte-Herbruggen  _et al._ (2005) T. Schulte-Herbruggen, A. Sporl, N. Khaneja, and S. J. Glaser, Optimal control-based efficient synthesis of building blocks of quantum algorithms: A perspective from network complexity towards time complexity, Phys. Rev. A 72, 042331 (2005).
  * Nielsen and Chuang (2011) M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information: 10th anniversary edition (Cambridge University Press, 2011) 10th ed.
  * Nielsen  _et al._ (2006a) M. A. Nielsen, M. R. Dowling, M. Gu, and A. C. Doherty, Quantum computation as geometry, Science 311, 1133 (2006a).
  * Nielsen  _et al._ (2006b) M. A. Nielsen, M. R. Dowling, M. Gu, and A. C. Doherty, Optimal control, geometry, and quantum computing, Phys. Rev. A 73, 062323 (2006b).
  * Khaneja  _et al._ (2005) N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Optimal control of coupled spin dynamics: design of nmr pulse sequences by gradient ascent algorithms, Journal of Magnetic Resonance 172, 296 (2005).
  * Maximov  _et al._ (2008) I. I. Maximov, Z. To$\check{s}$ner, and N. C. Nielsen, Optimal control design of nmr and dynamic nuclear polarization experiments using monotonically convergent algorithms, J. Chem. Phys. 128, 184505 (2008).
  * Maday and Turinici (2003) Y. Maday and G. Turinici, New formulations of monotonically convergent quantum control algorithms, J. Chem. Phys. 118, 8191 (2003).
  * Zhu and Rabitz (1998) W. Zhu and H. Rabitz, A rapid monotonically convergent iteration algorithm for quantum optimal control over the expectation value of a positive definite operator, J. Chem. Phys. 109, 385 (1998).
  * Shi and Rabitz (1990) S. Shi and H. Rabitz, Quantum mechanical optimal control of physical observables in microsystems, J. Chem. Phys. 92, 364 (1990).
  * Khaneja  _et al._ (2001) N. Khaneja, R. Brockett, and S. J. Glaser, Time optimal control in spin systems, Phys. Rev. A 63, 032308 (2001).
  * Boozer (2012) A. D. Boozer, Time-optimal synthesis of su(2) transformations for a spin-1/2 system, Phys. Rev. A 85, 012317 (2012).
  * Hegerfeldt (2013) G. C. Hegerfeldt, Driving at the quantum speed limit: Optimal control of a two-level system, Phys. Rev. Lett. 111, 260501 (2013).
  * Carlini  _et al._ (2006) A. Carlini, A. Hosoya, T. Koike, and Y. Okudaira, Time-optimal quantum evolution, Phys. Rev. Lett. 96, 60503 (2006).
  * Meiss (1992) J. D. Meiss, Symplectic maps, variational principles, and transport, Rev. Mod. Phys. 64, 795 (1992).
  * Bunimovich (1995) L. A. Bunimovich, Variational principle for periodic trajectories of hyperbolic billiards, Chaos 5, 349 (1995).
  * Lan and Cvitanovic (2004) Y. Lan and P. Cvitanovic, Variational method for finding periodic orbits in a general flow, Phys. Rev. E 69, 016217 (2004).
  * Ghoussoub and Moameni (2007) N. Ghoussoub and A. Moameni, Selfdual variational principles for periodic solutions of hamiltonian and other dynamical systems, Comm. in PDE 32, 771 (2007).
  * Dong and Lan (2014) C. Dong and Y. Lan, A variational approach to connecting orbits in nonlinear dynamical systems, Phys. Lett. A 378, 705 (2014).
  * Wang  _et al._ (2018) D. Wang, P. Wang, and Y. Lan, Accelerated variational approach for searching cycles, Phys. Rev. E 98, 042204 (2018).
  * Carlini  _et al._ (2007) A. Carlini, A. Hosoya, T. Koike, and Y. Okudaira, Time-optimal unitary operations, Phys. Rev. A 75, 042308 (2007).
  * Carlini  _et al._ (2008) A. Carlini, A. Hosoya, T. Koike, and Y. Okudaira, Time optimal quantum evolution of mixed states, J. Phys. A: Math. Theor. 41, 045303 (2008).
  * Carlini  _et al._ (2014) A. Carlini, A. Mari, and V. Giovannetti, Time-optimal thermalization of single-mode gaussian states, Phys. Rev. A 90, 052324 (2014).
  * Wang  _et al._ (2015) X. Wang, M. Allegra, K. Jacobs, S. Lloyd, C. Lupo, and M. Mohseni, Quantum brachistochrone curves as geodesics: Obtaining accurate minimum-time protocols for the control of quantum systems, Phys. Rev. Lett. 114, 170501 (2015).
  * Wang  _et al._ (2017) X. Wang, M. Allegra, K. Jacobs, S. Lloyd, C. Lupo, and M. Mohseni, Time-optimal quantum control via differential geometry, in  _Advances in Photonics of Quantum Computing, Memory, and Communication X_, Vol. 10118, edited by Z. U. Hasan, P. R. Hemmer, H. Lee, and A. L. Migdall, International Society for Optics and Photonics (SPIE, 2017) pp. 53 – 59\.
  * Wootters (1998) W. K. Wootters, Entanglement of formation of an arbitrary state of two qubits, Phys. Rev. Lett. 80, 2245 (1998).
  * Bennett  _et al._ (1996) C. H. Bennett, H. J. Bernstein, S. Popescu, and B. Schumacher, Concentrating partial entanglement by local operations, Phys. Rev. A 53, 2046 (1996).
  * Margolus and Levitin (1998) N. Margolus and L. B. Levitin, The maximum speed of dynamical evolution, Physica D: Nonlinear Phenomena 120, 188 (1998).
  * Levitin and Toffoli (2009) L. B. Levitin and T. Toffoli, Fundamental limit on the rate of quantum dynamics: The unified bound is tight, Phys. Rev. Lett. 103, 160502 (2009).

## Appendix A

The variation of the Hamiltonian $H$ could now be taken as multivariate change
with respect to $\\{\xi_{j}\\}_{j=1,2,...m}$. Let’s respectively check
different parts of the action integral Eq. (2). First, check one term of the
integrand

$\displaystyle I$ $\displaystyle=$
$\displaystyle\langle\phi|H|\psi\rangle+c.c.$ (37) $\displaystyle=$
$\displaystyle\sum_{j}\langle\phi|A_{j}|\psi\rangle\xi_{j}(t)+c.c.$
$\displaystyle=$ $\displaystyle\sum_{j}\xi_{j}(t)D_{j}\,,$

where $D_{j}=\langle\phi|A_{j}|\psi\rangle+c.c.$ are real numbers, thus the
Hamiltonian variation of the term I is

$\delta I=\sum\limits_{j}D_{j}\delta\xi_{j}(t).$ (38)

Next, let’s consider the second part

$II=g(\bm{\xi},\lambda,\lambda^{\prime})=\lambda(\frac{1}{2}\sum\limits_{j=1}^{m}\xi_{j}^{2}-\omega^{2})+\sum\limits_{j=n+1}^{m}\lambda^{{}^{\prime}}_{j}(\xi_{j}-Q_{i})\,.$

Let $G_{j}=\frac{\delta g(\bm{\xi},\lambda)}{\delta\xi_{j}}$. We obtain the
variation of the II

$\delta II=\sum\limits_{j}G_{j}\delta\xi_{j}.$ (39)

Third, before treating first term of the Eq. (2), we need to do some
preparation. Here, we define:

$F_{jk}=\langle A_{j}A_{k}+A_{k}A_{j}\rangle-2\langle A_{j}\rangle\langle
A_{k}\rangle\,,$ (40)

which is a symmetric tensor, and the energy variance could be written as:

$\displaystyle\Delta E^{2}$ $\displaystyle=$ $\displaystyle\langle
H^{2}\rangle-\langle H\rangle^{2}$ (41) $\displaystyle=$
$\displaystyle\sum_{jk}\xi_{j}\xi_{k}\langle
A_{j}A_{k}\rangle-\xi_{j}\xi_{k}\langle A_{j}\rangle\langle A_{k}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{2}\sum_{jk}F_{jk}\xi_{j}\xi_{k}\,,$

leading to

$\delta\Delta E^{2}=\delta\dfrac{\bm{\xi}\cdot
F\cdot\bm{\xi}}{2}=\sum\limits_{jk}F_{jk}\xi_{k}\delta\xi_{j}\,,$ (42)

which gives

$\displaystyle\delta\frac{1}{\Delta E}$ $\displaystyle=$
$\displaystyle-\frac{1}{\Delta E^{2}}\frac{1}{2\Delta E}\delta\Delta E^{2}$
(43) $\displaystyle=$ $\displaystyle-\frac{1}{2\Delta
E^{3}}\sum_{jk}F_{jk}\xi_{k}\delta\xi_{j}\,.$

So, the variation of the first term of Eq. (2) is

$\displaystyle\delta\frac{\sqrt{\langle\dot{\psi}|(1-P)|\dot{\psi}\rangle}}{\Delta
E}$
$\displaystyle=\sqrt{\langle\dot{\psi}|(1-P)|\dot{\psi}\rangle}\delta\dfrac{1}{\Delta
E}{}$ (44) $\displaystyle=-\dfrac{1}{2\Delta
E^{2}}\sum\limits_{jk}F_{jk}\xi_{k}\delta\xi_{j}\,,$

where Eq. (3) and (41) have been used. With Eqs. (38), (39)and (44), the
variation with respect to $\xi_{j}$ gives

$-\dfrac{1}{2\Delta E^{2}}\sum\limits_{k}F_{jk}\xi_{k}+D_{j}+G_{j}=0$ (45)

where

$G_{j}=\left\\{\begin{aligned} &\lambda\xi_{j},&(j\leq n).\\\
&\lambda\xi_{j}-\lambda^{{}^{\prime}}_{j},&(n<j\leq m).\end{aligned}\right.$
(46)

## Appendix B

Multiplyingv $\langle\psi|A_{j}$ on both sides of Eq. (4) results in

$\dfrac{i}{2\Delta E^{2}}(\langle A_{j}\dot{H}\rangle-\langle
A_{j}\rangle\langle\dot{H}\rangle)+i(\langle A_{j}H\rangle-\langle
A_{j}\rangle\langle H\rangle)\dfrac{d}{dt}\dfrac{1}{2\Delta
E^{2}}\\!-i\frac{d}{dt}\langle\psi|A_{j}|\phi\rangle+\langle\psi|[A_{j},H]|\phi\rangle\\!=\\!0.$
(47)

The derivative of $\dfrac{1}{2\Delta E^{2}}$ is

$\dfrac{d}{dt}\dfrac{1}{2\Delta E^{2}}=-\dfrac{1}{\Delta
E^{3}}\dfrac{1}{2\Delta E}\dfrac{d}{dt}\Delta E^{2}.$ (48)

Therefore, with Eq. (40) and Eq. (41) and using
$\bm{\xi}\cdot\dot{F}\cdot\bm{\xi}=\langle[H,H^{2}]\rangle=0$, we have

$\dfrac{d}{dt}\Delta E^{2}=\bm{\xi}\cdot
F\cdot\dot{\bm{\xi}}=\tilde{\bm{F}}\cdot\dot{\bm{\xi}},$ (49)

where $\tilde{\bm{F}}=\bm{\xi}\cdot F$. And the Eq.(48) can be rewritten as

$\dfrac{d}{dt}\dfrac{1}{2\Delta E^{2}}=-\dfrac{1}{2\Delta
E^{2}}\dfrac{\tilde{\bm{F}}\cdot\dot{\bm{\xi}}}{\Delta E^{2}}.$ (50)

Substituting Eq. (50) into Eq. (47), we arrive at

$\dfrac{d}{dt}D_{j}=\dfrac{1}{2\Delta
E^{2}}(\sum\limits_{k}F_{jk}\dot{\xi}_{k}-\dfrac{\tilde{\bm{F}}\cdot\dot{\bm{\xi}}}{\Delta
E^{2}}\tilde{F_{j}})+\sum_{k}C_{jk}D_{k}\,,$ (51)

where $C_{jk}$ is a real number and satisfies

$[H,A_{j}]=-i\sum\limits_{k}C_{jk}A_{k}.$ (52)

##  Appendix C  $\dfrac{d}{dt}\bm{D}\cdot\bm{\xi}=0$

$\displaystyle\dfrac{d}{dt}\bm{D}\cdot\bm{\xi}$
$\displaystyle=\bm{D}\cdot\dfrac{d\bm{\xi}}{dt}+\dfrac{d\bm{D}}{dt}\cdot\bm{\xi}$
$\displaystyle=\sum_{j}\dot{\xi}_{j}D_{j}+\dfrac{\xi_{j}}{2\Delta
E^{2}}\left(\sum_{k}F_{jk}\dot{\xi}_{k}-\dfrac{\tilde{F}\cdot\dot{\xi}}{\Delta
E^{2}}\tilde{F_{j}}\right)+\sum_{k}\xi_{j}C_{jk}D_{k}$
$\displaystyle=\dot{\bm{\xi}}\cdot\bm{D}+\dfrac{1}{2\Delta
E^{2}}(\tilde{\bm{F}}\cdot\dot{\bm{\xi}}-\tilde{\bm{F}}\cdot\dot{\bm{\xi}}\dfrac{\tilde{\bm{F}}\cdot{\bm{\xi}}}{\Delta
E^{2}})+\bm{\xi}\cdot\bm{C}\cdot\bm{D}$
$\displaystyle=\dot{\bm{\xi}}\cdot\bm{D}-\dfrac{\tilde{\bm{F}}\cdot\dot{\bm{\xi}}}{2\Delta
E^{2}}+\bm{\xi}\cdot\bm{C}\cdot\bm{D}\,,$ (53)

where $\tilde{\bm{F}}\cdot\bm{\xi}=2\Delta E^{2}$. As $D_{j}=\dfrac{1}{2\Delta
E^{2}}\sum_{k}F_{jk}\xi_{k}-\lambda\xi_{j}-\lambda_{j}^{\prime}$, which gives

$\dfrac{d}{dt}\bm{D}\cdot\bm{\xi}=\dfrac{\tilde{\bm{F}}\cdot\dot{\bm{\xi}}}{2\Delta
E^{2}}-\dfrac{\tilde{\bm{F}}\cdot\dot{\bm{\xi}}}{2\Delta
E^{2}}+\bm{\xi}\cdot\bm{C}\cdot\bm{D}=0$ (54)

with $\bm{\xi}\cdot\dot{\bm{\xi}}=0$,
$\bm{\lambda}^{\prime}\cdot\dot{\bm{\xi}}=\bm{\lambda}^{\prime}\cdot\bm{\dot{u}}=0$
and $\bm{\xi}\cdot\bm{C}=0$ (Seeing Appendix E).

##  Appendix D $\dot{\bm{F}}=\bm{C}\cdot\bm{F}+\bm{F}\cdot\bm{C}^{T}$

Acorrding to the Eq. (40):

$F_{jk}=\langle A_{j}A_{k}+A_{k}A_{j}\rangle-2\langle A_{j}\rangle\langle
A_{k}\rangle\,,$ (55)

we have

$\displaystyle\dot{F}_{jk}=$ $\displaystyle i\langle
HA_{j}A_{k}+HA_{k}A_{j}-A_{j}A_{k}H-A_{k}A_{j}H\rangle-2i\langle
HA_{j}-A_{j}H\rangle\langle A_{k}\rangle-2i\langle A_{j}\rangle\langle
HA_{k}-A_{k}H\rangle$ $\displaystyle=$ $\displaystyle
i\langle[H,A_{j}]A_{k}+A_{k}[H,A_{j}]+A_{j}[H,A_{k}]+[H,A_{k}]A_{j}\rangle-2i\langle[H,A_{j}]\rangle\langle
A_{k}\rangle-2i\langle A_{j}\rangle\langle[H,A_{k}]\rangle$ $\displaystyle=$
$\displaystyle\sum_{l}C_{jl}(\langle A_{l}A_{k}+A_{k}A_{l}\rangle-2\langle
A_{l}\rangle\langle A_{k}\rangle)+\sum_{l}C_{kl}(\langle
A_{j}A_{l}+A_{l}A_{j}\rangle-2\langle A_{j}\rangle\langle A_{l}\rangle)$
$\displaystyle=$ $\displaystyle\sum_{l}C_{jl}F_{lk}+\sum_{l}C_{kl}F_{lj}\,,$

then

$\dot{\bm{F}}=\bm{C}\cdot\bm{F}+\bm{F}\cdot\bm{C}^{T}\,.$ (56)

##  Appendix E $\bm{C}\cdot\bm{\xi}=\bm{\xi}\cdot\bm{C}=0$

Here we demonstrate that $\bm{C}\cdot\bm{\xi}=\bm{\xi}\cdot\bm{C}=0$ is
universally valid with respect to spin algebra. In the case of
$H\in\mathcal{A}$ i.e. $\bm{\xi}=\bm{u}$, it yields
$\bm{C}\cdot\bm{u}=\bm{u}\cdot\bm{C}=0$.

###  E.1 one qubit

Consider $A_{i},A_{j}$ in a qubit and suppose that
$\\{A_{j}\\}_{j=1\,,\cdots\,,M}$ satisfy the spin algebra

$[A_{n}\,,A_{j}]=\sum_{k}\dfrac{1}{\sqrt{TrI}}~{}i\epsilon_{njk}A_{k}\,,$ (57)

and $I$ is $2^{q}$ dimensional unit matrix (q is the number of qubits) and
$A_{k}=\frac{1}{\sqrt{2}}\sigma_{k}$ for $q=1$. Hence

$\sum_{k}C_{jk}A_{k}=i[H,A_{j}]=\sum_{n}i\xi_{n}[A_{n},A_{j}]=\sum_{n,k}-\dfrac{1}{\sqrt{TrI}}\epsilon_{njk}\xi_{n}A_{k}$

which yields $C_{jk}=\sum_{n}-\dfrac{1}{\sqrt{TrI}}\epsilon_{njk}\xi_{n}$.
Thus $\bm{\xi}\cdot\bm{C}=-\bm{C}\cdot\bm{\xi}=0$.

###  E.2 many qubits with and without entanglement

We will prove that $\bm{\xi}\cdot\bm{C}=\bm{C}\cdot\bm{\xi}=0$ remains valid.
Below, we denote the unit matrix $I_{0}$ with $\sigma_{0}$ and extend our
index to include 0.

Lemma $C_{jk}=\sum_{i}\Gamma_{ijk}\xi_{i}$, where $\Gamma$ is an antisymmetric
tensor.

Proof First let us check some properties of Pauli matrices. With the
commutation and anti-commutation relations
$\\{\sigma_{i},\sigma_{j}\\}=2\delta_{ij}I(i,j=1,2,3)$, it is easy to check
even if $i,j=0,1,2,3$

$\sigma_{i}^{(\alpha)}\sigma_{j}^{(\alpha)}=\sum_{k}(\theta_{ijk}^{(\alpha)}+i\epsilon_{ijk}^{(\alpha)})\sigma_{k}^{(\alpha)}\,,$
(58)

where $\sigma_{0}=I_{0}$ which is two-dimensional unit matrix and $\alpha$ is
the order index of qubits,

$\epsilon_{ijk}^{(\alpha)}=\left\\{\begin{aligned}
&1,~{}ijk\in\\{123,231,312\\}\\\ &-1,~{}ijk\in\\{321,213,132\\}\\\
&0,~{}\mathrm{otherwise}\end{aligned}\right.\,,$ (59)

and

$\theta_{ijk}^{(\alpha)}=\left\\{\begin{aligned}
&1,~{}&\mathrm{one~{}index~{}is~{}0,~{}the~{}other~{}two~{}equal}\\\
&0,~{}&\mathrm{otherwise}\end{aligned}\right.\,,$ (60)

which is totally symmetric with respect to any permutation of its subscripts.
Now consider the commutator of N-fold entanglement
$[A_{\mu},A_{\nu}]=[\dfrac{\sigma_{\mu_{1}}^{(1)}\sigma_{\mu_{2}}^{(2)}...\sigma_{\mu_{N}}^{(N)}}{\sqrt{TrI}},\dfrac{\sigma_{\nu_{1}}^{(1)}\sigma_{\nu_{2}}^{(2)}...\sigma_{\nu_{N}}^{(N)}}{\sqrt{TrI}}]$.
Here the factor $\dfrac{1}{\sqrt{TrI}}$ is introduced to guarantee that
$TrA_{j}^{2}=1$ and
$\\{\mu_{j},\nu_{j},\gamma_{j}\\}_{j=1,2,...,N}=\\{0,1,2,3\\}$. Without loss
of generality,

$TrI[A_{\mu},A_{\nu}]\\!\\!=\\!\\!\prod_{j}\sigma_{\mu_{j}}^{(j)}\sigma_{\nu_{j}}^{(j)}-\prod_{j}\sigma_{\nu_{j}}^{(j)}\sigma_{\mu_{j}}^{(j)}\\!=\\!\prod_{j}\\!\sum_{\gamma_{j}}(\theta_{\mu_{j}\nu_{j}\gamma_{j}}^{(j)}+i\epsilon_{\mu_{j}\nu_{j}\gamma_{j}}^{(j)})\sigma_{\gamma_{j}}^{(j)}-\\!\prod_{j}\\!\sum_{\gamma_{j}}(\theta_{\nu_{j}\mu_{j}\gamma_{j}}^{(j)}+i\epsilon_{\nu_{j}\mu_{j}\gamma_{j}}^{(j)})\sigma_{\gamma_{j}}^{(j)}\,.$
(61)

Define the polynomial coefficient of the $k$-th power of $\epsilon$ as
$c_{k}(\mu,\nu,\gamma)$ , which satisfies

$\sum_{k}c_{k}p_{k}(\mu,\nu,\gamma)=\prod_{j}(\theta_{\mu_{j}\nu_{j}\gamma_{j}}^{(j)}+i\epsilon_{\mu_{j}\nu_{j}\gamma_{j}}^{(j)})-\prod_{j}(\theta_{\mu_{j}\nu_{j}\gamma_{j}}^{(j)}-i\epsilon_{\mu_{j}\nu_{j}\gamma_{j}}^{(j)})\,,$
$[A_{\mu},A_{\nu}]=\sum_{\gamma}P_{\mu\nu\gamma}A_{\gamma}\,,$ (62)

where $p_{k}(\mu,\nu,\gamma)$ is the polynomial function of $\epsilon$ and
$A_{\gamma}=\dfrac{1}{\sqrt{TrI}}\prod\limits_{j}\sigma_{\gamma_{j}}^{(j)}$,
$P_{\mu\nu\gamma}=\dfrac{1}{\sqrt{TrI}}\sum\limits_{k}c_{k}p_{k}(\mu,\nu,\gamma)$
is a pure imaginary number and holds the series of odd powers of $\epsilon$
(It is easy to check that $c_{2k}=1-(-1)^{2k}\equiv 0,k\in\mathbb{N}$). Here
we have utilized the symmetry of $\theta$ and the anti-symmetry of $\epsilon$.

Now commutate two of the indices ${\mu,\nu,\gamma}$ of $P$. Notice that $P$ is
written as series of odd powers of $\epsilon$ and each pair of
$\\{\mu_{j},\nu_{j}\\}$ anti-commutes. Hence P has the same anti-symmetry as
$\epsilon$. Given the definition of $\bm{C}$

$\sum\limits_{\gamma}C_{\nu\gamma}A_{\gamma}=i[H,A_{\nu}]\,,$ (63)

we have $C_{jk}=\sum\limits_{n}iP_{njk}\xi_{n}$.

According to the lemma,
$\bm{C}\cdot\bm{\xi}=\sum\limits_{k}\sum\limits_{n}iP_{njk}\xi_{n}\xi_{k}=-\sum\limits_{k}\sum\limits_{n}iP_{kjn}\xi_{n}\xi_{k}=-\bm{C}\cdot\bm{\xi}$.
Thus $\bm{C}\cdot\bm{\xi}=0$. Furthermore, $\bm{C}$ is anti-symmetric in
consistency with $P$. Hence
$\bm{\xi}\cdot\bm{C}=\sum\limits_{j}\xi_{j}C_{jk}=-\sum\limits_{j}\xi_{j}C_{kj}=-\bm{C}\cdot\bm{\xi}=0$.

##  Appendix F The geometric properties of spin vector evolution

Lemma: $\bm{D}^{\prime}\cdot\langle\bm{\sigma}\rangle=0$ for single qubit.

Proof: From Eq. (26)in the main text, we have
$Re(\langle\phi^{\prime}|\psi\rangle)=Re(\langle\phi|\psi\rangle)=0$.

Considering the projection property of Pauli matrices

$\bm{\sigma}\cdot\hat{\bm{p}}=2|\psi\rangle\langle\psi|-I\\\ $ (64)

where $|\psi\rangle$ is a normalized eigenvector of $\bm{\sigma}$ in the
direction of $\hat{\bm{p}}=\langle\bm{\sigma}\rangle$, and $I$ is the unit
matrix, then we have that

$\displaystyle\bm{D}^{\prime}\cdot\langle\bm{\sigma}\rangle$
$\displaystyle=2Re(\langle\phi^{\prime}|\bm{\sigma}|\psi\rangle)\cdot\langle\bm{\sigma}\rangle$
$\displaystyle=2Re(\langle\phi^{\prime}|\bm{\sigma}\cdot\hat{\bm{p}}|\psi\rangle)$
$\displaystyle=2Re(2\langle\phi^{\prime}|\psi\rangle\langle\psi|\psi\rangle-\langle\phi^{\prime}|\psi\rangle)$
$\displaystyle=0\,.$ (65)

Similarly, we can easily verify that

$\bm{D}\cdot\langle\bm{\sigma}\rangle=2Re(\langle\phi|\bm{\sigma}|\psi\rangle)\cdot\langle\bm{\sigma}\rangle=0.$
(66)

## Appendix G

In this section, we’re going to elaborate on the rotation matrix
$R_{\bm{\Omega}^{\prime}}(\theta)=R^{-1}_{x}(\Psi)R_{y}(\theta)R_{x}(\Psi)$ in
the single qubit case. According to III.1, the $\bm{\Omega}$ lies in the $z$
direction and perpendicular to the $\bm{\omega}_{eff}$, and
$R_{\bm{\Omega}^{\prime}}(\theta)$ depicts the rotation of $\psi(t)$ arounds a
fixed axis with $\bm{\Omega}^{\prime}=\bm{\omega}_{eff}-\bm{\Omega}$.
Obviously, the $\bm{\omega}_{eff}$ lies on the $x-y$ plane, and the direction
of $\bm{\omega}_{eff}$ is taken as the $-y$ axis for simplicity. Then the
rotation $R_{\bm{\Omega}^{\prime}}(\theta)$ can be depicted as follow, first,
rotating $\psi$ around $x$-axis by $\Psi$ to make $\bm{\Omega}^{\prime}$ and
$y$-aixs coincide $(R_{x}(\Psi))$, then, rotating $\psi$ around $y$-axis by
$\theta$ $(R_{y}(\theta))$, finally, rotating $\psi$ around $x$-axis back by
$\Psi$ $(R^{-1}_{x}(\Psi))$, where $\Psi$ is the angle between
$\Omega^{\prime}$ and $\bm{\omega}_{eff}$ and
tan$(\Psi)=\frac{|\bm{\Omega}|}{|\bm{\omega}_{eff}|}$.