# S3HQA: A Three-Stage Approach for Multi-hop Text-Table Hybrid Question
Answering

Fangyu Lei1,2 Xiang Li1,2 Yifan Wei1,2
Shizhu He1,2 Yiming Huang1,2 Jun Zhao1,2 Kang Liu1,2
1The Laboratory of Cognition and Decision Intelligence for Complex Systems
Institute of Automation  Chinese Academy of Sciences
2School of Artificial Intelligence  University of Chinese Academy of Sciences
{leifangyu2022, lixiang2022, weiyifan2021<EMAIL_ADDRESS>
{shizhu.he, jzhao<EMAIL_ADDRESS>
###### Abstract

Answering multi-hop questions over hybrid factual knowledge from the given
text and table (TextTableQA) is a challenging task. Existing models mainly
adopt a retriever-reader framework, which have several deficiencies, such as
noisy labeling in training retriever, insufficient utilization of
heterogeneous information over text and table, and deficient ability for
different reasoning operations. In this paper, we propose a three-stage
TextTableQA framework S3HQA, which comprises of _retriever_ , _selector_ , and
_reasoner_. We use a _retriever with refinement training_ to solve the noisy
labeling problem. Then, a _hybrid selector_ considers the linked relationships
between heterogeneous data to select the most relevant factual knowledge. For
the final stage, instead of adapting a reading comprehension module like in
previous methods, we employ a _generation-based reasoner_ to obtain answers.
This includes two approaches: a row-wise generator and an LLM prompting
generator (first time used in this task). The experimental results demonstrate
that our method achieves competitive results in the few-shot setting. When
trained on the full dataset, our approach outperforms all baseline methods,
ranking first on the HybridQA
leaderboard.111https://codalab.lisn.upsaclay.fr/competitions/7979.

##  1 Introduction

Question answering systems devote to answering various questions with the
evidence located in the structured knowledge base (e.g., table) (Pasupat and
Liang, 2015; Yu et al., 2018) or unstructured texts (Rajpurkar et al., 2016).
Considering that many questions need to utilize multiple sources of knowledge
jointly in real-world applications, the hybrid form of question answering over
texts and tables (TextTableQA) has been proposed and attracted more and more
attention (Chen et al., 2020b, a; Zhu et al., 2021; Chen et al., 2021; Zhao et
al., 2022; Wang et al., 2022a). Fact reasoning (Chen et al., 2020a, b) is a
critical question type of TextTableQA. It requires jointly using multiple
evidence from tables and texts to reasoning the answers with different
operations, such as correlation (e.g., multi-hop) and aggregation (e.g.,
comparison). Hyperlinks among some table cells and linked passages are
essential resources to establish their relationship and support the retrieval
and reasoning for multi-hop questions. As shown in Figure 1, answering a
complex question Q1 requires jointly reasoning from textual evidence (P1) to
table evidence ([R2, Place]) and then to other table evidence ([R2, Athlete]).

Figure 1: The examples of HybridQA. Figure 2: An overview of S3HQA framework.
The retrieval stage is divided into two steps. The hybrid selector considers
the linked relationships between heterogeneous data to select the most
relevant factual knowledge.

Existing methods consist of two main stages: _retriever_ and _reader_ (Chen et
al., 2020b; Feng et al., 2022). The _retriever_ filters out the cells and
passages with high relevance to the question, and then the _reader_ extracts a
span from the retrieval results as the final answer. However, current methods
with two stages still have three limitations as follows.

1) Noisy labeling for training retriever. Existing retrieval methods usually
ignore the weakly supervised answer annotation (Chen et al., 2020b; Wang et
al., 2022b; Feng et al., 2022). For the Q2 of Figure 1, we cannot know the
specific location of the hybrid evidence, only given the final answer "1960".
Therefore, there is a lot of pseudo-true evidence labeled (Marked in green)
automatically by string matching, which introduces a lot of evidence noise.

2) Insufficient utilization of heterogeneous information. After retrieval,
existing methods selected a particular cell or passage for reading to extract
the final answer (Chen et al., 2020b; Wang et al., 2022b). As for Q1 in Figure
1, previous models were more likely to choose P1 or the coordinates [R2,Place]
to extract the answer. However, these methods seldomly used the hybrid
information of table schema and cell-passage hyperlinks, which is the key
factor in answering multi-hop questions.

3) Deficient ability for different reasoning operations. Previous methods
(Eisenschlos et al., 2021; Kumar et al., 2021; Wang et al., 2022b) mainly used
an extraction module to obtain answers, which cannot support knowledge
reasoning that requires comparison, calculation, and other operations.

In this paper, we propose a three-stage approach S3HQA to solve the above
problems. (1) Retriever with Refinement Training, we propose a two-step
training method, splitting the training data into two parts, so that the noise
in the retrieval phase can be alleviated. (2) Hybrid Selector has been
proposed and selects supporting facts with different granularity and resources
depending on the question type. By considering the hybrid data of tables and
text, this paper proposes a hybrid selection algorithm that can effectively
utilize the heterogeneous information of tables and passages. (3) Generation-
based reasoner utilizes a generation-based model for addressing different
question types. The model allows better aggregation of information on the
input side, which not only have better multi-hop reasoning capabilities but
also be able to handle comparison and counting questions. Furthermore, we are
the first to use the LLM in-context learning approach for table-text hybrid
question-answering tasks.

We evaluate our proposed model on the challenging TextTableQA benchmark
HybridQA. The empirical results show that our approach outperforms all the
existing models222We released the source code at
https://github.com/lfy79001/S3HQA.

##  2 Our Approach

###  2.1 Problem Definition

Given a natural language question
$\mathcal{Q}=\left\\{q_{i}\right\\}_{i=1}^{|\mathcal{Q}|}$ and a table
$\mathcal{T}$ with $\left\langle\mathcal{H},\mathcal{R}\right\rangle$,
$\mathcal{H}$ indicates the table headers, and
$\mathcal{R}=\left\\{r_{i}\right\\}_{i=1}^{|\mathcal{R}|}$ indicates the rows
with number $|\mathcal{R}|$. Each row $r_{i}$ is consists of $N$ cells
$r_{i}=\left\\{c_{ij}\right\\}_{j=1}^{N}$. The header’s number is also $N$.
Some cells have a linked passage $\mathcal{P}_{ij}$. Our goal aims to generate
the answer $\mathcal{A}$ with model $\Theta$, which is a span from table cells
or linked passage or a derivation result of counting questions.

###  2.2 Retriever with Refinement Training

The retriever aims to perform initial filtering of heterogeneous resources.
However, accurately labeling the location of answers consumes high labeling
costs. For TextTableQA data, the answer $\mathcal{A}$ usually appears in
multiple locations, which makes it difficult for us to generate precise
retrieval labels. We use a two-step training method, with a row-based
retriever and a passage-based retriever for each step.

Inspired by Kumar et al. (2021), the retrieval has two steps. First, we divide
the data $\mathcal{D}$ into two folds according to the string matching labels
$G_{i}$. Specifically, for a question-answer instance, the answer
$\mathcal{A}$ appears one time as $\mathcal{D}_{1}$, and the instance whose
answer $\mathcal{A}$ appears multiple times as $\mathcal{D}_{2}$. Take the
example in Figure 1, Q1, Q3 belongs to $\mathcal{D}_{1}$ while Q2 belongs to
$\mathcal{D}_{2}$. The data is organized in the form of
$\mathrm{[CLS]}q_{1}q_{2}...q_{|Q|}\mathrm{[SEP]c_{i1}c_{i2}...c_{iN}[SEP]}$
or $\mathrm{[CLS]}q_{1}q_{2}...q_{|Q|}\mathrm{[SEP]p_{ij}[SEP]}$.

In the first step, we only use $\mathcal{D}_{1}$ to train a model
$\Theta_{1}$, which data are noiseless. Then in the second step, we use the
trained weight $\Theta_{1}$ to train the model $\Theta_{2}$. For the input
$x$, the loss function is:

$L(\Theta_{2},x,\mathcal{R})=\sum_{z\in\mathcal{R}}-q(z)\log{p_{\Theta_{1}}(z|x)}$

where $q(z)=p_{\Theta_{1}}(z|x,z\in\mathcal{R})$ is the probability
distribution given by the model restricted to candidate rows $\mathcal{R}$
containing the answer span, taken here as a constant with zero gradients
(Eisenschlos et al., 2021).

Meanwhile, we use a passage-based retriever to enhance the performance of a
row-based retriever (PassageFilter). Specifically, we use the passage-based
retriever to obtain a prediction score of passage relevance. Based on this
score, we reorder the input of the row-based retriever. It avoids the
limitation on input sequence length imposed by the pre-trained model.

###  2.3 Hybrid Selector

This module needs to combine the results of the two granularity retrievers. As
for this task, we consider the question type and the relationships between the
table and linked passages essential. As shown in Figure 2, the hybrid selector
chooses the appropriate data source from the two retrieval results depending
on question types.

Specifically, for general _bridge_ multi-hop questions, we use a single row
and its linked passage. While for _comparison/count_ questions, we consider
multiple rows and further filter the related sentences, delete the linked
paragraphs with the low scores. This not only enables the generation module to
obtain accurate information, but also prevents the introduction of a large
amount of unrelated information. The selector algorithm outputs a mixed
sequence with high relevance based on the relationship between the question,
the table, and the passages. The algorithm is shown in Algorithm 1.

Algorithm 1 Hybrid Selector Algorithm.

0: question $\mathcal{Q}$, table rows $\mathcal{R}$, linked passages
$\mathcal{P}$, row-based retriever $\Theta_{R}$, passage-based retriever
$\Theta_{P}$, selector target row count $N_{S}$

0: generator input $\mathcal{S}$_Get the row/passage ordered list by relevant
scores_

1: $\mathcal{O_{R}}\leftarrow sort(\Theta_{R}(\mathcal{Q},\mathcal{R}))$

2: $\mathcal{O_{P}}\leftarrow sort(\Theta_{P}(\mathcal{Q},\mathcal{P}))$

3: $\mathrm{p}^{\mathrm{type}}\leftarrow Classification(Q)$

4: if $\mathrm{p}^{\mathrm{type}}=bridge$ then

5: if $\mathcal{O_{P}}[0]$ in $\mathcal{O_{R}}[0]$ then

6: $\mathcal{S}\leftarrow\mathcal{Q}+\mathcal{O_{R}}[0]$

7: else

8: $\mathcal{S}\leftarrow\mathcal{Q}+\mathcal{O_{R}}[0]+\mathcal{O_{P}}[0]$

9: end if

10: else

