# NLI4CT: Multi-Evidence Natural Language Inference for Clinical Trial Reports

Maël Jullien1, Marco Valentino2, 1, Hannah Frost1, 3, Paul O’Regan3, Donal
Landers3, André Freitas1, 2
Department of Computer Science, University of Manchester, United Kingdom1
Idiap Research Institute, Switzerland2
Digital Experimental Cancer Medicine Team, Cancer Research UK Manchester
Institute3
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

How can we interpret and retrieve medical evidence to support clinical
decisions? Clinical trial reports (CTR) amassed over the years contain
indispensable information for the development of personalized medicine.
However, it is practically infeasible to manually inspect over 400,000+
clinical trial reports in order to find the best evidence for experimental
treatments. Natural Language Inference (NLI) offers a potential solution to
this problem, by allowing the scalable computation of textual entailment.
However, existing NLI models perform poorly on biomedical corpora, and
previously published datasets fail to capture the full complexity of inference
over CTRs.

In this work, we present a novel resource to advance research on NLI for
reasoning on CTRs. The resource includes two main tasks. Firstly, to determine
the inference relation between a natural language statement, and a CTR.
Secondly, to retrieve supporting facts to justify the predicted relation. We
provide NLI4CT, a corpus of 2400 statements and CTRs, annotated for these
tasks. Baselines on this corpus expose the limitations of existing NLI models,
with 6 state-of-the-art NLI models achieving a maximum F1 score of 0.627. To
the best of our knowledge, we are the first to design a task that covers the
interpretation of full CTRs. To encourage further work on this challenging
dataset, we make the corpus, competition leaderboard, website, and code to
replicate the baseline experiments available at: https://github.com/ai-
systems/nli4ct

##  1 Introduction

Clinical trials are research studies performed to test the efficacy and safety
of novel treatments, and they are indispensable for the progression of
experimental medicine Avis et al. (2006). CTRs are documents that detail the
methodology and results of a particular trial. Clinical practitioners use the
information in these CTRs to design personalized and targeted interventions,
matching therapeutic agents to each patient’s tumor biomarkers. However, there
are over 400,000 CTRs available, being continually published at an ever-
increasing rate Bastian et al. (2010). Consequently, it has become unfeasible
to manually carry out comprehensive evaluations of all the relevant literature
when designing treatment protocols DeYoung et al. (2020). Natural Language
Inference (NLI) offers a potential solution for the large-scale interpretation
and retrieval of medical evidence, connecting the latest evidence to support
personalized care Sutton et al. (2020).

In this paper, we propose two different tasks, on breast cancer CTRs. Firstly
to determine the inference relation between a natural language statement, and
a CTR, shown in Figure 1. Secondly, to retrieve supporting facts from the
CTR(s) to justify the predicted relation demonstrated in Figure 1. With
respect to healthcare, these inference tasks are a move towards automatic
claim verification and evidence extraction, over biomedical corpora.

Figure 1: We propose two tasks for reasoning on clinical trial data expressed
in natural language. Firstly, to predict the entailment of a Statement and a
CTR premise, and secondly, to extract evidence to support the label.

These proposed tasks represent several challenges for existing NLI systems.
The inference task requires substantial amounts of quantitative reasoning and
numerical operations, as seen in Figure 1, on which many current state-of-the-
art NLI models perform poorly Ravichander et al. (2019); Galashov et al.
(2019). Additionally, much of the inference for these tasks revolves around
domain-specific terminology, also shown in Figure 1. Many state-of-the-art
(SOTA) NLI models fail to effectively surmount the word distribution shift
from general domain corpora to biomedical corpora through transfer learning
Lee et al. (2019). This phenomenon is well-represented in the proposed
dataset. This task differentiates itself from other biomedical NLI tasks by
considering the full scope of CTRs, not limited to results and interventions,
instead including eligibility and Adverse events (AE) as well DeYoung et al.
(2020). This task also provides a significant variance in the types of
inference required to solve each instance, with minimal, if any, repetition in
inference chains.

The contributions of this paper are as follows:

  1. 1.

The definition of a novel benchmark (NLI4CT), including two main tasks
requiring reasoning over CTRs, incorporates several of the fundamental
challenges currently faced by modern NLI systems over clinical trial text.

  2. 2.

Publishing a newly publicly-available corpus containing 2,400 expert-annotated
entailment relations, with associated CTRs from
ClinicalTrials.gov111https://clinicaltrials.gov/ct2/home, labels, and
extracted evidence lists.

  3. 3.

An extensive empirical evaluation of state-of-the-art NLI models demonstrating
current limitations and challenges involved in the proposed tasks.
Specifically, we test 6 SOTA NLI models on our dataset and report a maximum F1
score of 0.644. We attribute these results to generalization problems induced
by the shift in distribution required for inference.

To the best of our knowledge, we are the first to design a task that covers
the full scope of CTRs, combining together the complexity of biomedical and
numerical NLI. Solving these tasks would be a significant advancement toward
efficient retrieval, synthesis, and inference of published literature to
support evidence-based experimental medicine.

##  2 Related Work

Table 1: Table of Clinical NLP datasets.

There is a multitude of expert-annotated resources for clinical NLP, examples
are shown in Table 1. The TREC 2021 Clinical Track Soboroff (2021) is a large-
scale information retrieval task on CTR data, with an emphasis on eligibility.
Evidence Inference 2.0 DeYoung et al. (2020) presents a Question-Answering
(QA) task and a span selection task on CTR data, specifically with the CTR
results. The MEDNLI Romanov and Shivade (2018) dataset contains an entailment
task, which utilizes the medical history notes of patients as premises. These
datasets are predominantly designed to test for biomedical language
understanding and reasoning. Whereas NLI4CT, additionally requires complex
numerical inference and common-sense reasoning to solve and is designed to
avoid repetitive inference patterns. Currently, neural architectures achieve
the best results on biomedical NLI datasets Gu et al. (2021); DeYoung et al.
(2020).

Neural models continue to perform poorly on quantitative reasoning and
numerical operations within NLI Ravichander et al. (2019); Galashov et al.
(2019). Prior works such as Lee et al. (2020), Shin et al. (2020), and Gu et
al. (2021) experiment with biomedical pre-training strategies. Kiritchenko et
al. (2010) presents ExaCT, an automatic information extraction system for
clinical trial variables such as drug dosage, sample size, and primary
outcomes from CTRS. While there are many proficient evidence retrieval models,
there are no existing systems that are able to effectively carry out
biomedical and numerical NLI simultaneously.

##  3 Clinical Trial Reports

Clinical trials are research studies on human volunteers used to evaluate an
intervention DeYoung et al. (2020). An intervention may be medical, surgical,
or a behavioral protocol. These interventions are often compared with a
standard treatment, placebo, or control protocol. The investigators use
various outcome measurements from the volunteers to ascertain the
effectiveness of the intervention, as well as the likelihood and severity of
AEs.

For NLI4CT we retrieved 1000 publicly available Breast cancer CTRs with
published results from ClinicalTrials.gov. This data is developed by the U.S.
National Library of Medicine and covered by the HIPAA Privacy Rule, protecting
personally identifiable information. We separate the CTRs into 4 sections:

  * •

Eligibility criteria \- A set of conditions for patients to be allowed to take
part in the clinical trial

  * •

Intervention \- Information concerning the type, dosage, frequency, and
duration of treatments being studied.

  * •

Results \- Number of participants in the trial, outcome measures, units, and
the results.

  * •

Adverse Events \- These are signs and symptoms observed in patients during the
clinical trial.

See examples of the types of information contained in these sections in Table
2.

Eligibility
Criteria |  Inclusion:
• HER2 + diagnosis
• 18 years of age or older
Exclusion:
• Presence of central nervous
system or brain metastases.
---|---
Intervention |  Gemcitabine 1500 mg/m2 body surface
area (BSA) intravenously (IV) over 30
minutes (+/- 5 minutes) on days 1 and 15 of each cycle.
Results |  Variable: Progression Free Survival
Number of Participants: 29
Median (months): 10.4 (5.6 to 15.2)
Adverse Events |  Total: 3/29
• Neutropenic Fever 1 / 29 (3.45%)
• Peripheral Neuropathy 1 / 29 (3.45%)
• Seizure/Syncope 1 / 29 (3.45%)
Table 2: Sample excerpts from each section in a CTR. Data provided by
ClinicalTrials.gov.

##  4 Task Definition

We define two different tasks on the NLI4CT dataset, Task 1, a textual
entailment task, and Task 2, an evidence selection task. Each instance in
NLI4CT contains a CTR premise and a statement. The CTR premise consists of one
of the 4 sections of a CTR, with a length ranging from 5-500 tokens, and the
statements are sentences ranging from 10-35 tokens (see example in Figure 1).
On average there are 7.74 pieces of relevant evidence that need to be selected
out of a total of 21.67 facts per CTR premise. There are two types of
instances in NLI4CT, single and comparison, in single instances the CTR
premise is the primary CTR, and the statement makes a claim about the
information in this CTR. Whereas in the comparison type the premise contains a
primary and a secondary trial, and the statement will make a claim comparing
and contrasting the two trials. To summarize:

##### Task 1

Requires determining the entailment relation between the CTR premise and the
statement, therefore the output is either an entailment or contradiction
label, as shown in Figure 1.

##### Task 2

Requires outputting a set of supporting facts, extracted from the CTR premise,
necessary to justify the label predicted in Task 1.

##  5 Annotation

Figure 2: Schematic of the NLI4CT annotation process.

###  5.1 Statement annotation

A group of domain experts, including clinical trial organizers from a major
cancer research center, took part in the annotation task. Annotators were
provided with a CTR section prompt, and two CTRs, a primary and a secondary,
as seen in step (B) of Figure 2. The annotation task was first to generate an
entailment statement, a short statement that makes an objectively true claim
about the contents of the prompted section of the trial(s), Step (C) of Figure
2 Annotators could choose to write a statement about the contents of the
primary trial or to compare both the primary and secondary trials. The
objective was to generate non-trivial statements, i.e. solving the inference
requires reading and understanding multiple rows of the CTR section. Non-
trivial statements typically included summarisation, comparisons, negation,
relations, inclusion, superlatives, aggregations, or rephrasing. However, non-
trivial sentences are not limited to these types, and any statement involving
understanding and reasoning was accepted.

Each CTR is divided into 4 sections (Table 1), and each of these sections
represents a set of facts. And for all annotated statements we select a subset
of these facts that support the label assigned to the statement, see Figure 1
and (D) in Figure 2. This collection of facts is designated as evidence. There
is always at least one piece of evidence relevant to the inference for any
given statement. In cases where negation is used, e.g. clinical trial A does
not have a placebo arm annotators are asked to provide the full CTR section as
evidence, as we believe this type of retrieval is more reflective of human
inference patterns.

Finally, we employ a negative rewriting strategy Chen et al. (2019), using the
previously generated entailment statement to write an objectively wrong
contradiction statement, as shown in Step (E) of Figure 2. Annotators are
encouraged to modify the words, phrases, or sentence structures while
retaining the sentence style/length of the original statement. This is done to
minimize the occurrence of stylistic or linguistic patterns pertaining to
either entailment or contradictory statements. We then collect evidence that
contradicts the claims made in the contradiction statement. An example
annotation form is available in the appendix.

###  5.2 Resulting Dataset

The NLI4CT dataset consists of 2,400 annotated statements with accompanying
labels, CTRs, and evidence. Split into 1700 training, 500 test, and 200
development instances. The two labels and 4 CTR sections prompts are equally
distributed across the dataset and its splits. $60\%$ of the instances are the
single type, with the remaining $40\%$ the comparison type.

##  6 Reasoning Challenges

###  6.1 Biomedical Reasoning

Acronyms In CTRs acronyms are used for concepts such as test names,
administration times and orders, treatment types, and disease names. Acronyms
are significantly more prevalent in clinical texts than in general-domain
texts and consistently disrupt NLP performance in this field Grossman Liu et
al. (2021); Shickel et al. (2017); Jiang et al. (2011); Moon et al. (2015);
Jimeno-Yepes et al. (2011); Pesaranghader et al. (2019); Jin et al. (2019); Wu
et al. (2015). Consider the following statement, "Patients with a positive
FISH test will be administered trastuzumab IV Q.D.". Depending on the
terminology in the premise, verifying all the claims in this statement
requires a system to know the association between FISH/Fluorescence In Situ
Hybridization, IV/intravenous infusion, and QD/once daily. Over 100,000
healthcare acronyms and abbreviations have been identified with 170,00
corresponding senses Grossman Liu et al. (2021). Future systems must overcome
this challenge to perform biomedical NLI.