11: $\mathcal{O_{PC}}\leftarrow\mathcal{P}[len(\mathcal{O_{P}})//2:]$

12:
$\mathcal{S}\leftarrow\mathcal{Q}+\mathcal{O_{R}}[0:N_{S}]-\mathcal{O_{PC}}$

13: end if

14: return $\mathcal{S}$

###  2.4 Generation-based Reasoner

The results of the selector take into account both two granularity. Unlike the
previous approaches, which were based on a span extraction module, we use a
generation-based model for answer prediction.

####  2.4.1 Row-wise generator

To generate an accurate answer string $\mathcal{A}=(a_{1},a_{2},...,a_{n})$
given the question $\mathcal{Q}$ and selection evidence $\mathcal{S}$, we
perform lexical analysis to identify the question type, such as counting or
comparison, by looking for certain keywords or comparative adjectives. We
utilize two special tags $\mathrm{\left\langle Count\right\rangle}$ and
$\mathrm{\left\langle Compare\right\rangle}$, which indicates the question
types.

We then use the results of the passage retriever to rank the passages in order
of their relevance, eliminating the impact of model input length limitations.
Finally, we train a Seq2Seq language model with parameters $\Theta$, using the
input sequence $\mathcal{Q},\mathcal{S}$ and the previous outputs $a_{<i}$ to
optimize the product of the probabilities of the output sequence
$a_{1},a_{2},...,a_{n}$:

$\mathcal{A}=argmax\prod_{i=1}^{n}P(a_{i}|a_{<i},\mathcal{Q},\mathcal{S};\Theta)$

| Table | Passage | Total
---|---|---|---
| Dev | Test | Dev | Test | Dev | Test
| EM | F1 | EM | F1 | EM | F1 | EM | F1 | EM | F1 | EM | F1
Unsupervised-QG (Pan et al., 2021) | - | - | - | - | - | - | - | - | 25.7 | 30.5 | - | -
HYBRIDER (Chen et al., 2020b) | 54.3 | 61.4 | 56.2 | 63.3 | 39.1 | 45.7 | 37.5 | 44.4 | 44.0 | 50.7 | 43.8 | 50.6
DocHopper (Sun et al., 2021) | - | - | - | - | - | - | - | - | 47.7 | 55.0 | 46.3 | 53.3
MuGER2 (Wang et al., 2022b) | 60.9 | 69.2 | 58.7 | 66.6 | 56.9 | 68.9 | 57.1 | 68.6 | 57.1 | 67.3 | 56.3 | 66.2
POINTR (Eisenschlos et al., 2021) | 68.6 | 74.2 | 66.9 | 72.3 | 62.8 | 71.9 | 62.8 | 71.9 | 63.4 | 71.0 | 62.8 | 70.2
DEHG (Feng et al., 2022) | - | - | - | - | - | - | - | - | 65.2 | 76.3 | 63.9 | 75.5
MITQA (Kumar et al., 2021) | 68.1 | 73.3 | 68.5 | 74.4 | 66.7 | 75.6 | 64.3 | 73.3 | 65.5 | 72.7 | 64.3 | 71.9
MAFiD (Lee et al., 2023) | 69.4 | 75.2 | 68.5 | 74.9 | 66.5 | 75.5 | 65.7 | 75.3 | 66.2 | 74.1 | 65.4 | 73.6
S3HQA | 70.3 | 75.3 | 70.6 | 76.3 | 69.9 | 78.2 | 68.7 | 77.8 | 68.4 | 75.3 | 67.9 | 75.5
Human | - | - | - | - | - | - | - | - | - | - | 88.2 | 93.5
Table 1: Performance of our model and related work on the HybridQA dataset.

####  2.4.2 LLM prompting generator

With the emergence of large language models, In-Context Learning Dong et al.
(2022) and Chain-of-Thought prompting Wei et al. (2022) have become two
particularly popular research topics in this field. In this paper, we
introduce a prompting strategy for multi-hop TextTableQA.

We utilize selection evidence $\mathcal{S}$ and apply LLM-based prompting. We
conducted experiments on both vanilla prompting and chain-of-thought prompting
in zero-shot and few-shot scenarios.

##  3 Experiment

###  3.1 Experiment Setup

Datasets We conduct experiments on HybridQA (Chen et al., 2020b). The detailed
statistics are shown in Appendix A. For evaluation, we followed the official
evaluation to report exact match accuracy and F1 score.
Implementation details The implementation details are shown in Appendix B. The
experimental results are the average of five times results.

###  3.2 Fully-supervised Results

Table 1 shows the comparison results between our models with previous typical
approaches on both development and test sets. It shows that our proposed S3HQA
works significantly better than the baselines in terms of EM and F1 on
HybridQA. The results indicate that S3HQA is an effective model for multi-hop
question answering over tabular and textual data. Specifically, it can
effectively handle multi-hop reasoning and make full use of heterogeneous
information.

However, we found that our approach was outperformed by the DEHG model (Feng
et al., 2022) in terms of F1 score on the Dev set. We speculate that this
might be because the DEHG approach uses their own Open Information Extraction
(OIE) tool.

Model | Dev
---|---
EM | F1
Zero-shot prompt
GPT3.5 direct | 33.1 | 50.5
GPT3.5 CoT | 52.9 | 66.6
Few-shot prompt (2-shot)
GPT3.5 direct | 57.1 | 68.8
GPT3.5 CoT | 60.3 | 72.1
Table 2: Performance Comparison of LLM-Prompting Method on Zero-Shot and Few-
Shot Scenarios for HybridQA Dataset.

###  3.3 LLM-prompting Results

We present our zero-shot and few-shot results in Table 2. "Direct" refers to a
simple prompting method where only the question, context, and answer are
provided to the model without any additional reasoning process. In contrast,
"CoT" involves a human-authored Chain-of-Thought reasoning process that
provides a more structured and logical way of prompting the model. The
experiments demonstrate that in-context learning used to prompt large language
models can achieve promising results. Specifically, utilizing the Chain-of-
Thought prompt method can significantly enhance the model’s performance.

However, it’s worth noting that there is still a performance gap compared to
fine-tuning the model on the full dataset (Table 1). Fine-tuning allows the
model to learn more specific information about the TextTableQA task, resulting
in better performance. Nevertheless, our results show that the LLM-prompting
method can be a useful alternative to fine-tuning, especially when there is a
limited amount of labeled data available.

###  3.4 Ablation Studies

We conduct ablation studies on the test set. We validate the effects of three
modules: _retriever with refinement training_ , _hybrid selector_ , and
_generation-based reasoner_. The retriever performs initial filtering of
heterogeneous resources; Selectors combined with hyperlinks further identify
the exact evidence needed to answer multi-hop questions; and the reasoner uses
the selection evidence to obtain the final answer.

Model | Top1
---|---
$\textrm{S}^{3}\textrm{HQA-Retriever}_{\textrm{DB}}$ | 88.0
$\textrm{S}^{3}\textrm{HQA-Retriever}_{\textrm{BE}}$ | 87.3
w/o Refinement training | 84.1
w/o PassageFilter | 85.3
$\textrm{Vanilla-Retriever}_{\textrm{BE}}$ | 82.0
Table 3: Ablation study of retrieval results. DB and BE denote models based on Deberta-base (He et al., 2020) and BERT-base-uncased (Devlin et al., 2018), respectively Model | EM | F1
---|---|---
S3HQA | 67.9 | 76.5
w/o hybrid selector | 65.0 | 74.9
w/o special tags | 67.2 | 76.0
BERT-large reader | 66.8 | 75.8
Table 4: Ablation study of S3HQA.

Effect of proposed retriever. As shown in the Table 3, under the setting of
using the BERT-base-uncased model, sing the BERT-base-uncased model setting,
the retriever with _refinement training_ achieved 87.2. When we use Deberta-
base, the top1 retrieval performance improved by 0.8$\%$. For _w/o refinement
training_ , we use the entire data directly for training, the top1 recall
drops about 3.2$\%$. For _w/o PassageFilter_ , we remove the mechanism, the
top1 recall drops about 3.2$\%$. For _Vanilla-Retriever_ , we use the row-
based retriever (Kumar et al., 2021) and remove all our mechanisms, the top1
score drops about 5.3$\%$. This shows that our model can solve the weakly
supervised data noise problem well.

Effect of hybrid selector. As shown in the Table 4, we removed the selector of
S3HQA and replaced it with the previous cell-based selector (Wang et al.,
2022b). This method directly uses the top1 result of the row retriever as
input to the generator. _w/o hybrid selector_ shows that the EM drops 2.9$\%$
and F1 drops 1.6$\%$, which proves the effectiveness of our selector approach.

Effect of reasoner. As shown in the Table 4, we design two baselines. _BERT-
large reader_ (Chen et al., 2020b; Wang et al., 2022b) uses BERT (Devlin et
al., 2018) as encoder and solves this task by predicting the start/end tokens.
_w/o special tags_ deletes the special tags. Both the two experiments
demonstrate our S3HQA reasoner performs the best for HybridQA task.

##  4 Related Work

The TextTableQA task (Wang et al., 2022a) has attracted more and more
attention. As for multi-hop type dataset, previous work used pipeline approach
(Chen et al., 2020b), unsupervised approach (Pan et al., 2021), multi-
granularity (Wang et al., 2022b), table pre-trained language model
(Eisenschlos et al., 2021), multi-instance learning (Kumar et al., 2021) and
graph neural network (Feng et al., 2022) to solve this task. As for numerical
reasoning task, which is quite different from multi-hop type dataset, there is
also a lot of work (Zhu et al., 2021; Zhao et al., 2022; Zhou et al., 2022;
Lei et al., 2022; Li et al., 2022; Wei et al., 2023) to look at these types of
questions. Unlike these methods, our proposed three-stage model S3HQA can
alleviate noises from weakly supervised and solve different types of multi-hop
TextTableQA questions by handling the relationship between tables and text.

##  5 Conclusion

This paper proposes a three-stage model consisting of retriever, selector, and
reasoner, which can effectively address multi-hop TextTableQA. The proposed
method solves three drawbacks of the previous methods: noisy labeling for
training retriever, insufficient utilization of heterogeneous information, and
deficient ability for reasoning. It achieves new state-of-the-art performance
on the widely used benchmark HybridQA. In future work, we will design more
interpretable TextTableQA models to predict the explicit reasoning path.

## Limitations

Since the multi-hop TextTableQA problem has only one dataset HybridQA, our
model has experimented on only one dataset. This may lead to a lack of
generalizability of our model. Transparency and interpretability are important
in multi-hop question answering. While our model achieves the best results,
the model does not fully predict the reasoning path explicitly and can only
predict the row-level path and passage-level path. In future work, we will
design more interpretable TextTableQA models.

## Acknowledgements

This work was supported by the National Key R&D Program of China
(2022ZD0160503) and the National Natural Science Foundation of China
(No.U1936207, No.61976211). This work was supported by the Strategic Priority
Research Program of Chinese Academy of Sciences (No.XDA27020100), the Youth
Innovation Promotion Association CAS, Yunnan Provincial Major Science and
Technology Special Plan Projects (No.202202AD080004) and CCF-DiDi GAIA
Collaborative Research Funds for Young Scholars.

## References

  * Bird (2006) Steven Bird. 2006.  Nltk: the natural language toolkit.  In _Proceedings of the COLING/ACL 2006 Interactive Presentation Sessions_ , pages 69–72.
  * Chen et al. (2020a) Wenhu Chen, Ming-Wei Chang, Eva Schlinger, William Yang Wang, and William W Cohen. 2020a.  Open question answering over tables and text.  In _International Conference on Learning Representations_.
  * Chen et al. (2020b) Wenhu Chen, Hanwen Zha, Zhiyu Chen, Wenhan Xiong, Hong Wang, and William Yang Wang. 2020b.  Hybridqa: A dataset of multi-hop question answering over tabular and textual data.  In _Findings of the Association for Computational Linguistics: EMNLP 2020_ , pages 1026–1036.
  * Chen et al. (2021) Zhiyu Chen, Wenhu Chen, Charese Smiley, Sameena Shah, Iana Borova, Dylan Langdon, Reema Moussa, Matt Beane, Ting-Hao Huang, Bryan R Routledge, et al. 2021\.  Finqa: A dataset of numerical reasoning over financial data.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 3697–3711.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_.
  * Dong et al. (2022) Qingxiu Dong, Lei Li, Damai Dai, Ce Zheng, Zhiyong Wu, Baobao Chang, Xu Sun, Jingjing Xu, and Zhifang Sui. 2022.  A survey for in-context learning.  _arXiv preprint arXiv:2301.00234_.
  * Eisenschlos et al. (2021) Julian Eisenschlos, Maharshi Gor, Thomas Mueller, and William Cohen. 2021.  Mate: Multi-view attention for table transformer efficiency.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 7606–7619.
  * Feng et al. (2022) Yue Feng, Zhen Han, Mingming Sun, and Ping Li. 2022.  Multi-hop open-domain question answering over structured and unstructured knowledge.  In _Findings of the Association for Computational Linguistics: NAACL 2022_ , pages 151–156.
  * He et al. (2020) Pengcheng He, Xiaodong Liu, Jianfeng Gao, and Weizhu Chen. 2020.  Deberta: Decoding-enhanced bert with disentangled attention.  _arXiv preprint arXiv:2006.03654_.
  * Kumar et al. (2021) Vishwajeet Kumar, Saneem Chemmengath, Yash Gupta, Jaydeep Sen, Samarth Bharadwaj, and Soumen Chakrabarti. 2021.  Multi-instance training for question answering across table and linked text.  _arXiv preprint arXiv:2112.07337_.
  * Lee et al. (2023) Sung-Min Lee, Eunhwan Park, Daeryong Seo, Donghyeon Jeon, Inho Kang, and Seung-Hoon Na. 2023.  Mafid: Moving average equipped fusion-in-decoder for question answering over tabular and textual data.  In _Findings of the Association for Computational Linguistics: EACL 2023_ , pages 2292–2299.
  * Lei et al. (2022) Fangyu Lei, Shizhu He, Xiang Li, Jun Zhao, and Kang Liu. 2022.  Answering numerical reasoning questions in table-text hybrid contents with graph-based encoder and tree-based decoder.  In _Proceedings of the 29th International Conference on Computational Linguistics_ , pages 1379–1390.
  * Lewis et al. (2020) Mike Lewis, Yinhan Liu, Naman Goyal, Marjan Ghazvininejad, Abdelrahman Mohamed, Omer Levy, Veselin Stoyanov, and Luke Zettlemoyer. 2020.  Bart: Denoising sequence-to-sequence pre-training for natural language generation, translation, and comprehension.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 7871–7880.
  * Li et al. (2022) Xiao Li, Yin Zhu, Sichen Liu, Jiangzhou Ju, Yuzhong Qu, and Gong Cheng. 2022.  Dyrren: A dynamic retriever-reranker-generator model for numerical reasoning over tabular and textual data.  _arXiv preprint arXiv:2211.12668_.
  * Pan et al. (2021) Liangming Pan, Wenhu Chen, Wenhan Xiong, Min-Yen Kan, and William Yang Wang. 2021\.  Unsupervised multi-hop question answering by question generation.  In _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 5866–5880.
  * Pasupat and Liang (2015) Panupong Pasupat and Percy Liang. 2015.  Compositional semantic parsing on semi-structured tables.  In _Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 1470–1480.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. 2019\.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32.
  * Rajpurkar et al. (2016) Pranav Rajpurkar, Jian Zhang, Konstantin Lopyrev, and Percy Liang. 2016.  Squad: 100,000+ questions for machine comprehension of text.  In _Proceedings of the 2016 Conference on Empirical Methods in Natural Language Processing_ , pages 2383–2392.
  * Sun et al. (2021) Haitian Sun, William W Cohen, and Ruslan Salakhutdinov. 2021.  End-to-end multihop retrieval for compositional question answering over long documents.
  * Wang et al. (2022a) Dingzirui Wang, Longxu Dou, and Wanxiang Che. 2022a.  A survey on table-and-text hybridqa: Concepts, methods, challenges and future directions.  _arXiv preprint arXiv:2212.13465_.
  * Wang et al. (2022b) Yingyao Wang, Junwei Bao, Chaoqun Duan, Youzheng Wu, Xiaodong He, and Tiejun Zhao. 2022b.  MuGER2: Multi-granularity evidence retrieval and reasoning for hybrid question answering.  In _Findings of the Association for Computational Linguistics: EMNLP 2022_ , pages 6687–6697, Abu Dhabi, United Arab Emirates. Association for Computational Linguistics.
  * Wei et al. (2022) Jason Wei, Xuezhi Wang, Dale Schuurmans, Maarten Bosma, Ed Chi, Quoc Le, and Denny Zhou. 2022.  Chain of thought prompting elicits reasoning in large language models.  _arXiv preprint arXiv:2201.11903_.
  * Wei et al. (2023) Yifan Wei, Fangyu Lei, Yuanzhe Zhang, Jun Zhao, and Kang Liu. 2023.  Multi-view graph representation learning for answering hybrid numerical reasoning question.  _arXiv preprint arXiv:2305.03458_.
  * Yu et al. (2018) Tao Yu, Rui Zhang, Kai Yang, Michihiro Yasunaga, Dongxu Wang, Zifan Li, James Ma, Irene Li, Qingning Yao, Shanelle Roman, et al. 2018.  Spider: A large-scale human-labeled dataset for complex and cross-domain semantic parsing and text-to-sql task.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pages 3911–3921.
  * Zhao et al. (2022) Yilun Zhao, Yunxiang Li, Chenying Li, and Rui Zhang. 2022.  Multihiertt: Numerical reasoning over multi hierarchical tabular and textual data.  In _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 6588–6600.
  * Zhou et al. (2022) Yongwei Zhou, Junwei Bao, Chaoqun Duan, Youzheng Wu, Xiaodong He, and Tiejun Zhao. 2022.  Unirpg: Unified discrete reasoning over table and text as program generation.  _arXiv preprint arXiv:2210.08249_.
  * Zhu et al. (2021) Fengbin Zhu, Wenqiang Lei, Youcheng Huang, Chao Wang, Shuo Zhang, Jiancheng Lv, Fuli Feng, and Tat-Seng Chua. 2021.  Tat-qa: A question answering benchmark on a hybrid of tabular and textual content in finance.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 3277–3287.

##  Appendix A HybridQA Dataset

HybridQA is a large-scale, complex, and multi-hop TextTableQA benchmark.
Tables and texts are crawled from Wikipedia. Each row in the table describes
several attributes of an instance. Each table has its hyperlinked Wikipedia
passages that describe the detail of attributes. It contains 62,682 instances
in the train set, 3466 instances in the dev set and 3463 instances in the test
set.

Split | Train | Dev | Test | Total
---|---|---|---|---
In-Passage | 35,215 | 2,025 | 20,45 | 39,285 (56.4%)
In-Table | 26,803 | 1,349 | 1,346 | 29,498 (42.3%)
Computed | 664 | 92 | 72 | 828 (1.1%)
Total | 62,682 | 3,466 | 3,463 | 69,611
Table 5: Data Split: In-Table means the answer comes from plain text in the
table, and In-Passage means the answer comes from certain passage.

##  Appendix B Implementation Details

###  B.1 Fully-supervised Setting

We utilize PyTorch (Paszke et al., 2019) to implement our proposed model.
During pre-processing, the input of questions, tables and passages are
tokenized and lemmatized with the NLTK (Bird, 2006) toolkit. We conducted the
experiments on a single NVIDIA GeForce RTX 3090.

In the retriever stage, we use BERT-base-uncased (Devlin et al., 2018) and
Deberta-base (He et al., 2020) to obtain the initial representations. For the
first step, batch size is 1, epoch number is 5, learning rate is 7e-6
(selected from 1e-5, 7e-6, 5e-6). The training process may take around 10
hours. For the second step, we use a smaller learning rate 2e-6 (selected from
5e-6, 3e-6, 2e-6), epoch number is 5. The training process may take around 8
hours. In the selector stage, target row count $N_{S}$ is 3. In the generator
stage, we use BART-large language model (Lewis et al., 2020), the learning
rate is 1e-5 (selected from 5e-5, 1e-5, 5e-6), batch size is 8, epoch number
is 10, beam size is 3 and max generate length is 20.

###  B.2 LLM-prompting Setting

We use the OpenAI GPT-3.5 (text-davinci-003) API model with the setting
$temperature=0$ in our experiments. For the few-shot setting, we use 2 shots.
To elicit the LLM’s capability to perform multi-hop reasoning, we use the text
"Read the following table and text information, answer a question. Let’s think
step by step." as our prompt.


# Near-monochromatic tuneable cryogenic niobium electron field emitter

C.W. Johnson1 A.K. Schmid1 M. Mankos2 R. Röpke3 N. Kerker3 E.K. Wong1 D.F.
Ogletree1 A.M. Minor1,4 A. Stibor1,3∗ 1Lawrence Berkeley National Lab,
Molecular Foundry, Berkeley, USA  2Electron Optica, Inc., Palo Alto, USA
3Institute of Physics and LISA+, University of Tübingen, Tübingen, Germany,
4Department of Materials Science and Engineering, University of California,
Berkeley, USA<EMAIL_ADDRESS>
###### Abstract

Creating, manipulating, and detecting coherent electrons is at the heart of
future quantum microscopy and spectroscopy technologies. Leveraging and
specifically altering the quantum features of an electron beam source at low
temperatures can enhance its emission properties. Here, we describe electron
field emission from a monocrystalline, superconducting niobium nanotip at a
temperature of 5.9 K. The emitted electron energy spectrum reveals an ultra-
narrow distribution down to 16 meV due to tunable resonant tunneling field
emission via localized band states at a nano-protrusion’s apex and a cut-off
at the sharp low-temperature Fermi-edge. This is an order of magnitude lower
than for conventional field emission electron sources. The self-focusing
geometry of the tip leads to emission in an angle of 3.2∘, a reduced
brightness of 5.0 $\times$ 108 A/(m2 sr V) and a stability of several hours.
This source will decrease the impact of lens aberration and enable new modes
in low energy electron microscopy, electron energy loss spectroscopy, and
high-resolution vibrational spectroscopy.

With the rise of novel instruments in electron microscopy Kruit  _et al._
(2016); Turner  _et al._ (2021), spectroscopy Küster  _et al._ (2021);
Krivanek  _et al._ (2014), quantum electronics and electron wave function
engineering Madan  _et al._ (2020), the demand for coherent beam sources
becomes increasingly relevant in various fields of physics. Significant
progress has been made in the generation of novel electron field emitters,
that are e.g. laser-triggered Ehberger  _et al._ (2015), have a single atom
apex Kuo  _et al._ (2006); Pooch  _et al._ (2018) or consist of a LaB6
nanowire Zhang  _et al._ (2022). The key performance goals are high brightness
combined with a narrow field emission electron energy distribution, resulting
in large longitudinal and transversal coherence. Such sources allow new
techniques for electron phase modulation Johnson  _et al._ (2020); Feist  _et
al._ (2020), beam guiding Hammer  _et al._ (2015), interferometry Schütz  _et
al._ (2014); Johnson  _et al._ (2021), quantum decoherence measurement Kerker
_et al._ (2020); Beierle  _et al._ (2018), novel microscopy diffraction modes
Latychevskaia  _et al._ (2017), quantum sensing Pooch  _et al._ (2017), or
even secure information transfer Röpke  _et al._ (2021). They lead to key
applications in superconducting scanning tunneling microscopy Huang  _et al._
(2020), quantum electron microscopy Kruit  _et al._ (2016); Turner  _et al._
(2021) or scanning Josephson spectroscopy Küster  _et al._ (2021). According
to the Fowler-Nordheim (FN) description Forbes  _et al._ (2013), the energy
distribution of a beam gets more narrow with lower field emitter temperature.
But the effect is negligible below room temperature and usually reported
energy widths of cold field emitter are
$\sim$$200\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ Peter W. Hawkes and Erwin
Kasper (2017).

In this Letter, we analyze the field emission properties of a monocrystalline
niobium (Nb) nanotip at a temperature of $5.9\text{\,}\mathrm{K}$, well below
the superconducting transition at $T_{c}=$9.3\text{\,}\mathrm{K}$$. A scanning
electron microscope image of the tip can be seen in Fig. 1 a-c). It is
demonstrated that a surface nano-protrusion (NP) can be generated on the tip,
leading to a localized quantum band state at the apex Binh  _et al._ (1992);
Gohda and Watanabe (2001). The nano-protrusion emission (NPE) is self-focusing
with an emission angle of $3.2\text{\,}\mathrm{\SIUnitSymbolDegree}$ due to
the high field enhancement at the NP apex with a small radius of curvature.
The center of the NPE total energy distribution spectral peak can be tuned
relative to the Fermi-energy ($E_{F}$) and cut-off by the sharp low-
temperature Fermi-edge by changing the electric field at the emitter. It leads
to an extremely narrow energy distribution full-width-half-maximum (FWHM) down
to $\Delta E=$16\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$. This is the first
monocrystalline Nb tip characterized at superconducting temperatures and the
smallest reported electron field emitter energy width. It is 10-20 times
narrower than the best emitters available in microscopy, leading to high
coherence, while showing long-time stability and high brightness. This immense
increase in performance of our beam source will be a game-changer in electron
microscopy and spectroscopy.