Synonyms and aliases These are often employed for drugs and gene names.
Statements in NLI4CT regularly refer to treatments and diagnoses using
different aliases than are present in the CTR premise, e.g. where a CTR may
use Trastuzumab in its intervention, the statement may refer to Herceptin or
any of the 4 other brand names of Trastuzumab. Similarly, where the
eligibility criteria may require an ERBB2+ tumor diagnosis from its patients,
the statement may refer to HER2 positive or any of the 20 other aliases for
the ERBB2 gene.

Taxonomic Relations Concepts such as diseases, treatments, and diagnostic
tests can be classified and structured in a taxonomic hierarchy. As shown in
Figure 1, to achieve the necessary inference systems must understand that
chemotherapy is a hypernym or a super-set of Cyclophosphamide treatments,
similar to Cancer being a hypernym of Carcinoma.

Domain Knowledge NLI4CT statements regularly test for domain expert knowledge.
This is done by describing characteristics or conditions which pertain to
certain biomedical grading systems or diagnoses. For example, if a patient is
confined to bed or a chair for more than 50% of waking hours, and the
inclusion criteria require a WHO score of < 2, the model must be able to infer
that the patient will not be eligible. Or if a patient has a positive FISH
test, and the inclusion criteria require the patient to have a Triple Negative
Breast Cancer diagnosis, they would also be ineligible.

###  6.2 Common sense Reasoning

NLI4CT requires common sense reasoning, particularly co-reference resolution
and general world knowledge, which remain challenging for existing NLI models
Emami et al. (2018); Trichelair et al. (2018). Co-referencing resolution is
often necessary to associate dosages, frequencies, and routes of
administration to specific drugs in the intervention. General world knowledge
is most often required for claims about the eligibility section, e.g. a child
is a person < 18 years old, and women over the age of 60 are not of
childbearing potential.

###  6.3 Numerical Reasoning

The type of numerical reasoning required to solve instances most closely
resembles Math Word Problems (MWP) Huang et al. (2016); Patel et al. (2021).
Prior works have shown that even elementary MWPs remain unsolved for current
NLP models Patel et al. (2021). NLI4CT requires models to compare dosages,
frequencies, and fractions/percentages often necessitating unit conversions.
Consider the example instance in Figure 3. Here the model will need to
aggregate the number of occurrences of every different cardiac event recorded
in the adverse events of cohort 1. Then it must construct and evaluate the
inequality. This combination of biomedical knowledge and numerical reasoning
is expected to be a significant challenge for existing models.

Figure 3: Schematic of a numerical inference chain to solve an NLI4CT
instance.

##  7 Experiments

###  7.1 Baseline Model Architectures

We establish a baseline performance on NLI4CT, to assess the strengths and
limitations of existing NLI models. We test a range of transformer-based
architectures; BERT-base Devlin et al. (2018), RoBERTa-base Liu et al. (2019),
GPT-2 Radford et al. (2019), T5-base Raffel et al. (2020), and Megatron-lm
Shoeybi et al. (2019). Additionally, we experiment with the bag-of-words
ranking function Okapi BM25 Fang et al. (2004); Robertson et al. (1995).

###  7.2 Biomedical Pretraining

Transformer-based models leverage pre-trained contextualized word embeddings
to model inference, and optimizing pre-training strategies can drastically
improve downstream performance Gururangan et al. (2020). Previous works have
repeatedly shown that pre-training general-domain language models on
biomedical corpora significantly improves performance on bio-medical tasks
Shin et al. (2020); Beltagy et al. (2019); Lee et al. (2020); Li et al.
(2016). Therefore, we also include BioBERT Lee et al. (2020), and BioMegatron
Shin et al. (2020) (versions of BERT-base and Megatron-lm that are pre-trained
on biomedical corpora) in our baselines to assess the effects of pre-training
on performance on our dataset.

###  7.3 Experimental Details

For the transformer-based model baselines we tokenize the statements and the
section(s) of the CTR(s) as indicated by the instance type and section prompt.
If the instance type is comparison we concatenate the primary and secondary
CTR section data together. Then we pass the tokenized statement and section
information to the model. It should be noted that for a significant portion of
the instances, particularly of the comparison type, the input exceeds the
maximum input length of 512 tokens. We employ a binary sequence
classification/regression head for all baseline models, which predicts either
Entailement or Contradiction. The models were trained on the training set for
10 epochs, and each model checkpoint was archived. For the Okapi BM25
similarity, we test threshold values from 0.1 - 0.9, in intervals of 0.1 on
the development set. We then record the results of these models on the test
set. Code to reproduce the experiments is available at: anonymous-link

##  8 Evaluation Metrics

Task 1 is a binary classification task, the statement is either labeled as
entailment or contradiction. We evaluate performance on task 1 by computing
the Precision, Recall, and Macro F1-score of the predicted labels against the
annotated gold labels (Entailment/Contradiction). For Task 2 the output is a
subset of the facts within the CTR premise. We can frame the evidence
selection problem as a ranking problem Jansen et al. (2021). Models assign a
score to each fact in a CTR premise, which can be used to rank the facts. In
this framework, solving the task requires the model to rank all the facts in
the gold evidence higher than the irrelevant facts. We evaluate performance on
this task using Mean Average Precision mAP @ K Liu (2011).

##  9 Experiment Results

###  9.1 Task 1

Model | F1 | Precision | Recall | Accuracy
---|---|---|---|---
BM25 | 0.627 | 0.490 | 0.872 | 0.482
BERT-base | 0.528 | 0.555 | 0.504 | 0.550
RoBERTa-base | 0.628 | 0.654 | 0.604 | 0.642
BioBERT | 0.644 | 0.633 | 0.656 | 0.638
BioMegatron | 0.644 | 0.585 | 0.716 | 0.604
GPT2 | 0.603 | 0.500 | 0.760 | 0.500
T5-base | 0.602 | 0.504 | 0.748 | 0.506

Table 3: Results from the NLI4CT Task 1 baseline on the test set.

Table 8 summarises the performance of the baseline models on Task 1. We report
the results from the model checkpoints with the highest F1 score on the test
set. However, we omit the results of models which exclusively predict a single
class. None of the models achieve an F1-score significantly above the random
baseline (0.667). BM25, T5, BERT-base and GPT2 also fail achieve over 0.5
accuracy on the test set. The best-performing models on the test set were
BioBERT, BioMegatron, and RoBERTa-base, each achieving over 0.6 accuracy and
F1 score. BERT base reported a significantly lower F1 than the other
baselines.

BioBERT, and BioMegatron both achieve over 0.95 accuracy and F1 score on the
training set after 10 epochs. As shown in Table 8, this does not translate to
an improvement in performance on the test, indicating drastic over-fitting.
BioBERT significantly outperforms its general-domain counterpart on the
training set, achieving a maximum F1 score of 0.96, in contrast to 0.65 by
BERT-base. Biomedical pre-training clearly impacts the models’ ability to
encode relevant information, limited to the training set. On average, model
complexity does not have a consistent impact on performance. Neither the over-
fitting biomedical models nor the models which failed to learn on the training
set exhibited any positive responses to changes in the learning rate. The
aggregation of these results indicates that the NLI4CT entailment task
represents a significant challenge for existing NLI models.

###  9.2 Evidence-only Baseline

We generated an additional set of baselines, using only the gold evidence as
the premise. This reduces the likelihood of the models focusing on irrelevant
or adversarial information and removes the need to identify evidence spans
within the CTR section. This did not result in any significant differences in
performance, with BioBERT achieving the highest F1 score of 0.667. This
implies that even after training for 10 epochs none of the models draw
inferences from the relevant evidence. The results for this baseline are
available in the appendix.

###  9.3 Statistical Artefacts

We generated a baseline with BERT-base, using only the statements, removing
any access to evidence or context from the CTRs. BERT-base achieved a maximum
F1 score of 0.626 and an accuracy of 0.606 on the test set. This is a
significantly higher F1 score than on the Task 1 baseline. Additionally, we
observed above-random accuracy (0.5) with several model checkpoints. This
indicates that the evaluated models might exclusively rely on the presence of
superficial statistical artefacts without learning the underlying rules of the
tasks.

Lexical and syntactic patterns such as token distributions, statement lengths,
and discriminative conditions that are disproportionately associated with a
particular class of a dataset can superficially inflate model performance
Herlihy and Rudinger (2021). These statistical artifacts are typically
introduced during the annotation Gururangan et al. (2018) and allow statement-
only classifiers to achieve significantly above-random results Poliak et al.
(2018); Tsuchiya (2018).

We support the statement-only baseline with qualitative analysis of the
dataset. We found that there was no significant difference between the
sentence length distributions of the two classes. Additionally, we find a 100%
overlap in the 15 most common tokens between the two classes in the training
set, and 93% overlap in the test set. However, several individual tokens were
not evenly distributed across classes in the test set. In particular the
distributions of not, more, group, and cohort. But these uneven distributions
were not present, if not completely reversed in the training set. Therefore
these do not entirely explain the above random accuracy observed.

###  9.4 Categorical Analysis

Numerical Reasoning We categorize the instances in our dataset into those that
require numerical or quantitative reasoning and those that require pure NLI.
We separate out all instances with statements that contain numbers or any of
the following tokens; higher, lower, fewer, more, less, number, same into the
Numerical reasoning category. This results in a 62% Numerical 38% NLI split of
the test set. We evaluate the best-performing checkpoints on these categories
and report the results in Figure 4. BioBERT reported the highest F1 score for
both Numerical and NLI categories, at 0.552 and 0.660 respectively. The
performance on the NLI instances is significantly higher than the performance
on the numerical instances, for all models other than BM25 and GPT2 which only
marginally report the same trend. Biomedical pretraining does not appear to
provide a significant advantage within either category. Overall these results
follow the trends observed in prior works, numerical reasoning remains a
challenging task for NLI models Ravichander et al. (2019); Galashov et al.
(2019); Peng et al. (2021); Patel et al. (2021).

Figure 4: Graph of the baseline model F1 scores on the NLI vs Numerical
instances in the test set.

Single vs Comparison We repeat this experiment, this time categorizing by
instance type. Comparison type instances have 2 CTR sections as a premise,
effectively doubling the amount of irrelevant and potentially adversarial
information. Therefore, we would expect comparison instances to be more
challenging than the single type. We report the results of this experiment in
Figure 5. Contrary to expectations the majority of the models there was no
significant difference in performance across the categories. BioMegatron
reports a +0.1 F1 on single type instances, and BM25, and T5 report +0.09 and
+0.07 F1 score respectively on the comparison type. As T5 failed to learn Task
1, even on the training set, and BM25 is not able to capture clues that go
beyond lexical overlaps, we hypothesize that this is in part due to the model
relying on superficial clues to predict the output label. Therefore, in many
cases, the evidence is just not being leveraged properly.

CTR Sections Finally, we categorize the instances in our dataset by CTR
section and report the baselines on these categories in Table 4. All models
show significant fluctuations in performance across the different sections.
Interestingly, Results and AEs are on opposite ends of the rankings, despite
arguably containing the most similar types of information. Additionally,
Eligibility is ranked 2nd despite having the highest average number of tokens.

Figure 5: Graph of the baseline model F1 scores on the Single vs Comparison
instances in the test set.

| F1 score
---|---
Model | AEs | Eligibility | Intervention | Results
BERT-base | 0.649 | 0.554 | 0.547 | 0.537
BioBERT | 0.619 | 0.667 | 0.588 | 0.628
BioMegatron | 0.661 | 0.576 | 0.626 | 0.528
BM25 | 0.158 | 0.490 | 0.529 | 0.574
GPT2 | 0.615 | 0.531 | 0.523 | 0.598
T5-base | 0.504 | 0.515 | 0.500 | 0.455
RoBERTa | 0.667 | 0.609 | 0.523 | 0.460
Average Rank | 1.86 | 2.29 | 2.86 | 3.00

Table 4: F1 score of the baseline models on the test set, categorized by CTR
section.

###  9.5 Task 2