Figure 1: a-c) SEM image of the FIB milled monocrystalline nanotip in three
different magnifications. d) On top of the R${}_{0}=$
$25\text{\,}\mathrm{n}\mathrm{m}$ radius Nb tip, an r${}_{0}=$
$1.3\text{\,}\mathrm{n}\mathrm{m}$-nano-protrusion is formed that is not
resolved by SEM. The nano-protrusion’s field emission is magnified by an
einzel lens. e) Delay line detector image of the beam’s angular distribution
with $V_{tip}=$ $-364\text{\,}\mathrm{V}$ ($V_{ext}=$ $0\text{\,}\mathrm{V}$),
yielding a mean angular spread of $\Delta\theta\approx$
$3.2\text{\,}\mathrm{\SIUnitSymbolDegree}$. f) Sketch of our electron field
emitter characterization setup. The tip and the bar are sandwiched between two
sapphire plates (visible in b)) and mounted onto a rotatable cryostat
perpendicular to the rotation axis, allowing it to be pointed at multiple
measurement devices. g) Four-point resistance vs. temperature measurement of
the tip’s bar, indicating superconducting emitter conditions.

The Nb nanotip and the NP on its apex were fabricated in five steps. First, a
rectangular monocrystalline wire being the base of the tip was cut out from a
larger [100] oriented Nb crystal in pieces of
250$\times$$250\text{\,}\mathrm{\SIUnitSymbolMicro m}$ cross-section with
electric discharge machining sur . Then this wire was spot-welded in the
center of a V-shaped polycrystalline Nb wire bar. As a next step, the
monocrystalline tip base was electrochemically etched down to a length of only
$\sim$$2\text{\,}\mathrm{m}\mathrm{m}$ in a KOH solution according to a
procedure described in Hou  _et al._ (2016). This led to the
$\sim$$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$-base shown in the SEM images
of Fig. 1 a-b). The tip is then further shaped by gallium ion beam milling in
a FIB Hel  to fabricate the conical nanotip in Fig. 1 c) with a radius
$R_{0}=$ $25\text{\,}\mathrm{nm}$. The most delicate step is the formation of
a NP ($<5\text{\,}\mathrm{n}\mathrm{m}$ Binh  _et al._ (1992)) on the tip
after being installed and cooled on the cryostat. This is performed through
several annealing cycles by ramping up a current of $4.85\text{\,}\mathrm{A}$
through the polycrystalline Nb wire bar until the tip is glowing at $\sim$950
∘C Lee  and simultaneously setting a $-3\text{\,}\mathrm{k}\mathrm{V}$
electrical bias on the extractor. By automating this annealing process with a
programmable current source, NP formation on our tip is reproducible and the
NPE is stable on the order of hours. However, there is still variability in
the emission properties of the NPs after each annealing process, making
several annealings necessary to get the optimal conditions.