For Task 2, the evidence selection task, we generate a baseline using SOTA
general purpose transformer-based ranking models Reimers (2020), namely,
DistilRoBERTa Sanh et al. (2019), Mpnet Song et al. (2020), and MiniLM Wang et
al. (2020) variants. We additionally test BM25. We embed the facts in the CTR
sections, and the statements using these pre-trained models and then compute
the cosine similarity. We evaluate using the mAP @ K Liu (2011), with K =
total number of sentences. We report the results in Table 5. There is no
significant difference between the mAP scores of the models, with BM25
producing the highest mAP score. All models achieve well above the random
baseline. This indicates that the limiting factor for SOTA models on NLI4CT is
not evidence selection, but rather biomedical, and numerical inference.

Model | mAP
---|---
Random Baseline | 0.673
all-MiniLM-L6-v2 | 0.777
all-mpnet-base-v2 | 0.760
all-distilroberta-v1 | 0.762
multi-qa-MiniLM-L6-cos-v1 | 0.777
multi-qa-mpnet-base-dot-v1 | 0.770
BM25 | 0.786
Table 5: Results from the NLI4CT Task 2 baseline on the test set.

##  10 Conclusions & Future work

We present 2 tasks, textual entailment of claims over clinical trial reports
(CTRs), and extracting evidence from CTRs to support or contradict these
claims. We provide a corpus of 2400 expert-annotated instances to test and
train models on these tasks. If models could be developed to solve these
tasks, it would progress the field towards the efficient retrieval and
synthesis of published literature to support evidence-based medicine.
Additionally, it would address several important NLU challenges related to
numerical reasoning, biomedical NLI, and inference over long texts.

Our baselines outline the weaknesses of existing NLI models, particularly with
regard to numerical inference. We tested 6 SOTA NLI models on Task 1 and
reported a maximum F1 score of 0.644. Additionally, we show that the limiting
factor for current models is not evidence selection, but rather leveraging
that evidence, and using numerical and biomedical inference for predictions.
Our corpus, competition leaderboard, website, and code to replicate the
baseline experiments are available at: anonymous-link.

##  11 Limitations

The most significant limitation of this work is that despite the fact that we
are dialoguing with a motivational scenario for medicine, these models are not
fit for medical application, this would require a technology clinical trial
and regulatory assessment.

Additionally, our evaluation does not include an interventional study, meaning
we do not perform perturbations/interventions on the test instances to verify
that the models learn the underlying causal structure of the tasks. We plan to
do this in future work.

Finally, our dataset contains fewer instances than other published biomedical
NLI datasets. This is a consequence of the time and resource-intensive nature
of expert-level annotation. The size of the training set may be problematic
for large models.

## References

  * Avis et al. (2006) Nancy E Avis, Kevin W Smith, Carol L Link, Gabriel N Hortobagyi, and Edgardo Rivera. 2006.  Factors associated with participation in breast cancer treatment clinical trials.  _J Clin Oncol_ , 24(12):1860–1867.
  * Bastian et al. (2010) Hilda Bastian, Paul Glasziou, and Iain Chalmers. 2010.  Seventy-five trials and eleven systematic reviews a day: how will we ever keep up?  _PLoS medicine_ , 7(9):e1000326.
  * Beltagy et al. (2019) Iz Beltagy, Kyle Lo, and Arman Cohan. 2019.  Scibert: A pretrained language model for scientific text.  _arXiv preprint arXiv:1903.10676_.
  * Chen et al. (2019) Wenhu Chen, Hongmin Wang, Jianshu Chen, Yunkai Zhang, Hong Wang, Shiyang Li, Xiyou Zhou, and William Yang Wang. 2019.  Tabfact: A large-scale dataset for table-based fact verification.  _CoRR_ , abs/1909.02164.
  * 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_.
  * DeYoung et al. (2020) Jay DeYoung, Eric P. Lehman, Benjamin E. Nye, Iain James Marshall, and Byron C. Wallace. 2020.  Evidence inference 2.0: More data, better models.  _ArXiv_ , abs/2005.04177.
  * Emami et al. (2018) Ali Emami, Paul Trichelair, Adam Trischler, Kaheer Suleman, Hannes Schulz, and Jackie Chi Kit Cheung. 2018.  The knowref coreference corpus: Removing gender and number cues for difficult pronominal anaphora resolution.  _arXiv preprint arXiv:1811.01747_.
  * Fang et al. (2004) Hui Fang, Tao Tao, and ChengXiang Zhai. 2004.  A formal study of information retrieval heuristics.  In _Proceedings of the 27th annual international ACM SIGIR conference on Research and development in information retrieval_ , pages 49–56.
  * Galashov et al. (2019) Alexandre Galashov, Jonathan Schwarz, Hyunjik Kim, Marta Garnelo, David Saxton, Pushmeet Kohli, S. M. Ali Eslami, and Yee Whye Teh. 2019.  Meta-learning surrogate models for sequential decision making.  _CoRR_ , abs/1903.11907.
  * Grossman Liu et al. (2021) Lisa Grossman Liu, Raymond H Grossman, Elliot G Mitchell, Chunhua Weng, Karthik Natarajan, George Hripcsak, and David K Vawdrey. 2021.  A deep database of medical abbreviations and acronyms for natural language processing.  _Scientific Data_ , 8(1):1–9.
  * Gu et al. (2021) Yu Gu, Robert Tinn, Hao Cheng, Michael Lucas, Naoto Usuyama, Xiaodong Liu, Tristan Naumann, Jianfeng Gao, and Hoifung Poon. 2021.  Domain-specific language model pretraining for biomedical natural language processing.  _ACM Transactions on Computing for Healthcare (HEALTH)_ , 3(1):1–23.
  * Gururangan et al. (2020) Suchin Gururangan, Ana Marasović, Swabha Swayamdipta, Kyle Lo, Iz Beltagy, Doug Downey, and Noah A Smith. 2020.  Don’t stop pretraining: adapt language models to domains and tasks.  _arXiv preprint arXiv:2004.10964_.
  * Gururangan et al. (2018) Suchin Gururangan, Swabha Swayamdipta, Omer Levy, Roy Schwartz, Samuel R Bowman, and Noah A Smith. 2018.  Annotation artifacts in natural language inference data.  _arXiv preprint arXiv:1803.02324_.
  * Herlihy and Rudinger (2021) Christine Herlihy and Rachel Rudinger. 2021.  Mednli is not immune: Natural language inference artifacts in the clinical domain.  _arXiv preprint arXiv:2106.01491_.
  * Huang et al. (2016) Danqing Huang, Shuming Shi, Chin-Yew Lin, Jian Yin, and Wei-Ying Ma. 2016.  How well do computers solve math word problems? large-scale dataset construction and evaluation.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 887–896.
  * Jansen et al. (2021) Peter Jansen, Mokanarangan Thayaparan, Marco Valentino, and Dmitry Ustalov. 2021\.  Textgraphs 2021 shared task on multi-hop inference for explanation regeneration.  In _Proceedings of the Fifteenth Workshop on Graph-Based Methods for Natural Language Processing (TextGraphs-15)_.
  * Jiang et al. (2011) Min Jiang, Yukun Chen, Mei Liu, S Trent Rosenbloom, Subramani Mani, Joshua C Denny, and Hua Xu. 2011.  A study of machine-learning-based approaches to extract clinical entities and their assertions from discharge summaries.  _Journal of the American Medical Informatics Association_ , 18(5):601–606.
  * Jimeno-Yepes et al. (2011) Antonio J Jimeno-Yepes, Bridget T McInnes, and Alan R Aronson. 2011.  Exploiting mesh indexing in medline to generate a data set for word sense disambiguation.  _BMC bioinformatics_ , 12(1):1–14.
  * Jin et al. (2019) Qiao Jin, Jinling Liu, and Xinghua Lu. 2019.  Deep contextualized biomedical abbreviation expansion.  _arXiv preprint arXiv:1906.03360_.
  * Kiritchenko et al. (2010) Svetlana Kiritchenko, Berry De Bruijn, Simona Carini, Joel Martin, and Ida Sim. 2010\.  Exact: automatic extraction of clinical trial characteristics from journal publications.  _BMC medical informatics and decision making_ , 10(1):1–17.
  * Lee et al. (2019) Jinhyuk Lee, Wonjin Yoon, Sungdong Kim, Donghyeon Kim, Sunkyu Kim, Chan Ho So, and Jaewoo Kang. 2019.  Biobert: a pre-trained biomedical language representation model for biomedical text mining.  _CoRR_ , abs/1901.08746.
  * Lee et al. (2020) Jinhyuk Lee, Wonjin Yoon, Sungdong Kim, Donghyeon Kim, Sunkyu Kim, Chan Ho So, and Jaewoo Kang. 2020.  Biobert: a pre-trained biomedical language representation model for biomedical text mining.  _Bioinformatics_ , 36(4):1234–1240.
  * Li et al. (2016) Jiao Li, Yueping Sun, Robin J Johnson, Daniela Sciaky, Chih-Hsuan Wei, Robert Leaman, Allan Peter Davis, Carolyn J Mattingly, Thomas C Wiegers, and Zhiyong Lu. 2016.  Biocreative v cdr task corpus: a resource for chemical disease relation extraction.  _Database_ , 2016.
  * Liu (2011) Tie-Yan Liu. 2011.  Applications of learning to rank.  _Learning to Rank for Information Retrieval_ , page 181–191.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.  _arXiv preprint arXiv:1907.11692_.
  * Moon et al. (2015) Sungrim Moon, Bridget McInnes, and Genevieve B Melton. 2015.  Challenges and practical approaches with word sense disambiguation of acronyms and abbreviations in the clinical domain.  _Healthcare informatics research_ , 21(1):35–42.
  * Patel et al. (2021) Arkil Patel, Satwik Bhattamishra, and Navin Goyal. 2021.  Are nlp models really able to solve simple math word problems?  _arXiv preprint arXiv:2103.07191_.
  * Peng et al. (2021) Shuai Peng, Ke Yuan, Liangcai Gao, and Zhi Tang. 2021.  Mathbert: A pre-trained model for mathematical formula understanding.  _arXiv preprint arXiv:2105.00377_.
  * Pesaranghader et al. (2019) Ahmad Pesaranghader, Stan Matwin, Marina Sokolova, and Ali Pesaranghader. 2019.  deepbiowsd: effective deep neural word sense disambiguation of biomedical text data.  _Journal of the American Medical Informatics Association_ , 26(5):438–446.
  * Poliak et al. (2018) Adam Poliak, Jason Naradowsky, Aparajita Haldar, Rachel Rudinger, and Benjamin Van Durme. 2018.  Hypothesis only baselines in natural language inference.  _arXiv preprint arXiv:1805.01042_.
  * Radford et al. (2019) Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, Ilya Sutskever, et al. 2019.  Language models are unsupervised multitask learners.  _OpenAI blog_ , 1(8):9.
  * Raffel et al. (2020) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, Peter J Liu, et al. 2020.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _J. Mach. Learn. Res._ , 21(140):1–67.
  * Ravichander et al. (2019) Abhilasha Ravichander, Aakanksha Naik, Carolyn Penstein Rosé, and Eduard H. Hovy. 2019.  EQUATE: A benchmark evaluation framework for quantitative reasoning in natural language inference.  _CoRR_ , abs/1901.03735.
  * Reimers (2020) Nils Reimers. 2020.  Pretrained models.  https://www.sbert.net/docs/pretrained_models.html.  Accessed: 2010-09-30.
  * Robertson et al. (1995) Stephen E Robertson, Steve Walker, Susan Jones, Micheline M Hancock-Beaulieu, Mike Gatford, et al. 1995.  Okapi at trec-3.  _Nist Special Publication Sp_ , 109:109.
  * Romanov and Shivade (2018) Alexey Romanov and Chaitanya Shivade. 2018.  Lessons from natural language inference in the clinical domain.  _arXiv preprint arXiv:1808.06752_.
  * Sanh et al. (2019) Victor Sanh, Lysandre Debut, Julien Chaumond, and Thomas Wolf. 2019.  Distilbert, a distilled version of BERT: smaller, faster, cheaper and lighter.  _CoRR_ , abs/1910.01108.
  * Shickel et al. (2017) Benjamin Shickel, Patrick James Tighe, Azra Bihorac, and Parisa Rashidi. 2017.  Deep ehr: a survey of recent advances in deep learning techniques for electronic health record (ehr) analysis.  _IEEE journal of biomedical and health informatics_ , 22(5):1589–1604.
  * Shin et al. (2020) Hoo-Chang Shin, Yang Zhang, Evelina Bakhturina, Raul Puri, Mostofa Patwary, Mohammad Shoeybi, and Raghav Mani. 2020.  Biomegatron: Larger biomedical domain language model.  _arXiv preprint arXiv:2010.06060_.
  * Shoeybi et al. (2019) Mohammad Shoeybi, Mostofa Patwary, Raul Puri, Patrick LeGresley, Jared Casper, and Bryan Catanzaro. 2019.  Megatron-lm: Training multi-billion parameter language models using model parallelism.  _arXiv preprint arXiv:1909.08053_.
  * Soboroff (2021) Ian Soboroff. 2021.  Overview of trec 2021.  In _30th Text REtrieval Conference. Gaithersburg, Maryland_.
  * Song et al. (2020) Kaitao Song, Xu Tan, Tao Qin, Jianfeng Lu, and Tie-Yan Liu. 2020.  Mpnet: Masked and permuted pre-training for language understanding.  _Advances in Neural Information Processing Systems_ , 33:16857–16867.
  * Sutton et al. (2020) Reed T Sutton, David Pincock, Daniel C Baumgart, Daniel C Sadowski, Richard N Fedorak, and Karen I Kroeker. 2020.  An overview of clinical decision support systems: benefits, risks, and strategies for success.  _NPJ digital medicine_ , 3(1):1–10.
  * Trichelair et al. (2018) Paul Trichelair, Ali Emami, Adam Trischler, Kaheer Suleman, and Jackie Chi Kit Cheung. 2018.  How reasonable are common-sense reasoning tasks: A case-study on the winograd schema challenge and swag.  _arXiv preprint arXiv:1811.01778_.
  * Tsuchiya (2018) Masatoshi Tsuchiya. 2018.  Performance impact caused by hidden bias of training data for recognizing textual entailment.  _arXiv preprint arXiv:1804.08117_.
  * Wang et al. (2020) Wenhui Wang, Furu Wei, Li Dong, Hangbo Bao, Nan Yang, and Ming Zhou. 2020.  Minilm: Deep self-attention distillation for task-agnostic compression of pre-trained transformers.  _CoRR_ , abs/2002.10957.
  * Wu et al. (2015) Yonghui Wu, Jun Xu, Yaoyun Zhang, and Hua Xu. 2015.  Clinical abbreviation disambiguation using neural word embeddings.  In _Proceedings of BioNLP 15_ , pages 171–176.