The setup to cool the tip and to characterize its field emission properties is
illustrated in Fig. 1 f). It allows a high-resolution measurement of the beam
energy spectra, the emission current, and the angular distribution. The
nanotip and the extractor are cooled by a closed-cycle liquid helium cryostat
Adv  and thermally isolated by a copper cooling shield with a small opening
for the beam path. The step in the four-point resistance measurement of the
V-shaped Nb-bar in Fig. 1 g) is used to calibrate the measured temperatures to
$T_{c}$ and reveals that the emitter is operating under superconducting
conditions at $5.9\text{\,}\mathrm{K}$. The cryostat is mounted inside an
ultra-high vacuum chamber with a pressure of
$\sim$$8\text{\times}{10}^{-11}\text{\,}\mathrm{T}\mathrm{o}\mathrm{r}\mathrm{r}$
on a rotational flange and a 3D-manipulator. It allows positioning and
aligning of the emitter towards three measurement units. The first one is a
hemispherical electron beam energy analyzer Sci  with a resolution of
$3\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$. Between the analyzer and the tip,
a double deflector is installed for beam alignment. A rotation of the cryostat
points the tip towards a Faraday-cup Kim  with a pico-ammeter Kei  for the
beam current measurements. Another rotation points the emitter towards a
microchannel plate delay line detector (DLD Roe ) with a single-electron
resolution to measure the beam’s angular distribution. A custom einzel lens
with a deflector is also installed at this position for magnification and
alignment. The analyzer, Faraday cup, and DLD operate in different regimes of
current and tip voltage, allowing a wide range of beam analysis. This realizes
a unique system that combines spectroscopy, beam current, spatial
distribution, and temperature control.

a)b)c)d)e)f) Figure 2:  a) Qualitative illustration showing NP’s resonant
states in the Coulomb barrier, allowing electrons an enhanced field emission
to the vacuum. Lower panel: Decreasing the field strength in resonant
tunneling emission can shift resonances originally above EF into the occupied
region, resulting in a tuneable sharp cut-off at low temperatures and an
ultra-narrow beam energy distribution. b) Experimental spectra (normalized to
the largest peak) showing the nano-protrusion (NP) emission peaks decreasing
in energy by increasing the extractor voltage $V_{ext}$ from 0 to
$180\text{\,}\mathrm{V}$ at a fixed tip voltage of $V_{tip}=$
$-475\text{\,}\mathrm{V}$. c) The narrowest energy distribution observed of
$16\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ after noise correction
($19\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ before). The small energy width
is due to the use of the sharp low-temperature Fermi edge to increasingly
truncate the higher energy portion of the NP’s discrete state emission by
applying an extraction voltage of 180, 120 and $20\text{\,}\mathrm{V}$,
pushing the resonance peak over EF. d) Energy spectra corresponding to the
color-matching encircled data points in e) with the tip voltage labeled on the
right side ($V_{ext}=0$). The spectra with increasing Vtip correspond to the
circles in e) from right to left. e) Fowler-Nordheim (FN) plots of a NP’s
emission current with decreasing $V_{tip}$ from $-600\text{\,}\mathrm{V}$ to
$-420\text{\,}\mathrm{V}$. The orange line is a theory expectation for a pure
FN behavior, confirming that the resonant emission process is different
compared to a FN emitter. f) Stability measurement of the Nb field emitter
over $6\text{\,}\mathrm{h}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{s}$. The
sudden increase after $1\text{\,}\mathrm{h}\mathrm{o}\mathrm{u}\mathrm{r}$ is
a familiar effect observed in other cold field emitters and can be explained
by the temporary presence of an adatom on the emission site. The inset shows
energy spectra taken at the beginning of the recording of f), after
$3\text{\,}\mathrm{h}$ and at the end. A reasonable peak broadening is
visible, presumably due to adatoms accumulating near the tips apex or small
structural changes that alter the boundary conditions for the NP resonant
emission.

The energy spectrum of FN cold field emission of electrons from a metal’s
Fermi sea to vacuum with respect to $E_{F}$ is given by the product of the
Fermi-Dirac distribution function and an exponential decay term $G(E)\propto
f(E)\exp(E/D),$ where $f(E)=(1+\exp(E/k_{B}T))^{-1}$ and $D=e\hbar
F/2\sqrt{2m\phi}$ being the energy scale parameter for field emission with the
work function $\phi$ and the field strength at the emission site
$F=|\mathbf{E}|$ Young (1959). For a tip, it can be determined by $F=\beta
V_{tip}$, where $\beta=(kR_{0})^{-1}$ with $R_{0}$ being the base radius and a
geometry factor typically taken as $k=5$ for nanotips. The FN emission current
density goes as $J\propto F^{2}\exp(-B\phi^{3/2}/F)$, where $B\approx$
$7\text{\,}\mathrm{V}\mathrm{/}\mathrm{e}\mathrm{V}^{3/2}\mathrm{n}\mathrm{m}$.
At cryogenic temperatures, $f$ is nearly a step function, causing the energy
spectrum width to scale as $\Delta E=\ln(2)D$ due to the decreasing tunnel
width through the Coulomb barrier with increasing $F$. The addition of a NP
with radius $r_{0}\ll R_{0}$ on the tip apex creates localized, discrete
electronic states that can act as intermediate levels for resonant tunneling
through the Coulomb barrier. The energy spectrum of NPE via these resonant
states can be approximated as
$G_{tot}(E)=G(E)\left(1+\sum_{n}R_{n}(E-E_{n}+\alpha F)\right),$ where the
resonant enhancement factors $R_{n}$ are single peaked functions (Lorentzian
distributions in analytic approximation Lin  _et al._ (1992)) corresponding to
each discrete state at energy $E_{n}$ that shifts linearly with the applied
field strength at the tip. The maxima of the resonant enhancement factors can
be $|R_{n}|_{max}\gg\exp(-(E_{n}-\alpha F)/D)$ for emission well below
$E_{F}$, this means that NPE with isolated bright resonant tunneling spectral
peaks can occur even before any FN field emission is measurable. Additionally,
at cryogenic temperatures, the width of the Fermi edge is much narrower than
the width of any $R_{n}$ peak, allowing NPE peaks to be ’turned off’ or
partially truncated by tuning over $E_{F}$. The resonant tunneling field
emission process is illustrated in Fig. 2 a). The electric field around a
biased NP on a nanotip autofocuses the NPE into an emission angle
$\sim$2$-6\text{\,}\mathrm{\SIUnitSymbolDegree}$ Saenz  _et al._ (1990), as
illustrated in Fig. 1 d). This is much smaller than emission from a nanotip
geometry without the NP which would be within
$\sim$20$-30\text{\,}\mathrm{\SIUnitSymbolDegree}$ Saenz  _et al._ (1990).
Using an einzel lens with a calibrated magnification and the DLD, we
determined the angular beam distribution in Fig. 1 e) as $\Delta\theta\approx$
$3.2\text{\,}\mathrm{\SIUnitSymbolDegree}$ (azimuthally averaged emission
angle). With a reasonable assumption that the NP is hyperboloidal, we estimate
the size of the emission radius on the NP from the measured values of $R_{0}$
and $\Delta\theta$, using the relation $r_{0}=R_{0}\phi\Delta\theta^{2}/2D$
Saenz  _et al._ (1990). Inserting in the appropriate values for our tip
$R_{0}=$ $25\text{\,}\mathrm{n}\mathrm{m}$, $\phi=$
$4.5\text{\,}\mathrm{e}\mathrm{V}$, $\Delta\theta=$
$3.2\text{\,}\mathrm{\SIUnitSymbolDegree}$, $D=$
$134\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ at $V_{tip}=$
$-364\text{\,}\mathrm{V}$, we find an emission site radius of $r_{0}=$
$1.3\text{\,}\mathrm{n}\mathrm{m}$.

Furthermore, we recorded the energy spectrum of the NPE with an energy
analyzer. By keeping the tip voltage constant and varying the extractor
voltage to tune the NPE with respect to $E_{F}$, we recorded the NPE energy
spectra in Fig. 2 b) from two distinct NP states. The broad lower energy peak
is shifted away from $E_{F}$ at a linear rate of $\sim$2
$\mathrm{meV}\text{\,}{\mathrm{V}}^{-1}$. The higher energy peak is initially
above $E_{F}$ and begins to cross $E_{F}$ at $V_{ext}=$
$90\text{\,}\mathrm{V}$. It has a $\Delta E=$
$34\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ FWHM at $V_{ext}=$
$150\text{\,}\mathrm{V}$ and completely passes $E_{F}$ at $V_{ext}=$
$180\text{\,}\mathrm{V}$ showing a $\Delta E=$
$100\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ FWHM main peak. Nominally similar
values of peak widths and energy shifts are regularly reproduced after each
tip anneal and NP formation.

The direct tunability of these discrete NP states by the tip or extractor
voltage allows us to reduce the energy width by truncating a portion of the
NPE spectrum, with the sharp (less than
$1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$), low-temperature Fermi-Dirac
function at $E_{F}$. It is significantly narrower than the FWHM of the NPE
from the discrete states, producing a total FWHM $\Delta E=$
$19\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ at a tip voltage of $V_{tip}=$
$-475\text{\,}\mathrm{V}$ with $V_{ext}=$ $20\text{\,}\mathrm{V}$. This is the
smallest $\Delta E$ value ever recorded for field-emitted electrons. The
result is presented in Fig. 2 c) showing the spectrum width decreasing at
lower $V_{ext}$ as the peak shifts toward $E_{F}$ from below. We use
$G_{tot}(E)$ convolved with a kernel, being a normalized Gaussian distribution
$g(E,dE)=\exp(-x^{2}/2dE^{2})/\sqrt{2\pi dE^{2}}$ to account for systematic
instrument noise (analyzer, voltage source, vibrations), as a fit for the
total energy distribution: $G_{fit}(E)=G_{tot}(E)*g(E)$. We took
$R_{n}(E)=r_{n}/((E-E_{n})^{2}+\Gamma_{n}^{2})$ as the form for the resonant
enhancement factors, where $r_{n}$ are the resonance amplitudes, and
$\Gamma_{n}$ are the resonance linewidths. With $r_{n}$, $E_{n}$,
$\Gamma_{n}$, and $dE$ as fitting parameters, we performed a least-squares-fit
to the experimental spectra. Using the fit parameters $r_{n}$, $E_{n}$,
$\Gamma_{n}$ in the function $G_{tot}$, we determined an
$8.2\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ system energy resolution. If it
is deconvolved from the measured energy distribution, our narrowest emission
peak has a FWHM of $16\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$.

In the next step, a Faraday cup was used to measure the total NPE current.
After further annealing, a NP was formed and the tip voltage was decreased
from $V_{tip}=$ $-600\text{\,}\mathrm{V}$ to $V_{tip}=$
$-420\text{\,}\mathrm{V}$ in steps of $10\text{\,}\mathrm{V}$ at zero
extraction voltage, and the emission current was recorded. The $I/V$-data is
shown in Fig. 2 e) in a FN representation, where the orange linear relation
between $1/V_{tip}$ and $\log\left(I/V_{tip}^{2}\right)$ would be expected
applying the FN theory. Clearly, our measured resonance emission data
indicates a strong deviation to FN. After this measurement, the tip was
pointed towards the energy analyzer for recording of the spectra at the same
tip voltages. It allows correlating the current data with the energy widths of
the emissions, as shown in Fig. 2 d). The inset in Fig. 2 e) plots the energy
FWHM $\Delta E$ as a function of the tip voltage, again with an orange line
reflecting the FN expectation. It demonstrates that resonant field emission
can provide a significantly narrower energy width than the FN process. After
further annealing, we performed a field emission stability measurement over
$6\text{\,}\mathrm{h}$, as shown in Fig. 2 f). It was started by the recording
of a spectrum at $V_{tip}=-$$480\text{\,}\mathrm{V}$ with an energy FWHM of
$69\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ at an emission current of about
$4.1\text{\,}\mathrm{n}\mathrm{A}$. Subsequently, the current was recorded for
$3\text{\,}\mathrm{h}$ before the tip was pointed again towards the energy
analyzer, where a spectrum with a FWHM of
$95\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ was recorded. After further
current measurement for $3\text{\,}\mathrm{h}$, a final spectrum was taken
with a width of $155\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$. The spectra are
shown in the inset of Fig. 2 f). Finally, the angular beam distribution (at a
lower tip voltage to avoid DLD damage) presented in Fig. 1 e) was recorded.

One important figure of merit for electron sources is the reduced brightness
$B_{r}=J/\Delta\Omega V_{tip}$, where $J=I/\pi r_{0}^{2}$ is the current
density, and $\Delta\Omega=\pi\Delta\theta^{2}/4$ is the solid angle of
emission. This is an invariant quantity along the beam path and a more apt way
to compare different emission sources than current alone. In the NP model, the
angle of emission goes as $\Delta\theta\propto\sqrt{V_{tip}}$. Since we are
unable to measure the angular distribution of the NPE at high currents with
the microchannel plate of the DLD, we can use the relation
$\Delta\theta_{2}=\Delta\theta_{1}\sqrt{V_{tip,2}/V_{tip,1}}$ to determine the
high current, $4.1\text{\,}\mathrm{n}\mathrm{A}$, angular distribution
$\Delta\theta_{2}$ at $V_{tip,2}=-$$480\text{\,}\mathrm{V}$ from the
measurement in Fig. 1 e) of $\Delta\theta_{1}=$
$3.2\text{\,}\mathrm{\SIUnitSymbolDegree}$ at
$V_{tip,1}=-$$364\text{\,}\mathrm{V}$. Then, using the determined NPE radius
$r_{0}=$ $1.3\text{\,}\mathrm{n}\mathrm{m}$ we find a reduced brightness of
our source to be $B_{r}=$
$5.0\text{\times}{10}^{8}\text{\,}\mathrm{A}\mathrm{/}\mathrm{(}\mathrm{{}^{2}}\mathrm{m}\mathrm{sr}\mathrm{V}\mathrm{)}$
at an energy FWHM of $69\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$. This
extremely high brightness combined with the low energy width outperforms
commercial cold field emission sources Börrnert  _et al._ (2018); Shao  _et
al._ (2018) and other NPE tips Nagaoka  _et al._ (2001). For comparison, the
commercial source with a monochromator in Mook and Kruit (2000) has a
comparable energy distribution of $61\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$
but only a reduced brightness of
$6.3\text{\times}{10}^{6}\text{\,}\mathrm{A}\mathrm{/}\mathrm{(}\mathrm{{}^{2}}\mathrm{m}\mathrm{sr}\mathrm{V}\mathrm{)}$.
We are working on improving the long-term stability and repeatability of the
emission to make this field emitter a commercially viable electron source.

Having field emitted electrons with such a narrow energy distribution coming
from a nanoscale spatial location raises the question about the ultimate limit
in the width of the energy spectrum. For instance, in FN field emission at
cryogenic temperatures the transverse and longitudinal uncertainties in
momentum are $\Delta p_{t}\approx\Delta p_{l}$ Peter W. Hawkes and Erwin
Kasper (2017). As an estimate, we set the size of the emission region to be
the uncertainty in the position $\Delta x\approx r_{0}=$
$1.3\text{\,}\mathrm{n}\mathrm{m}$ and we approximate the uncertainty in the
momentum from the total energy distribution of the $\Delta
p\approx\sqrt{2m\Delta E}=$
$4.3\text{\times}{10}^{-16}\text{\,}\mathrm{e}\mathrm{V}\mathrm{s}\mathrm{/}\mathrm{n}\mathrm{m}$
from our measured $\Delta E=$ $16\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$.
Then we have $\Delta x\Delta p=1.7\times\frac{\hbar}{2}$, with the caveat that
we used a calculated value of $\Delta x$. It reveals that the uncertainty
product of our source is already very close, within a factor of two, to the
Heisenberg limit.

In summary, we realized a superconducting monocrystalline niobium nanotip
field emitter at a temperature of $5.9\text{\,}\mathrm{K}$ with an
unprecedented ultra-narrow tuneable energy distribution down to
$16\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ FWHM, a small emission angle of
$3.2\text{\,}\mathrm{\SIUnitSymbolDegree}$ and a brightness of up to
$5.0\text{\times}{10}^{8}\text{\,}\mathrm{A}\mathrm{/}\mathrm{(}\mathrm{{}^{2}}\mathrm{m}\mathrm{sr}\mathrm{V}\mathrm{)}$.
This was achieved by the fabrication of a nano-protrusion at the tip apex,
leading to a localized quantum state that mediates resonant tunneling through
the Coulomb barrier beyond the Fowler-Nordheim model. Although the tip
emission varies slightly after each preparation cycle, this low energy
distribution was measured several times following tip preparation and could be
reproduced with a completely new ion milled tip. The energy width is in the
same range as the recently developed flat surface source based on near-
threshold photoemission from single-crystal Cu(100) at $35\text{\,}\mathrm{K}$
with $\Delta E=$11.5\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ Karkare  _et
al._ (2020). However, the emission area of our source is orders of magnitude
smaller, providing potential advantages in spatial coherence and brightness.
It opens up new fields such as high-resolution vibrational spectroscopy in the
electron microscope Krivanek  _et al._ (2014) with sub-meV resolution when
combined with a monochromator. It will improve electron energy loss
spectroscopy (EELS), decrease the impact of chromatic aberrations in low-
voltage scanning electron microscopes and enhance the accuracy in the study of
band gaps or defects to the sub-nanometer level. The small spatial dimension
of the beam emission area in the nanometer regime combined with the narrow
energy distribution yields a high longitudinal and transversal coherence of
this source. This will benefit emerging techniques such as multipass-
transmission electron microscopy Juffmann  _et al._ (2017) and quantum
electron microscopy Kruit  _et al._ (2016) as well as quantum information
science applications with coherent electrons Röpke  _et al._ (2021) or
metrology with electron interferometers Pooch  _et al._ (2017).

###### Acknowledgements.

We thank Alexei Fedorov for his fruitful assistance. Work at the Molecular
Foundry was supported by the Office of Science, Office of Basic Energy
Sciences, of the U.S. Department of Energy under Contract No. DE-
AC02-05CH11231. This material is based upon work supported by the U.S.
Department of Energy, Small Business Innovation Research (SBIR) / Small
Business Technology Transfer (STTR) program under Award Number DE-SC-0019675.
We also acknowledge support by the Deutsche Forschungsgemeinschaft through the
research grant STI 615/3-1.

## References

  * Kruit  _et al._ (2016) P. Kruit, R. Hobbs, C.-S. Kim, Y. Yang, V. Manfrinato, J. Hammer, S. Thomas, P. Weber, B. Klopfer, C. Kohstall, T. Juffmann, M. Kasevich, P. Hommelhoff, and K. Berggren, Ultramicroscopy 164, 31 (2016).
  * Turner  _et al._ (2021) A. E. Turner, C. W. Johnson, P. Kruit, and B. J. McMorran, Physical Review Letters 127, 110401 (2021).
  * Küster  _et al._ (2021) F. Küster, A. M. Montero, F. S. M. Guimarães, S. Brinker, S. Lounis, S. S. P. Parkin, and P. Sessi, Nature Communications 12, 1108 (2021).
  * Krivanek  _et al._ (2014) O. L. Krivanek, T. C. Lovejoy, N. Dellby, T. Aoki, R. Carpenter, P. Rez, E. Soignard, J. Zhu, P. E. Batson, M. J. Lagos, _et al._ , Nature 514, 209 (2014).
  * Madan  _et al._ (2020) I. Madan, G. M. Vanacore, S. Gargiulo, T. LaGrange, and F. Carbone, Applied Physics Letters 116, 230502 (2020).
  * Ehberger  _et al._ (2015) D. Ehberger, J. Hammer, M. Eisele, M. Krüger, J. Noe, A. Högele, and P. Hommelhoff, Physical Review Letters 114, 227601 (2015).
  * Kuo  _et al._ (2006) H.-S. Kuo, I.-S. Hwang, T.-Y. Fu, Y.-C. Lin, C.-C. Chang, and T. T. Tsong, Japanese Journal of Applied Physics 45, 8972 (2006).
  * Pooch  _et al._ (2018) A. Pooch, M. Seidling, N. Kerker, R. Röpke, A. Rembold, W. T. Chang, I. S. Hwang, and A. Stibor, Physical Review A 97, 013611 (2018).
  * Zhang  _et al._ (2022) H. Zhang, Y. Jimbo, A. Niwata, A. Ikeda, A. Yasuhara, C. Ovidiu, K. Kimoto, T. Kasaya, H. T. Miyazaki, N. Tsujii, _et al._ , Nature nanotechnology 17, 21 (2022).
  * Johnson  _et al._ (2020) C. W. Johnson, J. S. Pierce, R. C. Moraski, A. E. Turner, A. T. Greenberg, W. S. Parker, and B. J. McMorran, Optics Express 28, 17334 (2020).
  * Feist  _et al._ (2020) A. Feist, S. V. Yalunin, S. Schäfer, and C. Ropers, Physical Review Research 2, 043227 (2020).
  * Hammer  _et al._ (2015) J. Hammer, S. Thomas, P. Weber, and P. Hommelhoff, Physical Review Letters 114, 254801 (2015).
  * Schütz  _et al._ (2014) G. Schütz, A. Rembold, A. Pooch, S. Meier, P. Schneeweiss, A. Rauschenbeutel, A. Günther, W. Chang, I. Hwang, and A. Stibor, Ultramicroscopy 141, 9 (2014).
  * Johnson  _et al._ (2021) C. W. Johnson, A. E. Turner, and B. J. McMorran, Physical Review Research 3, 043009 (2021).
  * Kerker  _et al._ (2020) N. Kerker, R. Röpke, L. M. Steinert, A. Pooch, and A. Stibor, New Journal of Physics 22, 063039 (2020).
  * Beierle  _et al._ (2018) P. J. Beierle, L. Zhang, and H. Batelaan, New Journal of Physics 20, 113030 (2018), publisher: IOP Publishing.
  * Latychevskaia  _et al._ (2017) T. Latychevskaia, W.-H. Hsu, W.-T. Chang, C.-Y. Lin, and I.-S. Hwang, Nature Communications 8, 14440 (2017).
  * Pooch  _et al._ (2017) A. Pooch, M. Seidling, M. Layer, A. Rembold, and A. Stibor, Applied Physics Letters 110, 223108 (2017).
  * Röpke  _et al._ (2021) R. Röpke, N. Kerker, and A. Stibor, New Journal of Physics 23, 023038 (2021).
  * Huang  _et al._ (2020) H. Huang, C. Padurariu, J. Senkpiel, R. Drost, A. L. Yeyati, J. C. Cuevas, B. Kubala, J. Ankerhold, K. Kern, and C. R. Ast, Nature Physics 16, 1227 (2020).
  * Forbes  _et al._ (2013) R. G. Forbes, A. Fischer, and M. S. Mousa, Journal of Vacuum Science & Technology B 31, 02B103 (2013).
  * Peter W. Hawkes and Erwin Kasper (2017) Peter W. Hawkes and Erwin Kasper, “Principles of Electron Optics, Volume 2 - 2nd Edition,”  (2017).
  * Binh  _et al._ (1992) V. T. Binh, S. T. Purcell, N. Garcia, and J. Doglioni, Physical Review Letters 69, 2527 (1992).
  * Gohda and Watanabe (2001) Y. Gohda and S. Watanabe, Physical Review Letters 87, 177601 (2001).
  * (25) Surface Preparation Laboratory B.V .
  * Hou  _et al._ (2016) J.-L. Hou, W.-T. Chang, C.-C. Shih, Y.-F. Yu, T.-Y. Fu, and I.-S. Hwang, Applied Physics Letters 108, 263107 (2016).
  * (27) Thermo Fisher Scientific, model FEI Helios G4 UX .
  * (28) Measured with an optical pyrometer, Leeds & Northrup Co., Cat. No. 8632-C .
  * (29) Advanced Research Systems, model DE-210 .
  * (30) ScientaOmicron DA20 R .
  * (31) Kimball Physics, model FC-71 .
  * (32) Keithley, model 237 .
  * (33) RoentDek, model DLD HEX100 .
  * Young (1959) R. D. Young, Physical Review 113, 110 (1959), publisher: American Physical Society.
  * Lin  _et al._ (1992) M. E. Lin, R. Reifenberger, and R. P. Andres, Physical Review B 46, 15490 (1992).
  * Saenz  _et al._ (1990) J. J. Saenz, N. Garcia, T. B. Vu, and H. De Raedt, in  _Scanning Tunneling Microscopy and Related Methods_, NATO ASI Series, edited by R. J. Behm, N. Garcia, and H. Rohrer (Dordrecht, 1990) pp. 409–441.
  * Börrnert  _et al._ (2018) F. Börrnert, J. Renner, and U. Kaiser, Microscopy and Microanalysis 24, 249 (2018).
  * Shao  _et al._ (2018) X. Shao, A. Srinivasan, W. K. Ang, and A. Khursheed, Nature communications 9, 1 (2018).
  * Nagaoka  _et al._ (2001) K. Nagaoka, H. Fujii, K. Matsuda, M. Komaki, Y. Murata, C. Oshima, and T. Sakurai, Applied Surface Science 182, 12 (2001).
  * Mook and Kruit (2000) H. Mook and P. Kruit, Ultramicroscopy 81, 129 (2000).
  * Karkare  _et al._ (2020) S. Karkare, G. Adhikari, W. A. Schroeder, J. K. Nangoi, T. Arias, J. Maxson, and H. Padmore, Physical Review Letters 125, 054801 (2020).
  * Juffmann  _et al._ (2017) T. Juffmann, S. A. Koppell, B. B. Klopfer, C. Ophus, R. M. Glaeser, and M. A. Kasevich, Scientific reports 7, 1 (2017).