###  11.1 Appendices

####  11.1.1 Pre-processing details

We extracted all CTR data from ClinicalTrials.gov on the 22${}^{n}d$ of
December 2022. The Eligibility criteria and intervention information is
retrieved from the Eligibility criteria and Intervention sections. The results
information is taken from the Primary Outcome section. The Adverse events
section is taken from the Serious Adverse Events section, not to be confused
with the adverse events section. We remove or replace symbols or tags produced
in the extraction, and segment the CTR data into individual lines. These are
the only changes made to the CTR data.

###  11.2 Tables and Graphs

| F1 score
---|---
Model | Single | Comparison
BERT-base | $0.564$ | $0.582$
BioBERT | $0.635$ | $0.619$
BioMegatron | $0.648$ | $0.550$
BM25 | $0.517$ | $0.431$
GPT-2 | $0.585$ | $0.550$
RoBERTa-base | $0.569$ | $0.576$
T5-base | $0.455$ | $0.527$
Table 6: F1 score of the baseline models on the test set, categorized by instance type. | F1 score
---|---
Model | Numerical | NLI
BERT-base | 0.493 | 0.664
BioBERT | 0.588 | 0.686
BioMegatron | 0.552 | 0.660
BM25 | 0.468 | 0.495
GPT-2 | 0.553 | 0.577
RoBERTa-base | 0.463 | 0.676
T5-base | 0.429 | 0.596
Table 7: F1 score of the baseline models on the test set, categorized by
numerical and NLI inference.

Model | F1 | Precision | Recall | Accuracy
---|---|---|---|---
BM25 | 0.622 | 0.489 | 0.856 | 0.480
BERT-base | 0.597 | 0.569 | 0.628 | 0.576
RoBERTa-base | 0.659 | 0.506 | 0.944 | 0.512
BioBERT | 0.667 | 0.589 | 0.768 | 0.616
BioMegatron | 0.596 | 0.543 | 0.660 | 0.552
GPT2 | 0.581 | 0.519 | 0.660 | 0.524
T5-base | 0.521 | 0.522 | 0.520 | 0.522

Table 8: Results from the NLI4CT Task 1 evidence-only baseline on the test set.  Epoch | F1 score | Accuracy
---|---|---
Pre-trained | 0.008 | 0.500
1 | 0.666 | 0.506
2 | 0.527 | 0.590
3 | 0.623 | 0.572
4 | 0.613 | 0.604
5 | 0.522 | 0.578
6 | 0.485 | 0.584
7 | 0.572 | 0.590
8 | 0.575 | 0.598
9 | 0.494 | 0.590
10 | 0.626 | 0.606
Table 9: BERT-base results on the Statement-only test set. Table 10:
Additional model implementation details for Task 1. Figure 6: Graph of the
BioBERT F1 score by epoch, on the test, training, and development set. Figure
7: Graph of BERT-base F1 score by epoch, on the test set, compared to Bio-
BERT. Figure 8: Distribution curve of the number of tokens per statement in
the test and training set. Figure 9: Page 1 of example annotation form. Figure
10: Page 2 of example annotation form.


# On the suitability of rigorous coupled-wave analysis for fast optical force
simulations

Bo Gao1, Henkjan Gersen1,2 and Simon Hanna1 1 H.H. Wills Physics Laboratory,
University of Bristol, BS8 1TL, Bristol, UK  $2$ iLoF - Intelligent Lab on
Fibre, Porto 4300-240, Portugal<EMAIL_ADDRESS>
###### Abstract

Optical force responses underpin nanophotonic actuator design, which requires
a large number of force simulations to optimize structures. Commonly used
computation methods, such as the finite-difference time-domain (FDTD) method,
are resource intensive and require large amounts of calculation time when
multiple structures need to be compared during optimization. This research
demonstrates that performing optical force calculations on periodic structures
using the rigorous coupled-wave analysis method is typically on the order of
10 times faster than FDTD with sufficient accuracy to suit optical design
purposes. Moreover, this speed increase is available on consumer grade laptops
with a CUDA-compatible GPU avoiding the need for a high performance computing
resource.

††:

Keywords: optical force, rigorous coupled-wave analysis (RCWA), finite-
difference time-domain (FDTD), metamaterials

##  1 Introduction

With the rapid development of nanofabrication techniques and increased
understanding of optical responses, fabricating fast optically driven
actuators becomes feasible. One of the challenges in practice is to design the
structure of a potential actuator in a given material system with the desired
optical response. This process is also known as inverse design[1], and it is
commonly carried out by comparing simulated responses for various structures.

A diverse range of conventional optimization algorithms have been used to
accelerate the search for optimal structures in inverse design, including
genetic algorithms[2], simulated annealing[3], gradient-search algorithms[4]
and topology optimization[5]. These optimization methods are search strategy
driven, because the structures to be simulated depend on decisions made during
the algorithm and cannot be foretold and prepared in advance. Recently, data-
driven algorithms involving machine learning[6] and neural networks[7] have
started to be used[8], where the simulations needed for the training sets are
independent of the optimization routines and can be prepared in advance and
reused.

Although the approaches mentioned above provide optimization strategies that
avoid exploring the entire structural parameter space, they still require a
large number of simulations of various structures. These simulations are
normally the most computationally expensive parts in inverse design routines
due to the high complexity of the structural parameters. As a result, the
suitability of any inverse design method strongly depends on the speed and
required computing resources to obtain an accurate simulation of an individual
structure.

The interest in optical forces for actuation control has exploded in recent
decades due to the development of optical trapping techniques[9]. As a tool
for manipulating the dynamic behaviour of nanophotonic structures[10, 11, 12],
optical forces have the potential for creating fast directly driven
actuators[13]. Although the magnitude of optical forces is usually small, of
the order of pN, and operate on the microscopic scale, research shows that
they can be manipulated and amplified to be observed at a macroscopic
scale[14] typically through the use of repeats of a designed motif[15].
Although many kinds of optical forces are considered in designing actuators
such as radiation pressure and gradient forces[16], we are particularly
interested in optimising so-called lateral optical forces, which typically act
in a direction that is not parallel to the propagation direction of the
incident light. Recent examples of lateral optical forces include the design
of plasmonic linear nanomotors[17], microscopic metavehicles[18], self-
stabilized lightsails[19], and diffractive solar sails[20, 21].

There are several methods available to calculate optical forces directly from
scattering theory, such as the discrete dipole approximation[22] and the
T-matrix method[23]. However, these methods are designed to calculate
scattering from single or simple geometric structure and are generally not
well-suited to deal with repeated complex structures or periodic boundary
conditions.

The finite-difference time-domain (FDTD), finite-difference frequency-domain
(FDFD) and finite element method (FEM) are commonly used for simulating
electromagnetic fields for periodic structures by discretising space and
adding consideration of boundary conditions. As a time-evolution based method,
FDTD can take a long time for a system to reach the steady state, rendering it
inefficient for inverse design studies used for force optimization. The FDFD
and FEM methods are capable of calculating the steady state directly, but
their solution requires the inversion of a huge sparse matrix whose size is
determined by the number of lattice sites required. This can be problematic
when the structure is 3D and requires a fine mesh to model[24].

The rigorous coupled-wave analysis (RCWA) method, as a type of Fourier modal
method which is optimal for monochromatic situations, has been proven to be an
extremely fast and efficient tool for studying reflectance, transmittance and
diffraction efficiencies of multilayered periodic structures[25]. RCWA is not
typically used in the field of nanophotonic inverse design due to its slow
convergence behaviour in dealing with structures much smaller than a
wavelength and high refractive index (RI) contrast materials. Although several
approaches based on RCWA[26, 27, 28] have recently demonstrated that these
difficulties can be overcome, the suitability of performing optical force
calculations using RCWA has, to the best of our knowledge, not been tested.

In this article, we present an optical force simulation approach for periodic
structures suitable for nanophotonic inverse design of optical actuators based
on the RCWA method. To illustrate our approach, we simulate an experimentally
studied structure[18] which acts as a benchmark for comparing force
calculations using the RCWA and FDTD methods. After describing the structure
studied and outlining the principle of optical force calculation using RCWA
and FDTD, we compare optical force calculations of the two methods using two
scenarios that are most commonly found in inverse design: changing of incident
wavelength and structural parameters. We then investigate the convergence and
accuracy for the calculation of optical forces using the two methods. Finally,
we demonstrate that the computational time usage of RCWA is a fraction of that
of FDTD whilst its hardware requirements are those within the availability of
consumer grade laptops.

##  2 Benchmark

In order to compare the computation performance of optical forces calculated
using RCWA and FDTD, a series of structures derived from the meta-vehicle
study of Andrén et al[18] were chosen as a benchmark. These structures consist
of an asymmetric repeating motif, as shown in figure 1, giving rise to a
component of optical force perpendicular to the incident field direction.

Figure 1:  (a) A schematic view of the optical system: linear _x_ -polarized
light propagating along the _z_ direction gets deflected due to the asymmetric
dimer inclusions. The resulting optical force $\bi{F}_{\rm opt}$ is decomposed
into a lateral optical force $\bi{F}_{\rm lof}$ in $x$ direction and a
radiation force $\bi{F}_{\rm rad}$ in $z$ direction. The unit cell structure
(b) and its top view (c) demonstrate the definition of coordinates and key
structural parameters. For details on dimensions and materials see further
details in the main text.

All the structures were illuminated with a normally incident plane-wave
propagating in the $z$ direction and linearly polarized in the $x$ direction
as shown schematically in figure 1(a). Due to the lack of mirror symmetry in
the $x$ direction of the system, the light deflected by the structure is
asymmetric and as a result exerts an optical force $F_{\text{opt}}$ with a
significant lateral component $F_{\text{lof}}$ in the $x$ direction and a
force component related to the radiation pressure $F_{\text{rad}}$ in the $z$
direction. As the structures of interest are periodic arrays of nanoantenna
inclusions, periodic boundary conditions in the $x$ and $y$ directions of the
unit cell are used in the simulation.

Figure 1(b) and figure 1(c) show the unit cell structure used in the
simulation. The unit cell has a periodicity of $\Lambda_{x}=\qty{950}{}$ and
$\Lambda_{y}=\qty{600}{}$ and a total thickness of $H=\qty{1000}{}$, and is
immersed in water (RI 1.33). The dimer nanoantenna inclusion is 200 in width
(w), 460 in height (h), 720 in total length (L) and with a 50 gap (g) that
separates the two blocks. The inclusion is made of polycrystalline Si (RI
3.45), whereas other parts of the unit cell consist of SiO2 (RI 1.45). The
incident wavelength in water $\lambda_{w}$ and the position of the gap centre
$x_{g}$ were modified and the corresponding optical forces exerted on the
structure were studied in the simulations that were performed.

To simplify the simulation, and because the FDTD code, MEEP[29], currently can
only calculate optical forces in a vacuum, we re-scaled all RIs and
wavelengths used in this study by the RI of water, such that the background
material is reduced to vacuum to allow a direct comparison between the two
methods.

##  3 Principle of RCWA and FDTD

Both the FDTD and RCWA methods solve Maxwell’s equations, either in the time-
domain for FDTD:

$\displaystyle\boldsymbol{\nabla}\times\bi{E}(\bi{r},t)$
$\displaystyle=-\partial_{t}\bi{B}(\bi{r},t),$ (1)
$\displaystyle\boldsymbol{\nabla}\times\bi{H}(\bi{r},t)$
$\displaystyle=\partial_{t}\bi{D}(\bi{r},t)+\bi{J}(\bi{r},t),$ (2)

or in the frequency-domain for RCWA:

$\displaystyle\boldsymbol{\nabla}\times\bi{E}(\bi{r},\omega)$
$\displaystyle=-\rmi\omega\bi{B}(\bi{r},\omega),$ (3)
$\displaystyle\boldsymbol{\nabla}\times\bi{H}(\bi{r},\omega)$
$\displaystyle=\rmi\omega\bi{D}(\bi{r},\omega)+\bi{J}(\bi{r},\omega),$ (4)

where $\bi{r}=(x,y,z)$ is the position vector, $\omega$ is the angular
frequency and $\bi{J}$ is the electric current density. The frequency-domain
equations are the Fourier transform of the time-domain equations. In both sets
of equations, the complexity of solving electric and magnetic fields $\bi{E}$
and $\bi{H}$ comes from the electric and magnetic flux densities $\bi{D}$ and
$\bi{B}$, as they rely on the complex relative permittivity
$\tilde{\epsilon}_{\rm r}(\bi{r})$ and permeability $\tilde{\mu}_{\rm
r}(\bi{r})$ distribution of the materials:

$\displaystyle\bi{D}(\bi{r})$ $\displaystyle=\epsilon_{0}\tilde{\epsilon}_{\rm
r}(\bi{r})\bi{E}(\bi{r}),$ (5) $\displaystyle\bi{B}(\bi{r})$
$\displaystyle=\mu_{0}\tilde{\mu}_{\rm r}(\bi{r})\bi{H}(\bi{r}),$ (6)

where $\mu_{0}$ and $\epsilon_{0}$ are the permeability and permittivity of
vacuum, respectively.

FDTD solves (1) and (2) by discretising both space and time and updates the
field values at each time step iteratively[30]. In order to obtain a second-
order accuracy for discretization, both time- and space-derivatives are
approximated with central finite differences, and different field components
are stored on different grid locations known as the Yee Lattice [30]. The Yee
Lattice is designed to allow the $\bi{E}(\bi{r})$ and $\bi{H}(\bi{r})$
components to leap-frog each others in time and space on each iteration.
During this process, the electric and magnetic fields are computed as a
function of time, so that all the secondary properties such as optical forces
and fluxes can be deduced from the primary fields. Since this is a time
evolution simulation, it is necessary to include the source current density in
the simulation region, and begin from a point in time where there are no
fields in space generally. Consequently, any steady state measurement has to
wait for the system to settle, which can be a lengthy process.

The RCWA method, on the other hand, is a semi-analytical method that solves
the frequency-domain equations (3) and (4). The current density term can be
dropped if the region of interest has no current source in it. By considering
a monochromatic incident plane-wave with angular frequency $\omega_{0}$, or
wavenumber $k_{0}$, and substituting (5) and (6) as well as a reduced magnetic
field $\bi{H}_{\rm r}(\bi{r})=-\rmi\sqrt{\mu_{0}/\epsilon_{0}}\bi{H}(\bi{r})$,
(3) and (4) are reduced to

$\displaystyle\boldsymbol{\nabla}\times\bi{E}(\bi{r})$
$\displaystyle=k_{0}\tilde{\mu}_{\rm r}(\bf{r})\bi{H}_{\rm r}(\bi{r}),$ (7)
$\displaystyle\boldsymbol{\nabla}\times\bi{H}_{\rm r}(\bi{r})$
$\displaystyle=k_{0}\tilde{\epsilon}_{\rm r}(\bi{r})\bi{E}(\bi{r}).$ (8)

Here, the reduced magnetic field $\bi{H}_{\rm r}(\bi{r})$ is introduced for
improved numerical robustness as it shares the same order of magnitude with
the electric field $\bi{E}(\bi{r})$.

As a Fourier modal method, RCWA also takes the Fourier transform of space in
the $x$ and $y$ directions, which are perpendicular to the beam propagation
direction. Hence the fields as well as the permittivity and permeability
distribution can be expressed as a Fourier series:

$\displaystyle\bi{E}(x,y,z)$
$\displaystyle=\sum_{m,n}{\bi{e}^{mn}(z)\exp[-\rmi(k_{x}^{mn}x+k_{y}^{mn}y)]},$
(9) $\displaystyle\bi{H}_{\rm r}(x,y,z)$
$\displaystyle=\sum_{m,n}{\bi{h}^{mn}(z)\exp[-\rmi(k_{x}^{mn}x+k_{y}^{mn}y)]},$
(10)

where $\bi{e}^{mn}(z)$ and $\bi{h}^{mn}(z)$ are $z$-direction analytical
amplitude functions of the corresponding Fourier mode. Taking the grating
diffraction equation with periodicity $\Lambda_{x}$ and $\Lambda_{y}$ into
account, the Fourier mode $k_{x}^{mn}$ and $k_{y}^{mn}$ are equivalent to the
$x$ and $y$ components of the wavevector of the $mn^{\rm th}$ order diffracted
plane wave:

$\displaystyle k_{x}^{mn}$ $\displaystyle=k_{{\rm in},x}-2\pi m/\Lambda_{x},$
(11) $\displaystyle k_{y}^{mn}$ $\displaystyle=k_{{\rm in},y}-2\pi
n/\Lambda_{y},$ (12)

where $k_{{\rm in},x}$ and $k_{{\rm in},y}$ are the $x$ and $y$ components of
the wavevector of the incident plane wave. Similarly, the relative
permittivity and permeability can be expressed in reciprocal space as
${\tilde{\epsilon}_{\rm r}}^{mn}(z)$ and ${\tilde{\mu}_{\rm r}}^{mn}(z)$ such
that their distribution only depends on the $z$ coordinate:

$\displaystyle{\tilde{\epsilon}_{\rm r}}(x,y,z)$
$\displaystyle=\sum_{m,n}{\epsilon_{\rm
r}^{mn}(z)\exp[-\rmi(k_{x}^{mn}x+k_{y}^{mn}y)]},$ (13)
$\displaystyle{\tilde{\mu}_{\rm r}}(x,y,z)$ $\displaystyle=\sum_{m,n}{\mu_{\rm
r}^{mn}(z)\exp[-\rmi(k_{x}^{mn}x+k_{y}^{mn}y)]}.$ (14)

Considering a structure can be divided into several layers which are uniform
in the $z$ direction, the distribution of $\tilde{\epsilon}_{\rm r}$ and
$\tilde{\mu}_{\rm r}$ in real space depends only on the $x$ and $y$
coordinates for each layer. Then, within each layer, after taking the Fourier
transform that eliminates the dependence on $x$ and $y$ coordinates, the
relative permittivity and permeability in each Fourier mode become constant.

Consequently, the partial derivatives of the $x$ and $y$ components are
reduced into algebraic products of a factor of $\rmi k_{x}^{mn}$ and $\rmi
k_{y}^{mn}$ for the $mn^{\rm th}$ plane wave mode, respectively. Meanwhile,
after taking the Fourier transform, the product of the permittivity and the
electric field (or the permeability and the magnetic field) becomes a
convolution of their Fourier transforms. In numeric computation, the total
number of plane wave modes is truncated by $N_{\rm modes}=N_{m}\times N_{n}$,
where $N_{m}$ and $N_{n}$ are the number of modes in the $x$ and $y$
directions, respectively. In doing so the convolution can be replaced by a
matrix multiplication by converting the set of constants
${\tilde{\epsilon}_{\rm r}}^{mn}$ and ${\tilde{\mu}_{\rm r}}^{mn}$ into
Toeplitz convolution matrices[31] $\lshad\epsilon_{\rm r}\rshad$ and
$\lshad\mu_{\rm r}\rshad$; the size of the Toeplitz convolution matrix is
$N_{\rm modes}\times N_{\rm modes}$.

Similarly, the $x$ and $y$ components of the wavevector for each plane wave
mode $k_{x}^{mn}$ and $k_{y}^{mn}$ can be expressed in the matrix forms
$\bf{K_{x}}$ and $\bf{K_{y}}$. The $x$ and $y$ components of the amplitude
functions for each Fourier mode defined in (9) and (10) can therefore be
written into column vectors $\bi{e}_{x}(z)$, $\bi{e}_{y}(z)$ and
$\bi{h}_{x}(z)$, $\bi{h}_{y}(z)$.

Eventually, (7) and (8) reduce to:

$\displaystyle\frac{\rmd}{\rmd z}\left[\begin{array}[]{c}\bi{e}_{x}(z)\\\
\bi{e}_{y}(z)\end{array}\right]$ $\displaystyle=\bf{M}$
$\displaystyle\left[\begin{array}[]{c}\bi{h}_{x}(z)\\\
\bi{h}_{y}(z)\end{array}\right],$ (19) $\displaystyle\frac{\rmd}{\rmd
z}\left[\begin{array}[]{c}\bi{h}_{x}(z)\\\ \bi{h}_{y}(z)\end{array}\right]$
$\displaystyle={\bf N}$ $\displaystyle\left[\begin{array}[]{c}\bi{e}_{x}(z)\\\
\bi{e}_{y}(z)\end{array}\right],$ (24)

where $\bf M$ and $\bf N$ are matrices with size $N_{\rm modes}\times N_{\rm
modes}$ and defined by:

$\displaystyle{\bf M}$
$\displaystyle=k_{0}\left[\begin{array}[]{@{}c@{}c@{}}{\bf
K_{x}}\lshad\epsilon_{\rm r}\rshad^{-1}{\bf K_{y}}&\lshad\mu_{\rm
r}\rshad-{\bf K_{x}}\lshad\epsilon_{\rm r}\rshad^{-1}{\bf K_{x}}\\\ {\bf
K_{y}}\lshad\epsilon_{\rm r}\rshad^{-1}{\bf K_{y}}-\lshad\mu_{\rm
r}\rshad&-{\bf K_{y}}\lshad\epsilon_{\rm r}\rshad^{-1}{\bf
K_{x}}\end{array}\right]~{},$ (27) $\displaystyle{\bf N}$
$\displaystyle=k_{0}\left[\begin{array}[]{@{}c@{}c@{}}{\bf
K_{x}}\lshad\mu_{\rm r}\rshad^{-1}{\bf K_{y}}&\lshad\epsilon_{\rm
r}\rshad-{\bf K_{x}}\lshad\mu_{\rm r}\rshad^{-1}{\bf K_{x}}\\\ {\bf
K_{y}}\lshad\mu_{\rm r}\rshad^{-1}{\bf K_{y}}-\lshad\epsilon_{\rm
r}\rshad&-{\bf K_{y}}\lshad\mu_{\rm r}\rshad^{-1}{\bf
K_{x}}\end{array}\right]~{}.$ (30)

By taking another derivative of (19) and (24) with respect to $z$, the
association between the electric and magnetic terms is decoupled:

$\displaystyle\frac{\rmd^{2}}{\rmd
z^{2}}\left[\begin{array}[]{c}\bi{e}_{x}(z)\\\
\bi{e}_{y}(z)\end{array}\right]-{\bf M}{\bf
N}\left[\begin{array}[]{c}\bi{e}_{x}(z)\\\ \bi{e}_{y}(z)\end{array}\right]=0,$
(35) $\displaystyle\frac{\rmd^{2}}{\rmd
z^{2}}\left[\begin{array}[]{c}\bi{h}_{x}(z)\\\
\bi{h_{y}}(z)\end{array}\right]-{\bf N}{\bf
M}\left[\begin{array}[]{c}\bi{h}_{x}(z)\\\ \bi{h}_{y}(z)\end{array}\right]=0,$
(40)