# Spatio-Temporal Decomposition of Sum-of-Squares Programs for the Region of
Attraction and Reachability

Vít Cibulka1,2, Milan Korda1,2 and Tomáš Haniš1 *This research was supported
by the Czech Science Foundation (GACR) under contracts No. GA19-18424S,
GA20-11626Y, and by the Grant Agency of the Czech Technical University in
Prague, grant No. SGS19/174/OHK3/3T/13. This work has also been supported by
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie Actions, grant agreement 813211 (POEMA) and by the AI
Interdisciplinary Institute ANITI funding, through the French “Investing for
the Future PIA3” program under the Grant agreement n∘ ANR-19-PI3A-0004.
1Department of Control Engineering, Faculty of Electrical Engineering, Czech
Technical University in Prague, The Czech Republic<EMAIL_ADDRESS><EMAIL_ADDRESS>CNRS, Laboratory for Analysis and Architecture of
Systems, Toulouse, France<EMAIL_ADDRESS>work has been submitted to the
IEEE for possible publication. Copyright may be transferred without notice,
after which this version may no longer be accessible.

###### Abstract

This paper presents a method for calculating Region of Attraction of a target
set (not necessarily an equilibrium) for controlled polynomial dynamical
systems, using a hierarchy of semidefinite programming problems (SDPs). Our
approach builds on previous work and addresses its main issue, the fast-
growing memory demands for solving large-scale SDPs.

The main idea in this work is in dissecting the original resource-demanding
problem into multiple smaller, interconnected, and easier to solve problems.
This is achieved by spatio-temporal splitting akin to methods based on partial
differential equations. We show that the splitting procedure retains the
convergence and outer-approximation guarantees of the previous work, while
achieving higher precision in less time and with smaller memory footprint.

##  I Introduction

This paper deals with stability and reachability analysis of nonlinear
dynamical systems. A prominent method for quantifying the stability of a
nonlinear system is based on calculating its Region Of Attraction (ROA) with
respect to a given target set, which is the topic of this work. Using a time-
reversal, one obtains the reachable set. Both of these objects are crucial in
assessing stability and safety of control system; as a concrete example, let
us mention the Flight Envelopes in the aerospace industry that directly
coincide with the ROA [1]. A more recent application would be the study of
closed-loop parameters variation, of both the controller and the controlled
plant, and its effect on the ROA. For example, given a vehicle on the road,
these parameters could be the distribution of mass of the vehicle or various
road conditions [2].

The dominant methods for computing ROA are based on Lyapunov functions whose
level sets provide inner approximations the ROA [3]. For polynomial systems,
the Lyapunov function can be found by solving semidefinite programming
problems (SDPs) [4] and [5]. These methods only work for autonomous systems
and ROA’s with respect to given equilibria.

This work expands the approach presented in [6], which is based on optimizing
over trajectories of a polynomial dynamical systems modelled as occupation
measures and is not limited to an autonomous system nor a stable equilibrium.
The contribution proposed in this paper is dissection of the original problem
into multiple interconnected problems of lower complexity. The original
formulation is discretized in both the time and state variables, aiming to
strike a trade-off between the original optimization-based approach [6] where
a single polynomial is defined on the entire state-space and various methods
based on extensive discretization such as the PDE-based approach [7], set
oriented methods [8] or transfer operator approaches [9]. By doing so, our
approach, based fully on convex optimization, preserves the outer
approximation and convergence guarantees of [6], both of which are explicitly
proven in this work, while significantly improving scalability of the method.

Structure of this paper Section II presents the problem statement and Section
III introduces the splitting procedure with the proofs of outer approximation
and convergence. Section IV states the practical sum-of-squares variant of the
problem which is then demonstrated in Section V on numerical examples. The
paper ends with a conclusion and future work discussion in sections VI and VII
respectively.

Notation The symbol $\mathbb{Z}_{j}$ denotes the set of consecutive integers
$\\{1,2,\dots,j\\}$. The Lebesgue measure (i.e., the volume) of a set is
denoted by $\lambda(A)$. The indicator function $I_{A}(x)$ of a set $A$ is the
function that takes the value $1$ for $x\in A$, and $0$ otherwise. The
interior of a set $A$ is denoted by $A^{\circ}$. The symbols $C(A)$ and
$C^{1}(A)$ denote respectively the sets of continuous and continuously
differentiable functions on $A$. The set of Borel measurable functions defined
on a set $A$ taking values in a set $B$ is denoted by $L(A;B)$

##  II Problem statement

Let us consider the nonlinear system with control

$\dot{x}(t)=f(t,x(t),u(t)),t\in[0,T],$ (1)

where $x(t)\in\mathbb{R}^{n}$ is the state vector, $u(t)\in\mathbb{R}^{m}$ is
the control input vector, $t$ is time, $T>0$ is the final time and $f$ is the
vector field, which is assumed to be polynomial in variables $x$ and $u$.

The state and control input are constrained by basic semialgebraic sets

$\displaystyle u(t)\in U:=$
$\displaystyle\\{u\in\mathbb{R}^{m}:g_{j}^{U}(u)\geq
0,j\in\mathbb{Z}_{n_{U}}\\},t\in[0,T],$ (2) $\displaystyle x(t)\in X:=$
$\displaystyle\\{x\in\mathbb{R}^{n}:g_{j}^{X}(x)\geq
0,j\in\mathbb{Z}_{n_{X}}\\},t\in[0,T],$ $\displaystyle x(T)\in X_{T}:=$
$\displaystyle\\{x\in\mathbb{R}^{n}:g_{j}^{X_{T}}(x)\geq
0,j\in\mathbb{Z}_{N_{X_{T}}}\\},$

where $g_{j}^{U}(u)$, $g_{j}^{X}(x)$, and $g_{j}^{X_{T}}(x)$ are polynomials.
The region of attraction (ROA) is then defined as

$\displaystyle X_{0}=\\{$ $\displaystyle x_{0}\in$ $\displaystyle
X:\exists\,u(\cdot)\in L([0,T];U)$ (3) s.t.
$\displaystyle\dot{x}=f(t,x(t),u(t))\text{ a.e. on}\;[0,T],$ $\displaystyle
x(0)=x_{0},\;x(t)\in X\;\forall\,t\in[0,T],\;x(T)\in X_{T}\\},$

where “a.e.” stands for “almost everywhere” with respect to the Lebesgue
measure.

##  III Time and state space splitting

It was shown in [6], that the ROA can be characterized by an infinite-
dimensional linear programming (LP) problem in the space of Borel measures or
by its dual in the space of continuous functions. These LPs can then be
approximated by SDPs, with guaranteed convergence. For space reason, the
entire exposition in this work will be in the “dual” setting of continuous
functions approximated by polynomials.

Let us restate the original problem from [6]:

$\displaystyle d^{\star}$ $\displaystyle=\inf\int_{X}w(x)d\lambda(x)$ , (4)
s.t. $\displaystyle(\mathcal{L}v)(t,x,u)\leq 0$
$\displaystyle,\forall(t,x,u)\in[0,T]\times X\times U$ $\displaystyle w(x)\geq
v(0,x)+1$ $\displaystyle,\forall x\in X$ $\displaystyle v(T,x)\geq 0$
$\displaystyle,\forall x\in X_{T}$ $\displaystyle w(x)\geq 0$
$\displaystyle,\forall x\in X,$

with variables $w(x)\in C(X)$ and $v(t,x)\in C([0,T]\times X)$.

Any minimizing sequence $(w_{k},v_{k})$ for (4) satisfies $w_{k}\geq
I_{X_{0}}(x)$ and $w_{k}\to I_{X_{0}}$ in $L_{1}$ as well as $\\{x\in
X:v_{k}(0,x)\geq 0\\}\supset X_{0}$ with convergence in terms of the volume
discrepancy tending to zero (see [6] for proof).

Let us now split the state space $X$ into $I$ closed subsets $X_{i}$

$X=\bigcup\limits_{i=1}^{I}X_{i}$ (5)

and the time interval $[0,T]$ into $K-1$ intervals $[T_{k},T_{k+1}]$

$[0,T]=\bigcup\limits_{k=1}^{K-1}[T_{k},T_{k+1}],$ (6)

where $K$ is the number of time splits (meaning $K-1$ intervals). It is
assumed that $X_{i}^{\circ}\cap X_{j}^{\circ}=\emptyset$ for $i\neq j$.

The function $w(x)$ will be split into $I$ functions $w_{i}(x)$

$w(x)=\begin{cases}w_{1}(x)&\text{for }x\in X_{1}\\\ \dots\\\
w_{i}(x)&\text{for }x\in X_{i}\\\ \dots\\\ w_{I}(x)&\text{for }x\in
X_{I}\end{cases}$ (7)

and $v(t,x)$ will be split into $I\cdotp(K-1)$ functions $v_{i,k}(t,x)$

$v(t,x)=\begin{cases}v_{1,1}(t_{1},x_{1})&\text{for }t\in[T_{1},T_{2}],x\in
X_{1}\\\ \dots\\\ v_{i,k}(t_{k},x_{i})&\text{for }t\in[T_{k},T_{k+1}],x\in
X_{i}\\\ \dots\\\ v_{I,K-1}(t_{K-1},x_{I})&\text{for }t\in[T_{K-1},T_{K}],x\in
X_{I}.\end{cases}$ (8)

Assuming that neighbouring subsets $X_{i}$ share boundaries, let us define the
set of indices of these neighbours as

$\displaystyle N_{X}:=\\{(a,b):X_{a}\cap X_{b}\not=\emptyset\\}.$ (9)

We are now ready to write the split version of (4)

$\begin{split}\begin{aligned}
&d_{\mathrm{s}}^{\star}=&&\inf\sum_{i}\int_{X_{i}}w_{i}(x)d\lambda(x)\\\
&\text{s.t.\,}&&\text{for all }i\in\mathbb{Z}_{I},k\in\mathbb{Z}_{K-1}\text{
and }\\\ &&&\hskip 85.35826pt(a,b)\in N_{X},&&x_{a,b}\in X_{a}\cap X_{b}\\\
&&&(\mathcal{L}v_{i,k})(t,x,u)\leq
0\quad\,\forall(t,x,u)\in[&&T_{k},T_{k+1}]\times X_{i}\times U\\\
&&&w_{i}(x)\geq v_{i,k}(0,x)+1\quad&&\forall x\in X_{i}\\\ &&&v_{i,K}(T,x)\geq
0\quad&&\forall x\in X_{T}\\\ &&&w_{i}(x)\geq 0\quad&&\forall x\in X_{i}\\\
&&&v_{i,k}(T_{k+1},x)\geq v_{i,k+1}(T_{k+1},x)\quad&&\forall x\in X_{i}\\\
&&&(v_{a,k}(t,x_{a,b})-v_{b,k}(t,x_{a,b}))\cdot
h_{a,b}^{\top}&&f(t,x_{a,b},u)\geq 0\end{aligned}\\\ \end{split}$ (10)