Therefore, within each layer, the two 3D partial differential equations in
real space (7) and (8) are rewritten into two sets of 1D second order ordinary
differential equations with respect to the $z$ coordinate.

In practice, only one set of the equations in (35) and (40) need to be solved
as the other one can be calculated directly by matrix multiplication; here, we
solve for (35). Since it is a wave equation in the matrix form, its solution
can be found analytically by using eigendecomposition[32]:

$\left[\begin{array}[]{c}\bi{e}_{x}(z)\\\ \bi{e}_{y}(z)\end{array}\right]={\bf
W}e^{-\boldsymbol{\kappa}z}\bi{c}^{+}+{\bf
W}e^{+\boldsymbol{\kappa}z}\bi{c}^{-},$ (41)

where $\boldsymbol{\kappa}$ and ${\bf W}$ are the diagonal eigenvalue matrix
and the corresponding eigenvector matrix of the matrix ${\bf MN}$.
$\bi{c}^{+}$ and $\bi{c}^{-}$ represent light that enters into the given layer
with forward and backward propagation directions, respectively. The magnetic
terms can be calculated from electric terms by using (24) giving:

$\displaystyle\left[\begin{array}[]{c}\bi{h}_{x}(z)\\\
\bi{h}_{y}(z)\end{array}\right]$ $\displaystyle=-{\bf
NW}\boldsymbol{\kappa}\rme^{-\boldsymbol{\kappa}z}\bi{c}^{+}+{\bf
NW}\boldsymbol{\kappa}\rme^{+\boldsymbol{\kappa}z}\bi{c}^{-}$ (44)
$\displaystyle=-{\bf V}\rme^{-\boldsymbol{\kappa}z}\bi{c}^{+}+{\bf
V}\rme^{+\boldsymbol{\kappa}z}\bi{c}^{-},$ (45)

where ${\bf V}={\bf NM}\boldsymbol{\kappa}$. Therefore, the overall solution
${\boldsymbol{\Psi}}_{i}(z)$ within the $i^{\rm th}$ layer with thickness
$L_{i}$ can be written as:

$\displaystyle{\boldsymbol{\Psi}}_{i}(z)$
$\displaystyle=\left[\begin{array}[]{@{}c@{}}\bi{e}_{x,i}(z)\\\
\bi{e}_{y,i}(z)\\\ \bi{h}_{x,i}(z)\\\ \bi{h}_{y,i}(z)\end{array}\right]$ (50)
$\displaystyle=\left[\begin{array}[]{@{}cc@{}}{\bf W}_{i}&{\bf W}_{i}\\\
{\bf-V}_{i}&{\bf
V}_{i}\end{array}\right]\left[\begin{array}[]{@{}c@{}c@{}}\rme^{-\boldsymbol{\kappa}_{i}z}&0\\\
0&\rme^{\boldsymbol{\kappa}_{i}z}\end{array}\right]\left[\begin{array}[]{@{}c@{}}\bi{c}_{i}^{+}\\\
\bi{c}_{i}^{-}\end{array}\right].$ (57)

The scattering matrix (S-matrix) can be constructed by considering the
boundary between the $i^{\rm th}$ and ${(i+1)}^{\rm th}$ layers:

$\boldsymbol{\Psi}_{i}(z=z_{i}+L_{i})=\boldsymbol{\Psi}_{i+1}(z=z_{i+1}),$
(58)

where $z_{i}$ and $z_{i+1}$ are the lower boundary coordinates of the
corresponding layer. The multiple reflection between each layer can be handled
by using the Redheffer star product to stack the S-matrices of each layer to
form a global S-matrix that represents the entire structure[32]. Since the
solution is expressed in Fourier modes and satisfies the plane-wave
decomposition, the diffraction efficiency of each plane wave mode can be found
directly by applying the corresponding (reflection and transmission) S-matrix
block on the incident light.

##  4 Optical force calculation

Optical forces $\bi{F}_{\rm opt}$ typically originate from absorption and
scattering of light momentum[33]. Direct calculation of optical forces
requires knowledge of the electromagnetic nearfield distribution in the
vicinity of the structure, which is usually performed by a surface integral of
the Maxwell stress tensor $\boldsymbol{\sigma}$ over a surface $S$ enclosing
the object or region of interest:

$\bi{F}_{\rm opt}=\oint_{S}\boldsymbol{\sigma}\boldsymbol{\cdot}\rmd\bi{A}.$
(59)

The components of the Maxwell stress tensor $\sigma_{\alpha\beta}$ in a vacuum
can be found easily by using the local electric and magnetic field components:

$\sigma_{\alpha\beta}=\epsilon_{0}E_{\alpha}E_{\beta}+\mu_{0}H_{\alpha}H_{\beta}-{\textstyle\frac{1}{2}}\delta_{\alpha\beta}(\epsilon_{0}|\bi{E}|^{2}+\mu_{0}|\bi{H}|^{2})~{},$
(60)

where $\alpha,\beta\in\\{x,y,z\\}$. This is the force calculation method used
in the FDTD simulation as the nearfield information is easy to obtain for the
discretization.

Optical forces $\bi{F}_{\rm opt}$ can also be calculated indirectly by
considering the conservation of energy and momentum. The average optical force
exerted on a scattering object in one time period of oscillation $T$ can be
deduced from the momentum difference of the scattered and incident fields:

$\langle\bi{F}_{\rm opt}\rangle=-\frac{\bi{p}_{\rm scat}-\bi{p}_{\rm
in}}{T}~{}.$ (61)

where $\bi{p}_{\rm scat}$ and $\bi{p}_{\rm in}$ are the momentum of scattered
and incident fields, respectively. Since the optical force is linear with the
incident power for a given structure, a dimensionless coefficient
$\boldsymbol{\mathcal{F}}_{\rm opt}$ is used to characterize the force
response resulting from scattering:

$\boldsymbol{\mathcal{F}}_{\rm opt}=\frac{\langle\bi{F}_{\rm
opt}\rangle}{P_{\rm in}/c},$ (62)

where $P_{\rm in}$ is the incident power.

Since the total field calculated by the RCWA method is a summation over all
the modes in the plane wave decomposition, the total momentum carried by the
scattered field is contributed to linearly by the momentum of each mode. The
linear momentum $\bi{p}^{mn}$ of the $mn^{\rm th}$ plane wave mode can be
found directly from its unit wave vector $\hat{\bi{k}}^{mn}$ and radiation
energy $\mathcal{U}^{mn}$:

$\bi{p}^{mn}=\frac{\mathcal{U}^{mn}n_{0}}{c}\hat{\bi{k}}^{mn}~{},$ (63)

where $n_{0}$ is the RI of the surrounding medium. Therefore, the
dimensionless force coefficient $\boldsymbol{\mathcal{F}}_{\rm opt}$ can be
found by:

$\displaystyle\boldsymbol{\mathcal{F}}_{\rm opt}$
$\displaystyle=\frac{(\mathcal{U}_{{\rm in}}\hat{\bi{k}}_{\rm
in}-\sum_{m,n}{\mathcal{U}^{mn}\hat{\bi{k}}^{mn})n_{0}/cT}}{P_{\rm in}/c}$
(64) $\displaystyle=n_{0}(\hat{\bi{k}}_{\rm
in}-\sum_{m,n}{a^{mn}\hat{\bi{k}}^{mn}}),$

where $a^{mn}=P^{mn}/P_{\rm in}$ is the corresponding diffraction efficiency
for each mode directly calculated from the RCWA method, $\mathcal{U}_{\rm in}$
is the radiation energy of incident light, and $P^{mn}$ is the power of the
$mn^{\rm th}$ plane wave mode. Similar to the optical force $\bi{F}_{\rm opt}$
shown in figure 1, the dimensionless force coefficient
$\boldsymbol{\mathcal{F}}_{\rm opt}$ can also be decomposed into a lateral
component $\mathcal{F}_{x}$ that corresponds to the lateral optical force and
a forward component $\mathcal{F}_{z}$ that corresponds to the radiation
pressure.

##  5 Implementation

The FDTD method was performed using MEEP[29], an open-source software package
implemented in C++ with a Python interface. It was parallelized using the
Message Passing Interface (MPI). The RCWA method was implemented in Python
using the PyTorch package for GPU acceleration of matrix operations with CUDA.

Most of simulations were carried out on a single node of BlueCrystal Phase4
(BC4), a high performance computing (HPC) platform at the University of
Bristol. The FDTD simulations on BC4 were performed by two Intel E5-2860 v4
CPU (28 physical cores in total), whereas the RCWA simulations on BC4 were
executed by one of the CPUs and an Nvidia P100 GPU with 16GB graphical memory.

Two consumer grade laptops were also used to compare the computing
performances of RCWA and FDTD with BC4. Laptop A has an 8 core Intel i7-11800H
CPU and an Nvidia GeForce RTX3080 GPU with 8GB graphical memory, while laptop
B has a 12 core Intel i7-1270P CPU and an Nvidia GeForce MX550 GPU with 2GB
graphical memory.

##  6 Results and Discussion

###  6.1 Force spectrum

In order to verify that our implementation of the RCWA method provides a
correct calculation of the force coefficients, the force spectrum of the meta-
vehicle structure[18] was calculated by varying the incident wavelength using
both the RCWA and FDTD approaches. Simulations were performed using the
structure shown in figure 1 with the parameters taken from the paper[18] and
listed in section 2. Figure 2(a) and figure 2(b) show the lateral and forward
optical force coefficients $\mathcal{F}_{x}$ and $\mathcal{F}_{z}$ versus
incident wavelength in water, computed for wavelengths ranging from 300nm to
1200nm. Results are superimposed from the FDTD and RCWA approaches showing
excellent agreement between the two methods, validating the RCWA
implementation.

Figure 2:  Dimensionless force coefficients (a) $\mathcal{F}_{x}$ and (b)
$\mathcal{F}_{z}$ calculated by FDTD and RCWA for the structure shown in
figure 1 with $x_{g}=\qty{65}{nm}$ and different incident wavelengths. The
dark grey background represents the region where only zero-order diffraction
exists, whereas the lighter one represents the region where second-order
diffraction exists. The red and blue dashed lines mark the two wavelengths in
water $\lambda_{w}=\qty{800}{}$ and $\lambda_{w}=\qty{900}{}$, respectively,
to be investigated in the convergence study .

A number of features of the spectra deserve comment. Due to the periodicity of
the structure, and taking the grating diffraction equation into account, in
the region where the wavelength in water $\lambda_{w}>\qty{950}{nm}$, the
lateral optical force coefficient $\mathcal{F}_{x}$ is constantly zero as only
the zero order diffraction is permitted and there is no lateral momentum
transfer in this region. This region is highlighted by the dark grey
background in the right hand region of figure 2. It can be seen that the
optical forces vary dramatically with the wavelength in water, needing a high
sampling density to fully demonstrate details in the spectrum. Any local
maxima and minima or narrow peaks and valleys appearing in the spectrum may be
an indication of interference effects within the structure or possible
resonance effects. For example, the two huge peaks in forward optical force
coefficient $\mathcal{F}_{z}$ spectrum at $\lambda_{w}=\qty{970.5}{}$ and
$\lambda_{w}=\qty{1006.5}{}$ are due to total (or almost total) reflection.
The two valleys located at $\lambda_{w}=\qty{984.0}{}$ and
$\lambda_{w}=\qty{1191.0}{}$ where $\mathcal{F}_{z}$ is nearly zero
corresponds to total transmission.

The wavelength used by Andrén et al[18] 1064 ($\lambda_{w}=\qty{800}{}$)
marked by the red dashed line appears to be a local minimum of
$\mathcal{F}_{z}$ which may result from an interference effect around the
inclusion layer. This can be seen in terms of thin film interference by taking
the effective refractive index[34] (ERI) averaged by volume fraction into
account. Since the inclusion layer has the height of $h=\qty{460}{nm}$ and ERI
of 2.86, the effective optical path length of the layer is 1316nm, which is
about 5/4 of the incident wavelength (1064nm). This indicates that the light
reflected by the inclusion layer has a phase difference of $\pi$ with the
incident light. The blue dashed line with $\lambda_{w}=\qty{900}{}$ marks a
place where the optical force coefficient spectrum is rather flat and shows no
sign of strong interference or resonance effects. This point will be examined
in more detail in the convergence studies (see below).

The force spectrum calculated by FDTD can be performed in a single simulation
provided that the light source is broadband and able to cover all wavelengths
of interest. This follows from Fourier transforming the time sequence data
[35]. However, one of the disadvantages of this technique is that when the
wavelength of interest is too far from the central wavelength of the light
source, the fraction of power in that wavelength is small and the error in the
force coefficients can be amplified due to machine truncation error.
Furthermore, as the time sequence of fields is truncated when the system
reaches the criterion for steady state, taking a Fourier transform also
produces a termination effect. This is the reason the FDTD curves in figure 2
show larger noise, especially in the zero-order diffraction region (dark grey
region). Although these fluctuations can be reduced by running the simulation
for a longer time or applying a window to the time sequence to reduce
termination error, the most reliable solution is to run the simulation
multiple times for each wavelength individually.

###  6.2 Varying gap centre

To further compare the consistency between RCWA and FDTD and their potential
behaviour in inverse design, a set of simulations that varies the position of
the gap centre $x_{g}$ (see figure 1(c)) was performed, while the incident
wavelength in water was kept unchanged at $\lambda_{w}=\qty{800}{nm}$. These
simulations, as an example of geometry optimization, were designed to find the
optimal gap position for the largest lateral optical force coefficient
$\mathcal{F}_{x}$. Figure 3 shows the calculated optical force coefficients
versus gap centre position ranging from $x_{g}=\qty{-335}{nm}$ to
$x_{g}=\qty{335}{nm}$. The total length of the inclusion and the width of the
gap in between was kept unchanged as $L=\qty{720}{nm}$ and $g=\qty{50}{nm}$.

Figure 3:  Force coefficients $\mathcal{F}_{x}$ and $\mathcal{F}_{z}$
calculated by FDTD and RCWA methods with different gap centre position $x_{g}$
while the total length of the inclusions and the gap width in between are kept
unchanged as 720nm and 50nm. The red and green dashed lines indicate gap
positions $x_{g}=\qty{65}{nm}$ and $x_{g}=\qty{100.5}{nm}$, which correspond
to the original structure[18] and a more optimal structure for lateral optical
force coefficient $\mathcal{F}_{x}$ response, respectively.

The calculated optical force coefficients from RCWA and FDTD show no
significant differences. The original structure ($x_{g}=\qty{65}{nm}$) from
Andrén et al[18] is indicated by the red dashed line with
$\mathcal{F}_{x}=0.44$. However, the optimal gap position was found at
$x_{g}=\qty{100.5}{nm}$ as indicated by the green dashed line with
$\mathcal{F}_{x}=0.55$.

In addition, an odd symmetry can be observed for $\mathcal{F}_{x}$ and an even
symmetry for $\mathcal{F}_{z}$ about the gap position $x_{g}=\qty{0}{nm}$,
which is to be expected due to the fact that, at this gap position, the two
blocks on either side of the gap have the same length of 335nm and the mirror
symmetry in the $x$ direction is regained. Consequently, the lateral optical
force coefficients $\mathcal{F}_{x}$ is 0 at this point, and for similar
reasons at the left and right end as well.

###  6.3 Convergence

Since the force spectrum study and the varying gap centre study confirm that
the optical force coefficients calculated from the RCWA and FDTD methods are
in good agreement, it is important to find the condition for each method where
they can work most efficiently. In each case this will be a compromise between
the need for the calculated optical force coefficient to be sufficiently
accurate whilst the computing resources are kept as small as possible.

The computing resources required and the accuracy of the simulation are
controlled directly by the number of lattice sites per (resolution) for FDTD,
and the total number of plane wave modes in each direction (number of modes)
for RCWA. Higher resolution or a higher number of modes indicate more precise
modelling of the structure, but require greater computing resources in terms
of memory and time usage.

Figure 4:  Convergence of $x$ and $z$ components of optical force,
coefficients $\mathcal{F}_{x}$ and $\mathcal{F}_{z}$, using (a) the FDTD
method and (b) the RCWA method for three models. The key parameters of the
three models are incident wavelength in water $\lambda_{w}$ and the $x$
coordinate of the gap centre position $x_{g}$, which can be found in table 1.
The black dashed lines represent the condition taken as a standard criterion
for convergence and the calculated optical force coefficients at these points
are compared in table 2.

Figure 4 shows convergence of the lateral optical force coefficient
$\mathcal{F}_{x}$ and the forward optical force coefficient $\mathcal{F}_{z}$
of three models as a function of the resolution for FDTD and the number of
modes for RCWA. The key parameters of the three models are the incident
wavelength in water $\lambda_{w}$ and the $x$ coordinate of the gap centre
position $x_{g}$, which are summarized in table 1. The “flat $\lambda_{w}$”
model, as marked by the blue dashed line in figure 2 with
$\lambda_{w}=\qty{900}{nm}$ and $x_{g}=\qty{65}{nm}$, represents a region of
the force spectrum where possible resonance and interference effects within
the structure are less significant and the force response versus wavelength is
relatively flat. The “optimal $x_{g}$” model, as indicated by the green dashed
line in figure 2 with $\lambda_{w}=\qty{800}{nm}$ and $x_{g}=\qty{100.5}{nm}$,
corresponds to the optimal gap position with largest lateral optical force
coefficient $\mathcal{F}_{x}$ as discussed previously in section 6.2. The
“original” model, as highlighted by the red dashed line in both figures with
$\lambda_{w}=\qty{800}{nm}$ and $x_{g}=\qty{65}{nm}$, refers to the parameters
used by Andrén et al[18].

Table 1:  Parameters used in the three convergence studies. The total length
of inclusion $L=\qty{720}{nm}$ and the gap width $g=\qty{50}{nm}$ are kept
unchanged.

  * Model | Wavelength $\lambda_{w}$ (nm) | Gap centre $x_{g}$ (nm)
---|---|---
original (red) | 800 |  065.0
flat $\lambda_{w}$ (blue) | 900 |  065.0
optimal $x_{g}$ (green) | 800 | 100.5

All the optical force coefficients are varying with a fractional difference
between neighbouring data points less than 1% at and after the 100^-1
resolution for FDTD or 41 modes on each direction ($41^{2}$ in total) for
RCWA, as marked by the black dashed lines in figure 4. Therefore the values at
100^-1 for FDTD and 41 modes for RCWA were taken to represent a standard
criterion for convergence across all our simulations, representing a
convenient compromise between possible accuracy and available computing
resources. The number of modes includes both positive and negative modes and
the zero order mode, and hence will always be an odd integer. Table 2 compares
the optical force coefficients taken at these standard points and their
relative differences ($|\text{RCWA}-\text{FDTD}|/|\text{FDTD}|\times 100\%$)
as well. It shows that the calculated optical force coefficients from the two
methods under standard conditions have good agreement and their differences
are all no larger than 2%.

Table 2:  Converged optical force coefficients computed using the FDTD and
RCWA methods, taken at the standard points 100^-1 or 41 modes as shown in
figure 4. The differences between the two methods in the three cases are also
shown and are all within 2%.

  * Model | FDTD | RCWA | Differences
---|---|---|---
$\mathcal{F}_{x}$: original | 0.4378 | 0.4396 | 0.41%
$\mathcal{F}_{z}$: original | 0.7536 | 0.7454 | 1.09%
$\mathcal{F}_{x}$: flat $\lambda_{w}$ | 0.2679 | 0.2643 | 1.33%
$\mathcal{F}_{z}$: flat $\lambda_{w}$ | 0.6169 | 0.6069 | 1.62%
$\mathcal{F}_{x}$: optimal $x_{g}$ | 0.5596 | 0.5519 | 1.39%
$\mathcal{F}_{z}$: optimal $x_{g}$ | 0.6761 | 0.6852 | 1.36%

Close examination shows that the fluctuation of the RCWA results for the
$\mathcal{F}_{z}$ curve of the “original” model after the chosen convergence
point is slightly larger than for the FDTD method. This variation cannot be
removed by simply adding more modes and it appears to be due to the Gibbs
phenomenon which occurs when taking Fourier transforms of discontinuous
systems. In the present case we are dealing with a discontinuous permittivity
and permeability at the boundary of different materials in our test structure.
This is the main difficulty for RCWA when dealing with structures with
characteristic sizes which are much smaller than the incident wavelength.
Resonance or interference effects in the system at certain wavelengths might
also amplify such fluctuations, which could explain the $\mathcal{F}_{z}$
curve of the “original” model in the RCWA method in figure 4.

Convergence of the near field in RCWA has previously been discussed and shown
to be much more problematic than for the far field[28]. This is because the
near field region relies on accurate calculation of both propagating
diffraction orders and evanescent modes, which are particularly susceptible to
the Gibbs phenomenon[28]. Thus the suitability of direct optical force
calculation using the near field results from RCWA using, for example, the
Maxwell stress tensor, is questionable. However, the influence of the Gibbs
phenomenon on far field properties such as diffraction efficiencies can be
controlled within a reasonable degree as seen in figure 4, and the consequent
fluctuations in the optical forces calculated indirectly are within 1%.
Therefore, we conclude that RCWA is suitable for indirect optical force
calculation using momentum conservation, as in (61) and (62).

Concerning the resolution limit for FDTD, any lower resolution than 100^-1 is
not a reasonable choice for the converged point because the 100^-1 resolution
only gives 5 lattice sites to represent the 50nm gap. Thus any lower
resolution cannot describe the geometry of the structure adequately.

###  6.4 Speed and computing resources

In order to test the suitability of the RCWA method for optical force
optimization, aside from the ability of yielding accurate results, the speed
and computing resource requirement must also be considered. During the
simulations mentioned above, we observed that the RCWA method is typically on
the order of 10 times faster than the FDTD method on BC4 for a given force
calculation, whilst its computing resource requirements are generally within
the availability of a consumer grade laptop.

Figure 5 shows the time usage on BC4 for the “original” model in the
convergence study with respect to memory usage. It shows that the time usage
of RCWA method for a single force evaluation is a fraction of that of the FDTD
method. Simulations of other cases perform similarly, _i.e._ , the time usage
at the standard points, 100^-1 resolution for FDTD and 41 modes for RCWA, take
roughly 400 seconds and 60 seconds on BC4, respectively.

Figure 5:  Computational time usage (seconds) versus memory usage (MB) for the
“original” model with parameters shown in table 1 using the FDTD method with
different resolution and the RCWA method with different number of modes
performed on BC4. The solid points represent the time and memory usage taken
at the standard points 100^-1 and 41 modes as shown in figure 4. Arrows
indicate increasing resolution or number of modes.

In the force spectrum study, the time usage of the FDTD simulation with a
broadband source is about 35 minutes on BC4, compared with roughly 201 minutes
(3.5 hours) for the RCWA method. RCWA is much slower than the FDTD in this
situation because it cannot take advantage of Fourier transforming the time
sequence data, as explained previously, and therefore it has to simulate for
each wavelength of interest individually. The 201 minutes time usage is
expected because there were 201 different wavelengths involved in figure 2,
and each RCWA simulation takes about 1 minute. On the other hand, the FDTD
time usage in figure 2 is greater than its typical value, _i.e._ 400 seconds,
because the broadband light source needs longer for the simulation to reach a
steady state. However, a disadvantage of Fourier transforming the time series
data from the FDTD simulation is that it introduces termination effects in the
force spectrum, and the simplest way to eliminate them is to run the FDTD
simulation for each wavelength individually, which would take about 28 hours
on BC4.

In the gap centre variation study (figure 3), for both RCWA and FDTD, each
geometry was investigated with a single separate simulation. As there were 21
different geometries in the study, the time usage for the FDTD method was
about 3.25 hours whilst the RCWA method only took about 21 minutes. Therefore
the average time usage for a single simulation in this case is about 550
seconds for FDTD and 60 seconds for RCWA. The slowing of the FDTD method
compared with the standard case, appears to be because some geometries take
longer to reach the criterion of steady state than others. Overall, these
timings indicate that for a typical inverse design problem requiring force
responses for a large number of different structures, the RCWA method is
around an order of magnitude more efficient than the FDTD method.

Given the apparent speed benefit of RCWA when running on a HPC platform, we
decided to test the feasibility of running the same code on more accessible,
consumer grade, equipment. Table 3 shows a time usage comparison of two
laptops with BC4 for our “original” model at the standard point (100^-1
resolution for FDTD and $41\times 41$ modes for RCWA). Single precision
(64bits) and double precision (128bits) data types were tested because the
graphics cards on BC4 are optimized for double precision whereas the graphics
cards of the most commonly available consumer grade laptops are configured for
predominantly single precision arithmetic.