where $h_{a,b}^{\top}$ is a normal vector of a shared boundary between two
neighbouring sets $X_{a}$ and $X_{b}$, and $x_{a,b}$ is a set of points on
said boundary such that $x_{a,b}\in X_{a}\cap X_{b}$. For simplicity we assume
that the normal vector $h_{a,b}$ is independent of $x$; the case of polynomial
or rational dependence of $h$ on $x$ can also be handled [10, Section 4.3].
The optimization variables in (10) are the continuously differentiable
functions of $v_{i,k}$, each defined some neighborhood of
$[T_{k},T_{k+1}]\times X_{i}$ and the continuous functions $w_{i}$, each
defined on $X_{i}$.

Let us now show that the modified problem (10) provides a guaranteed outer
approximation of $X_{0}$ which can be defined as

$\bar{X}_{v,0}:=\\{x:v(0,x)\geq 0\\}.$ (11)

###### Theorem 1

For any pair $(v,w)$ feasible in (10), it holds that $v(0,\cdot)\geq 0$ on
$X_{0}$ and $\bar{X}_{v,0}\supset X_{0}$.

###### Proof:

We first need to show that the discontinuous function $v(t,x(t))$ is
decreasing along the system trajectories. That is, given two time instants
$t_{\alpha}\leq t_{\beta}$, we want to show that

$v(t_{\beta},x(t_{\beta}))\leq v(t_{\alpha},x(t_{\alpha})).$ (12)

If we $v$ were to be differentiable, this follows by simply integrating the
first constraint of (4) along a trajectory.

Let us now show that (12) holds even for the discontinuous $v(t,x(t))$ as
defined in (8). The first constraint of (10) ensures that the values of $v$
decrease whenever the trajectory resides in the interior of one of the sets
$X_{i}$. We therefore need to argue only about what happens on the boundary of
these sets.

##### Time splits

The result for time splits follows immediately from the fifth constraint of
(10).

##### State-space splits

Let us first assume that the state space $X$ is split into two parts,
$X_{\alpha}$ and $X_{\beta}$ by a hyper-plane with normal vector $h$, pointing
from $X_{\alpha}$ to $X_{\beta}$, so that

$\displaystyle h^{\top}(x_{\beta}-x_{\alpha})\geq 0,$ (13)

for $x_{\alpha}\in X_{\alpha}$ and $x_{\beta}\in X_{\beta}$.

The function $v(t,x(t))$, now split between $X_{\alpha}$ and $X_{\beta}$, is
defined as

$v(t,x(t))=\begin{cases}v_{\alpha}(t,x(t))&\text{for }t\in[0,T],x(t)\in
X_{\alpha}\\\ v_{\beta}(t,x(t))&\text{for }t\in[0,T],x(t)\in
X_{\beta}.\end{cases}$ (14)

Let $x_{0}\in X_{0}\cap X_{\alpha}$, $t_{\alpha}\in[0,T]$ and $u(\cdot)$ be
given. Let $x(\cdot|x_{0})$ be the trajectory starting at $t_{\alpha}$
generated by $u(\cdot)$ and suppose that $x(t\mid x_{0})\in X$, $u(t)\in U$
for $t\in[t_{\alpha},T]$. Assume further that this trajectory crosses from
$X_{\alpha}$ to $X_{\beta}$ at the crossing time

$\tau=\inf_{t}\\{t_{\alpha}\leq t\mid x(t|x_{0})\in X^{\circ}_{\beta}\\}\leq
T$ (15)

and assume that this trajectory stays in $X_{\beta}$ for $t\in[\tau,T]$. At
the crossing point $x(\tau)$, it holds

$h^{\top}f(\tau,x(\tau),u(\tau))\geq 0.$

The last constraint of (10) implies that

$v_{\alpha}(\tau,x(\tau))\geq v_{\beta}(\tau,x(\tau))$ (16)

whereas the first constraint implies

$\frac{\mathrm{d}}{\mathrm{d}t}v_{\alpha}(t,x(t))\leq 0,t\in[t_{\alpha},\tau)$
(17)

and

$\frac{\mathrm{d}}{\mathrm{d}t}v_{\beta}(t,x(t))\leq 0,t\in(\tau,t_{\beta}].$
(18)

Let us now calculate the value of $v_{\beta}(t_{\beta},x(t_{\beta}))$:

$v_{\beta}(t_{\beta},x(t_{\beta}))=v_{\alpha}(t_{\alpha},x(t_{\alpha}))+\int_{t_{\alpha}}^{\tau}\frac{\mathrm{d}}{\mathrm{d}t}v_{\alpha}(t,x(t))\mathrm{d}t\\\
+[v_{\beta}(\tau,x(\tau))-v_{\alpha}(\tau,x(\tau))]+\int_{\tau}^{t_{\beta}}\frac{\mathrm{d}}{\mathrm{d}t}v_{\beta}(t,x(t))\mathrm{d}t.$
(19)

By inspecting (16), (17) and (17), we can see that the last three summands are
nonpositive and we get the inequality $v_{\beta}(t_{\beta},x(t_{\beta}))\leq
v_{\alpha}(t_{\alpha},x(t_{\alpha}))$, which is equivalent to (12), recalling
the definition of $v(t,x(t))$ in (14). The procedure for the negative
trajectory direction is analogous. We note that this analysis encompasses the
subtle case of the trajectory sliding on the boundary between the sets
$X_{\alpha}$ and $X_{\beta}$.

By induction, we can prove the inequality for arbitrary splitting of the
state-space and time axis by considering a sequence of crossing times
associated to a given trajectory. Therefore $v(t_{\beta},x(t_{\beta}))\leq
v(t_{\alpha},x(t_{\alpha}))$ for any $0\leq t_{\alpha}\leq t_{\beta}\leq T$.
By setting $t_{\alpha}=0,t_{\beta}=T$, and using the constraints of (10), we
get

$\displaystyle v(T,x(T))$ $\displaystyle\leq v(0,x_{0})$ (20) $\displaystyle
0\leq v(T,x(T))$ $\displaystyle\leq v(0,x_{0})$ $\displaystyle 0$
$\displaystyle\leq v(0,x_{0})$

for any $x_{0}\in X_{0}$ as desired. This also implies that
$x_{0}\in\bar{X}_{v,0}$ and hence $\bar{X}_{v,0}\supset X_{0}$.

∎

##  IV SOS representation

We can now obtain the SDP representation of (10) by applying Putinar’s
Positivstellensatz [11]. For example, given polynomials $c(x)$ and $g(x)$ the
inequality

$c(x)\geq 0\quad\text{for }x\in\\{x:g(x)\geq 0\\}$ (21)

is implied by

$c(x)=q(x)+s(x)g(x),$ (22)

where $q(x)$ and $s(x)$ are sum-of-squares polynomials. The condition that a
polynomial $s$ of degree $2d$ is sum-of-squares is in turn equivalent to
$s(x)=m(x)^{\top}Wm(x)$, $W\succeq 0$, where $m(x)$ is a basis of polynomials
up degree $d$ and hence this constraint is SDP representable.

The SOS approximation of (10) reads

$\displaystyle\inf\sum_{i}\text{ }\mathbf{w}_{i}^{\top}l_{i}$ (23) s.t.
$\displaystyle\text{for all }i\in\mathbb{Z}_{I},k\in\mathbb{Z}_{K-1}\text{ and
}(a,b)\in N_{X}$
$\displaystyle-(\mathcal{L}v_{i,k})(z)=q_{i,k}(z)+\bm{{s^{\tau}_{i,k}}}(z)^{\top}\bm{{g^{\tau}_{k}}}(t)$
$\displaystyle\hskip
62.59596pt+\bm{{s^{X}_{i,k}}}(z)^{\top}\bm{{g^{X}_{i,k}}}(x)+\bm{{s^{U}_{i,k}}}(z)^{\top}\bm{{g^{U}}}(u)$
$\displaystyle
w_{i}(x)-v_{i,k}(0,x)-1=q_{0_{i,k}}(x)+\bm{{s^{0}_{i,k}}}(x)^{\top}\bm{{g_{i}^{X}}}(x)$
$\displaystyle
v_{i,K}(T,x)=q^{T}_{i}(x)+\bm{{s^{X_{T}}_{i}}}(x)^{\top}\bm{{g^{X_{T}}_{i}}}(x)$
$\displaystyle
w_{i}(x)=q^{w}_{i}(x)+\bm{{s^{w}_{i}}}(x)^{\top}\bm{{g^{X}_{i}}}(x)$
$\displaystyle v_{i,k}(T_{k+1},x)-v_{i,k+1}(T_{k+1},x)=$ $\displaystyle\hskip
95.3169ptq_{{i,k}}^{\tau}(x)+\bm{{s^{t}_{i,k}}}(x)^{\top}\bm{{g^{X}_{i}}}(x)$
$\displaystyle(v_{a,k}(t,x_{a,b})-v_{b,k}(t,x_{a,b}))=$ $\displaystyle\hskip
16.78702ptq^{1}_{k,a,b}(z)+\sum\nolimits_{j=1}^{n_{X}}s^{1}_{j,k,a,b}(z)h_{a,b}^{\top}f(t,x_{a,b},u)$
$\displaystyle(v_{b,k}(t,x_{a,b})-v_{a,k}(t,x_{a,b}))=$ $\displaystyle\hskip
16.78702ptq^{2}_{k,a,b}(z)-\sum\nolimits_{j=1}^{n_{X}}s^{2}_{j,k,a,b}(z)h_{a,b}^{\top}f(t,x_{a,b},u)$

where $z=[t,x,u]^{\top}$, $w_{i}(x)$ and $v_{i,k}(t,x)$ are polynomials,
$\mathbf{w}_{i}$ is a vector of coefficients of $w_{i}(x)$ and $l_{i}$ is a
vector of Lebesgue measure moments indexed with respect to the same basis as
the coefficients of $w_{i}$. The decision variables in the problem are the
polynomials $v_{i,k}$ and $w_{i}$ as well as sum-of-squares multipliers $q$,
$s$ and $\bm{{s}}$. The symbols $\bm{g_{i}^{X}}$, $\bm{g_{i}^{X_{T}}}$,
$\bm{g_{i}^{U}}$ and $\bm{g_{k}^{\tau}}$ denote the column vectors of
polynomials describing the sets $X_{i}$, $X_{T}\cap X_{i}$, $U$ and
$[T_{k},T_{k+1}]$ in that order. The degrees of all polynomial decision
variables is chosen such that the degrees of all polynomials appearing in (23)
do not exceed a given relaxation order $d$. This is a design parameter
controlling the accuracy of the approximation.

Given the picewise polynomial functions $(w^{d},v^{d})$ of degree $d$
constructed from a solution to (23) as in (7) and (8), the outer approximation
to the ROA is defined by

$X_{d}=\\{x\mid v^{d}(0,x)\geq 0\\}.$

Convergence of the SDP approximations holds under the classical Archimedianity
assumption, e.g., [6, Assumption 3].

###### Theorem 2

For each $d\in\mathbb{N}$, we have $X_{d}\supset X_{0}$. If in addition the
algebraic description of each element of the space-time partition used in (23)
satisfies the Archimedianity condition, then
$\lim_{d\to\infty}\lambda(X_{d}\setminus X_{0})=0$.

###### Proof:

The result follows from [6, Theorem 6] by taking all functions defined on the
partition equal, i.e., $v_{i,k}=v$ and $w_{i}=w$ for some polynomials $v$ and
$w$; this leads to the setting of [6, Theorem 6]. ∎

###  IV-A Practical implications

The ROA with splits is expected to improve accuracy of the original
formulation by allowing one to trade off the degree of the polynomials for
number of splits. By increasing the degree $d$, the size of the SDP will
increase with the rate of the binomial coefficient ${n+d}\choose{n}$. By
fixing $d$ and splitting the state-space into $\lambda$ cells, the SDP size
will grow only linearly with the number of the cells with rate
$\lambda{n+d\choose{n}}$. The numerical results in the following section
suggest that even the computation time grows linearly. A full quantitative
understanding of this tradeoff is beyond the scope of the present work; the
first step in this direction would be the convergence rate bounds of [12].

##  V Numerical examples

This section presents numerical examples, showcasing the performance
difference between the proposed method and the original approach from [6].

The first example in Section V-A shows the influence of the split positions on
the resulting ROA. It is shown that by having the splits exactly at the
boundaries of the ROA, we can retrieve the theoretical indicator function
$I_{X_{0}}(x)$.

The second example Section V-B benchmarks the algorithm on a Brockett
integrator, which mimics a kinematic model of a nonholonomic system (it can be
shown that three-dimensional nonholonomic vehicle with two inputs can be
transformed into the Brockett integrator [13]).

Finally, Section V-C presents a comparison of computational demands of the
proposed method and the original one from [6].

All the examples were implemented in MATLAB [14] with the use of YALMIP [15].
The YALMIP’s sum-of-squares package [16] was used for rapid prototyping; for
larger examples, the SDPs were assembled using a custom routine. All SDP’s
were solved by MOSEK [17].

###  V-A Univariate cubic dynamics

This example shows that one can find the ROA with a very low degree
polynomials by correctly positioning the splits.

The system in question is defined as

$\dot{x}=x(x-0.5)(x+0.5)$ (24)

with the state space $X=[-1,1]$, the target set $X_{T}=[-0.01,0.01]$ and
terminal time $T=100$. The analytic solution of the ROA is $X_{0}=[-0.5,0.5]$.

In Fig. 1, we can see a comparison between the original method (without
splits) and multiple calculations with splits, going from the inside of the
real ROA to the outside. The ROA estimates here are given by $\\{x:w(x)\geq
1\\}$, which follows from (11) and (10). We can observe that the ROA estimates
get more precise, the closer the splits are to the real ROA. Let us define the
estimate of $I_{X_{0}}(x)$ as $\bar{I}_{X_{0}}(x)$, which takes $1$ on
$\\{x:w(x)\geq 1\\}$ and $0$ otherwise. We can observe from Fig. 1, that for
the exact split it holds that $\bar{I}_{X_{0}}=I_{X_{0}}$ and we obtain the
theoretically optimal estimate. This example shows that our method can be used
in an iterative manner with splits along an inner approximation of the ROA
(such as the one in [18]) as a starting point.

Figure 1: Univariate cubic dynamics, approximations of the indicator function
with degree 8 polynomials and two splits. The ROA is given by $\\{x:w(x)\geq
1\\}$. The approximation is more precise for splits that are close to the
bounds of the ROA (red), and gives the exact ROA when the splits are exactly
at the bounds (blue). Both split polynomials are more precise than the green,
non-split one.

###  V-B Brockett integrator

The Brockett integrator is defined according to [19] as

$\displaystyle\dot{x}_{1}$ $\displaystyle=u_{1}$ (25)
$\displaystyle\dot{x}_{2}$ $\displaystyle=u_{2}$ $\displaystyle\dot{x}_{3}$
$\displaystyle=u_{1}x_{2}-u_{2}x_{1}.$

With $X=\\{x\in\mathbb{R}^{3}:||x||_{\infty}\leq 1\\}$, $X_{T}=\\{0\\}$,
$U=\\{u\in\mathbb{R}^{2}:||u||_{2}\leq 1\\}$, and $T=1$. As was stated before,
this system usually serves as a benchmark for nonholonomic control strategies,
because it is the simplest system for which there exists no continuous control
law which would make the origin asymptotically stable [19].

We shall use the system for calculation of the controlled ROA, which can be
computed analytically [20] as

$\displaystyle\mathcal{T}(x)=\frac{\theta\sqrt{x_{1}^{2}+x_{2}^{2}+2|x_{3}|}}{\sqrt{\theta+\sin^{2}\theta-\sin\theta\cos\theta}},$
(26)

where $\theta=\theta(x)$ is the unique solution in $[0,\pi)$ to

$\displaystyle\frac{\theta-\sin\theta\cos\theta}{\sin^{2}\theta}(x_{1}^{2}+x_{2}^{2})=2|x_{3}|.$
(27)

The Fig. 2 shows that given a fixed time for the calculation, the proposed
approach is always better than the original one and that the split-method
approaches the real volume much faster, although neither of the two methods
reached the real volume, due to memory constraints.

A visual example of the difference between the two methods can be seen in the
Fig. 3 where two ROA’s with the same computation time are compared, with one
being calculated by the original method and the other by the proposed method.
We can see that there is a notable difference between the two approximations,
and that the better approximation is done by the lower-degree polynomials.

Figure 2: Comparison of various methods on the Brockett integrator. The blue
line shows the original algorithm without splits and with increasing degree of
the approximation $d_{\text{o}}\in\\{6,8,10,12\\}$. The other lines show the
performance of the proposed approach, each having a fixed degree
$d_{\text{s}}\in\\{6,8\\}$ with increasing numbers of cells $X_{i}$, which are
denoted as numbers next to the datapoints. The volumes were estimated using
Monte-Carlo methods. Figure 3: Sliced ROA of the Brockett integrator (green)
and two slices of its approximations with similar computation time. Blue
approximation is by single degree 10 polynomial (754s) and the red is by
sixteen polynomials of degree 8 (753s). The red, lower-degree, approximation
is visibly closer to the real ROA.

###  V-C Performance and scalability

####  V-C1 Problem size

First, we shall investigate the accuracy of the algorithm with increasing size
of the SDP. The problem size is measured as the number of nonzero elements in
the $A$ matrix of the SDP

$\displaystyle\min\quad c^{\top}x$ (28) $\displaystyle\text{s.t. }Ax=b,\quad
x\in\mathcal{K}$

for variable $x\in\mathbb{R}^{n}$, convex cone $\mathcal{K}$ and data
$A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$, and $c\in\mathbb{R}^{n}$.

We can clearly see in Fig. 4 that the split versions are always more precise
than the non-split version with the same memory footprint.

Figure 4: Brockett integrator - Given SDP size, the split problem always give
better results. The volumes were estimated using Monte-Carlo methods.

####  V-C2 Computation time

The problem size increases linearly with the number of cells (as was explained
in IV-A), but the computation time does not necessarily have to follow the
same pattern. In this case, however, the computation time also showed linear
growth as can be seen in the Figures 5 and 6 for the Brockett integrator and
the Double integrator respectively. The Double integrator is defined as

$\dot{x}_{1}=x_{2},\dot{x}_{2}=u$

with $X=[-0.7\times 0.7]\times[-1.2\times 1.2]$, $X_{T}=\\{0\\}$ and $T=1$.
See [6, 9.3] for more details. In the Fig. 6, the variables to be split were
chosen randomly and the splits were always halving the largest interval of the
randomly selected variable. This was done in order to ensure that the linear
growth is not simply a fortunate result of a particular split order.

Figure 5: Brockett integrator - The computation time increases linearly with
the number of cells $X_{i}$. Figure 6: Double integrator - The computation
time grows linearly with the number of randomly chosen cells $X_{i}$.

##  VI Conclusion

We showed that the convex SOS-based calculation of ROA can be extended by
splitting the time and state space to achieve better accuracy with lower
degree polynomials, while keeping the outer-approximation guarantees. The
method is faster and provides more accurate estimates of the ROA. It was also
empirically demonstrated that the computation time grows linearly with
increasing number of cells.

##  VII Future work

The sparse problem structure can be exploited by ADMM-like algorithm, where
the ADMM could alternate between solving small independent SDPs, and enforcing
the flow constraints between them. This would allow for parallelization and
significantly decrease the required computation time.

The splits can be found in an optimization-based fashion by differentiating
the SDP [21] with respect to the location (or more general parametrization of
the splits), thereby increasing the accuracy while keeping the memory
requirements constant. This also leads to the option of creating the splits in
a less naive fashion than the one presented here. Using inner-approximation of
the ROA lends itself as an obvious candidate, due to the behaviour shown in
V-A.

The same splitting approach can also be applied to other SOS-based algorithms,
e.g., to the optimal control problem [22], where the increased precision would
lead to more precise estimates of the optimal control sequences and also allow
for control of higher dimensional systems. When applied to the problem of the
maximum control invariant set estimation [23], temporal splitting is no longer
possible but could be replaced by introducing multiple discount factors.

Finally, a natural way to decrease the computational complexity further is to
combine the approach with the sparsity-based decomposition method of [24],
where the time and state-space splitting would be applied to the subsystems
that cannot decomposed using the method of [24].

## References

  * [1] G. Yuan and Y. Li, “Determination of the flight dynamic envelope via stable manifold,” _Measurement and Control_ , vol. 52, no. 3-4, pp. 244–251, feb 2019.
  * [2] M. I. Masouleh and D. J. N. Limebeer, “Region of attraction analysis for nonlinear vehicle lateral dynamics using sum-of-squares programming,” _Vehicle System Dynamics_ , vol. 56, no. 7, pp. 1118–1138, dec 2017.
  * [3] H. K. Khalil, _Nonlinear Systems_. Upper Saddle River, NJ: Prentice Hall, 2002.
  * [4] G. Chesi, A. Garulli, A. Tesi, and A. Vicino, “LMI-based computation of optimal quadratic lyapunov functions for odd polynomial systems,” _International Journal of Robust and Nonlinear Control_ , vol. 15, no. 1, pp. 35–49, 2004.
  * [5] B. Tibken, “Estimation of the domain of attraction for polynomial systems via LMIs,” in _Proceedings of the 39th IEEE Conference on Decision and Control (Cat. No.00CH37187)_. IEEE, 2000\.
  * [6] D. Henrion and M. Korda, “Convex computation of the region of attraction of polynomial control systems,” _IEEE Transactions on Automatic Control_ , vol. 59, no. 2, pp. 297–312, 2014.
  * [7] I. M. Mitchell and C. J. Tomlin, “Overapproximating reachable sets by hamilton-jacobi projections,” _Journal of Scientific Computing_ , vol. 19, no. 1/3, pp. 323–346, 2003.
  * [8] M. Dellnitz, G. Froyland, and O. Junge, “The algorithms behind GAIO — set oriented numerical methods for dynamical systems,” in _Ergodic Theory, Analysis, and Efficient Simulation of Dynamical Systems_. Springer Berlin Heidelberg, 2001, pp. 145–174.
  * [9] K. Wang, “Transfer operator-based approach for domain of attraction computation and experimental data analysis,” Ph.D. dissertation.
  * [10] M. Korda, D. Henrion, and J.-B. Lasserre, “Moments and convex optimization for analysis and control of nonlinear partial differential equations,” _arXiv preprint arXiv:1804.07565_ , 2018.
  * [11] M. Putinar, “Positive polynomials on compact semi-algebraic sets,” _Indiana University Mathematics Journal,_ , vol. 42, pp. 969–984, 1993.
  * [12] M. Korda, D. Henrion, and C. N. Jones, “Convergence rates of moment-sum-of-squares hierarchies for optimal control problems,” _Systems & Control Letters_, vol. 100, pp. 1–5, 2017.
  * [13] A. P. Aguiar and A. Pascoal, “Stabilization of the extended nonholonomic double integrator via logic-based hybrid control,” _IFAC Proceedings Volumes_ , vol. 33, no. 27, pp. 351–356, 2000, 6th IFAC Symposium on Robot Control (SYROCO 2000), Vienna, Austria, 21-23 September 2000. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1474667017379545
  * [14] MATLAB, _version 9.9.0 (R2020b)_. Natick, Massachusetts: The MathWorks Inc., 2020.
  * [15] J. Löfberg, “Yalmip : A toolbox for modeling and optimization in matlab,” in _In Proceedings of the CACSD Conference_ , Taipei, Taiwan, 2004\.
  * [16] ——, “Pre- and post-processing sum-of-squares programs in practice,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 5, pp. 1007–1011, 2009.
  * [17] M. ApS, _The MOSEK optimization toolbox for MATLAB manual. Version 9.0._ , 2019\. [Online]. Available: http://docs.mosek.com/9.0/toolbox/index.html
  * [18] M. Korda, D. Henrion, and C. N. Jones, “Inner approximations of the region of attraction for polynomial dynamical systems,” Oct. 2012.
  * [19] R. W. Brockett, “Asymptotic stability and feedback stabilization,” in _Differential Geometric Control Theory_. Birkhauser, 1983, pp. 181–191.
  * [20] J.-B. Lasserre, D. Henrion, C. Prieur, and E. Trélat, “Nonlinear optimal control via occupation measures and lmi-relaxations,” Mar. 2007.
  * [21] A. Agrawal, S. Barratt, S. Boyd, E. Busseti, and W. M. Moursi, “Differentiating through a cone program,” Apr. 2019.
  * [22] J. B. Lasserre, D. Henrion, C. Prieur, and E. Trélat, “Nonlinear optimal control via occupation measures and LMI-relaxations,” _SIAM Journal on Control and Optimization_ , vol. 47, no. 4, pp. 1643–1666, jan 2008.
  * [23] M. Korda, D. Henrion, and C. N. Jones, “Convex computation of the maximum controlled invariant set for polynomial control systems,” _SIAM Journal on Control and Optimization_ , vol. 52, no. 5, pp. 2944–2969, 2014.
  * [24] C. Schlosser and M. Korda, “Sparse moment-sum-of-squares relaxations for nonlinear dynamical systems with guaranteed convergence,” _arXiv preprint arXiv:2012.05572_ , 2020.