Table 3:  Time usage comparison of FDTD and RCWA method on the original model
($\lambda_{w}=\qty{800}{}$, $x_{g}=\qty{65}{nm}$) with 100^-1 resolution
(FDTD) and $41\times 41$ modes (RCWA) on different machines.

  * Brand | BlueCrystal Phase4 (HPC) | Laptop A | Laptop B
---|---|---|---
CPU: cores | Intel E5-2680 v4: 28 cores | Intel i7-11800H: 8 cores | Intel i7-1270P: 12 cores
GPU: memory | Nvidia Tesla P100: 16 GB | Nvidia GeForce RTX3080: 8 GB | Nvidia GeForce MX550: 2GB
Time FDTD (seconds) |  0401.0 | 1701.8 | 1563.0
Time RCWA-single (seconds) |  0050.4 |  0049.6 |  0054.6
Time RCWA-double (seconds) |  0058.6 |  0092.2 | —

The single precision RCWA results show that the time usage on both laptops is
about the same as BC4 and could be even faster if its GPU had better
performance. Although the double precision RCWA method on Laptop A is 1.5
times slower, it is still about 20 times faster than the FDTD method on the
laptop and more than 4 times faster than the FDTD method on BC4. Therefore, it
is clear that the RCWA method can also work on consumer grade computers and
laptops with a good performance. This strengthens the suitability of RCWA for
optical force optimization, because access to HPC systems tends to be
restricted or expensive, while CUDA capable graphics cards are readily
available.

Single precision variables are not recommended for use in practice, as they
bring higher truncation and rounding error which would potentially be
amplified by the interference effects discussed above. Figure 5 also indicates
that the RCWA method is more memory-intensive than the FDTD method, and its
memory usage is mainly the graphical memory. The missing data in table 3 is
due to the fact that running 41 modes in double precision requires 2070MB
graphical memory which is larger than the graphical memory capacity (2GB) of
Laptop B. The available graphical memory of a GPU limits the largest number of
modes in each direction that can be used on it, and it is important to
estimate the maximum graphical memory usage during the computation. The
graphical memory usage (in bytes) is $768N^{4}$ for double precision and
$384N^{4}$ for single precision for a given number of modes, $N$, on each
direction. This suggests an upper limit for double precision work of $N=39$
for the 2GB GPU in Laptop B and $N=57$ for the 8GB GPU in Laptop A.

Therefore, the RCWA method is more suitable for fast optical force simulation
than FDTD on a consumer grade laptop with a CUDA capable graphics card and a
graphical memory larger than 2GB, which opens up opportunities for optical
force based inverse design and geometry optimization without the need of HPC
resources.

##  7 Conclusion

This paper shows that optical force coefficients can be calculated using the
RCWA method giving results which are consistent with the FDTD method for a
fraction of the computer time. The FDTD method converged when the resolution
reached 100^-1, with a calculation time about 400 seconds on the local HPC
system. On the other hand, the RCWA method converged when the number of modes
in each direction reached 41, and took about 60 seconds on the same HPC
system. The differences in optical force coefficients for the benchmark case
between RCWA and FDTD methods are no larger than 2%.

The RCWA method is capable of running on both HPC and consumer grade laptops
within a reasonable time as long as the laptop has a CUDA-compatible GPU and
graphical memory larger than 2GB (for double precision and 41 modes on each
direction). RCWA is typically about 10 times faster than FDTD on HPC system
such as BC4, and about 20 times faster on consumer grade laptops. The FDTD
method could be made faster by using more HPC nodes, but such facilities are
not generally available. Therefore, the RCWA method is very suitable for
studies that include calculating optical force responses for a large number of
distinct structures. Furthermore, in the consideration of green computing,
using the RCWA method on personal laptops could save considerable resources
and energy compared to using the HPC, and is therefore the more sustainable
option.

BG acknowledges the China Scholarship Council (CSC) and University of Bristol
for a postgraduate scholarship (award number 202008060217), and Advanced
Computing Research Centre (ACRC) at the University of Bristol for High
Performance Computing platform BlueCrystal phase4 (BC4) support. HG
acknowledges support through an Impact Acceleration Account co-funded by EPSRC
and Bristol Nano Dynamics (grant number EP/R511663/1).

## References

## References

  * [1] Sean Molesky, Zin Lin, Alexander Y. Piggott, Weiliang Jin, Jelena Vucković, and Alejandro W. Rodriguez.  Inverse design in nanophotonics.  Nature Photonics, 12:659–670, 2018.
  * [2] Peter R. Wiecha, Arnaud Arbouet, Christian Girard, Aurélie Lecestre, Guilhem Larrieu, and Vincent Paillard.  Evolutionary multi-objective optimization of colour pixels based on dielectric nanoantennas.  Nature Nanotechnology, 12(2):163–169, 2017.
  * [3] Yaxin Xie, Menghao Liu, Tianhua Feng, and Yi Xu.  Compact disordered magnetic resonators designed by simulated annealing algorithm.  Nanophotonics, 9(11):3629–3636, 2020.
  * [4] Momchil Minkov, Ian A. D. Williamson, Lucio C. Andreani, Dario Gerace, Beicheng Lou, Alex Y. Song, Tyler W. Hughes, and Shanhui Fan.  Inverse design of photonic crystals through automatic differentiation.  ACS Photonics, 7(7):1729–1741, 2020.
  * [5] Mingkun Chen, Jiaqi Jiang, and Jonathan A. Fan.  Design space reparameterization enforces hard geometric constraints in inverse-designed nanophotonic devices.  ACS Photonics, 7(11):3141–3151, 2020.
  * [6] Peter R. Wiecha and Otto L. Muskens.  Deep learning meets nanophotonics: A generalized accurate predictor for near fields and far fields of arbitrary 3d nanostructures.  Nano Letters, 20(1):329–338, 2020.
  * [7] Dianjing Liu, Yixuan Tan, Erfan Khoram, and Zongfu Yu.  Training deep neural networks for the inverse design of nanophotonic structures.  ACS Photonics, 5(4):1365–1369, 2018.
  * [8] Sunae So, Trevon Badloe, Jaebum Noh, Jorge Bravo-Abad, and Junsuk Rho.  Deep learning enabled inverse design in nanophotonics.  Nanophotonics, 9(5):1041–1057, 2020.
  * [9] Arthur Ashkin.  Optical trapping and manipulation of neutral particles using lasers.  Proceedings of the National Academy of Sciences, 94(10):4853–4860, 1997.
  * [10] Hang Li, Yongyin Cao, Lei-Ming Zhou, Xiaohao Xu, Tongtong Zhu, Yuzhi Shi, Cheng-Wei Qiu, and Weiqiang Ding.  Optical pulling forces and their applications.  Advances in Optics and Photonics, 12:288, 2020.
  * [11] D. B. Phillips, M. J. Padgett, S. Hanna, Y. L.D. Ho, D. M. Carberry, M. J. Miles, and S. H. Simpson.  Shape-induced force fields in optical trapping.  Nature Photonics, 8(5):400–405, 2014.
  * [12] D. B. Phillips, G. M. Gibson, R. Bowman, M. J. Padgett, S. Hanna, D. M. Carberry, M. J. Miles, and S. H. Simpson.  An optically actuated surface scanning probe.  Opt. Express, 20(28):29679–29693, 2012.
  * [13] Lip Ket Chin, Yuzhi Shi, and Ai-Qun Liu.  Optical Forces in Silicon Nanophotonics and Optomechanical Systems: Science and Applications.  Advanced Devices & Instrumentation, 2020:1–14, 2020.
  * [14] Hernando Magallanes and Etienne Brasselet.  Macroscopic direct observation of optical spin-dependent lateral forces and left-handed torques.  Nature Photonics, 12(8):461–464, 2018.
  * [15] Michael O’Donnell and Simon Hanna.  Optical forces on patterned particles.  In Complex Light and Optical Forces XIII, volume 10935, pages 110–115. Proc. SPIE, 2019.
  * [16] H. Cai, K. J. Xu, A. Q. Liu, Q. Fang, M. B. Yu, G. Q. Lo, and D. L. Kwong.  Nano-opto-mechanical actuator driven by gradient optical force.  Applied Physics Letters, 100(1):013108, 2012.
  * [17] Yoshito Y. Tanaka, Pablo Albella, Mohsen Rahmani, Vincenzo Giannini, Stefan A. Maier, and Tsutomu Shimura.  Plasmonic linear nanomotor using lateral optical forces.  Science Advances, 6(45):1–7, 2020.
  * [18] Daniel Andrén, Denis G. Baranov, Steven Jones, Giovanni Volpe, Ruggero Verre, and Mikael Käll.  Microscopic metavehicles powered and steered by embedded optical metasurfaces.  Nature Nanotechnology, 16:970–974, 2021.
  * [19] Ramon Gao, Michael D. Kelzenberg, Yonghwi Kim, Ognjen Ilic, and Harry A. Atwater.  Optical characterization of silicon nitride metagrating-based lightsails for self-stabilization.  ACS Photonics, 9(6):1965–1972, 2022.
  * [20] Grover A. Swartzlander.  Theory of radiation pressure on a diffractive solar sail.  J. Opt. Soc. Am. B, 39(9):2556–2563, 2022.
  * [21] Anirban S. Swakshar, Seongsin M. Kim, and Grover A. Swartzlander.  Broadband radiation pressure on a small period diffractive film.  Opt. Express, 30(25):45279–45287, 2022.
  * [22] Lin Ling, Fei Zhou, Lu Huang, and Zhi-Yuan Li.  Optical forces on arbitrary shaped particles in optical tweezers.  Journal of Applied Physics, 108(7):073110, 2010.
  * [23] Paolo Polimeno, Rosalba Saija, Cristian Boschi, Maria Iatì, and Onofrio M. Marago.  Optical forces in the t-matrix formalism.  Atti della Accademia Peloritana dei Pericolanti, 19:A2, 2019.
  * [24] Nathan J. Champagne, James G. Berryman, and H.Michael Buettner.  Fdfd: A 3d finite-difference frequency-domain code for electromagnetic induction tomography.  Journal of Computational Physics, 170(2):830–848, 2001.
  * [25] M. G. Moharam, Eric B. Grann, Drew A. Pommet, and T. K. Gaylord.  Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings.  J. Opt. Soc. Am. A, 12(5):1068–1076, 1995.
  * [26] Ilia M. Fradkin, Sergey A. Dyakov, and Nikolay A. Gippius.  Fourier modal method for the description of nanoparticle lattices in the dipole approximation.  Phys. Rev. B, 99:075310, 2019.
  * [27] Sandeep Inampudi, Mohammad Mahdi Salary, Samad Jafar-Zanjani, and Hossein Mosallaei.  Rigorous space-time coupled-wave analysis for patterned surfaces with temporal permittivity modulation.  Opt. Mater. Express, 9(1):162–182, 2019.
  * [28] Martin Weismann, Dominic F.G. Gallagher, and Nicolae C. Panoiu.  Accurate near-field calculation in the rigorous coupled-wave analysis method.  Journal of Optics, 17(12):125612, 2015.
  * [29] Ardavan F. Oskooi, David Roundy, Mihai Ibanescu, Peter Bermel, J.D. Joannopoulos, and Steven G. Johnson.  Meep: A flexible free-software package for electromagnetic simulations by the fdtd method.  Computer Physics Communications, 181(3):687–702, 2010.
  * [30] Allen Taflove.  Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology.  Artech House, 2013.
  * [31] William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery.  Numerical recipes 3rd edition: The art of scientific computing.  Cambridge university press, 2007.
  * [32] Raymond C Rumpf.  Improved formulation of scattering matrices for semi-analytical methods that is consistent with convention.  Progress In Electromagnetics Research B, 35:241–261, 2011.
  * [33] A. Ashkin.  Applications of laser radiation pressure.  Science, 210(4474):1081–1088, 1980.
  * [34] R. Magnusson and D. Shin.  Diffractive optical components.  In Robert A. Meyers, editor, Encyclopedia of Physical Science and Technology (Third Edition), pages 421–440. Academic Press, New York, third edition, 2003.
  * [35] Wei-Hua Guo, Wei-Jun Li, and Yong-Zhen Huang.  Computation of resonant frequencies and quality factors of cavities by fdtd technique and pade approximation.  IEEE Microwave and Wireless Components Letters, 11(5):223–225, 2001\.