# CAPABILITY ITERATION NETWORK FOR ROBOT PATH PLANNING

Buqing Nie1, Yue Gao∗2, Yidong Mei3 and Feng Gao4 1Buqing Nie is with
Department of Automation, Shanghai Jiao Tong University, Shanghai, P.R. China,
<EMAIL_ADDRESS>Gao is with MoE Key Lab of Artificial Intelligence
and Department of Automation, Shanghai Jiao Tong University, Shanghai, P.R.
China<EMAIL_ADDRESS>Mei is with Department of Automation,
Shanghai Jiao Tong University, Shanghai, P.R. China,
<EMAIL_ADDRESS>Gao is with State Key Laboratory of Mechanical
System and Vibration, Shanghai Jiao Tong University, Shanghai, P.R. China,
<EMAIL_ADDRESS>author.

Abstract

Path planning is an important topic in robotics. Recently, value iteration
based deep learning models have achieved good performance such as Value
Iteration Network(VIN). However, previous methods suffer from slow convergence
and low accuracy on large maps, hence restricted in path planning for agents
with complex kinematics such as legged robots. Therefore, we propose a new
value iteration based path planning method called Capability Iteration
Network(CIN). CIN utilizes sparse reward maps and encodes the capability of
the agent with state-action transition probability, rather than a convolution
kernel in previous models. Furthermore, two training methods including end-to-
end training and training capability module alone are proposed, both of which
speed up convergence greatly. Several path planning experiments in various
scenarios, including on 2D, 3D grid world and real robot with different map
sizes are conducted. The results demonstrate that CIN has higher accuracy,
faster convergence and lower sensitivity to random seed compared to previous
VI-based models, hence more applicable for real robot path planning.

Key Words
path planning, value iteration, value iteration network, capability iteration
network

##  I Introduction

Path planning is an important direction for autonomous robot, whose objective
is to compute a path to the goal region that not only avoids collisions with
obstacles but also satisfies kinematic constraints of the robot [1, 2].
Various classic methods including $A^{*}$[3], Potential Fields[4] and Rapidly
exploring Random Tree[5] are proposed. Recently, Deep Reinforcement
Learning(DRL) has been utilized to solve path planning problems without
handcraft parameters tuning according to the robot topological structures[6],
for example on robotic manipulator[7, 8], drilling robot[9], planetary rover
[10], and long-range navigation task[11].

DRL has the benefits of representing complex state space and learning a
network which enables function approximation of the features and future
rewards values[12]. However, most existing DRL based path planning models lack
planning modules which makes them rely on large amounts of training data and
difficult to generalize for unseen environments[13]. This makes DRL models
difficult to be applied in real scenarios, like planning paths for real
robots. A recent work, _Value Iteration Network_(VIN)[13] combines recurrent
convolution neural networks and max-pooling to emulate Value Iteration(VI)[14]
algorithm, which enables VIN to learn and plan paths for unseen mazes.
However, VIN is particularly challenging on training when given large maps. In
addition, how to utilize VIN to solve path planning for robots with complex
kinematics remains a challenging task.

For RL models, sparse reward is an important property required for convergence
and generalization[15, 16, 17]. Current VIN method utilizes convolution layers
to compute reward map $R(s,a)$, which is not sparse and may decrease the
accuracy of VIN because of its defects, especially on large maps. Besides, VIN
utilizes concatenation and a global convolution kernel $P_{s^{\prime}}^{a}$
(probability distribution for the next state only based on the action to take)
instead of the state conditioned transition probability $P_{ss^{\prime}}^{a}$
in the Value Iteration algorithm.

Figure 1: Structure of the _Capability Iteration Network_ (CIN). CIN utilizes
sparse reward maps and state-conditioned transition probabilities
$P_{ss^{\prime}}^{a}$ instead of reward generated by CNN and global transition
probabilities $P_{s^{\prime}}^{a}$ in VIN.

In this paper, a new VI-based model _Capability Iteration Network_ (CIN) is
proposed. In CIN, sparse reward maps instead of reward maps generated by CNN
are utilized, which limit reward maps influencing the accuracy of the model,
especially on large maps. Besides, as illustrated in fig.1, a capability
module is utilized to encode each state $s$ with a state-conditioned
transition probability $P_{ss^{\prime}}^{a}$, which is similar to the
environment dynamics utilized in model based RL. These state-conditioned
transition probabilities in CIN represents the capability of the agent, which
differentiates itself from previous VI based models that utilize
$P_{s^{\prime}}^{a}$ without the knowledge of the current state. CIN is
trained to learn the capability of the agent, which is relevant to the local
region and independent with the global map. Therefore, CIN is able to be
trained on local maps with faster learning speed, In the experiment section,
several experiments demonstrate that CIN can be applied to scenarios including
2D mazes, large-scale 3D terrain maps and path planning for real hexapod
robots.

The main contributions of this paper are as follows:

  * •

A new VI-based model called CIN is proposed, which greatly improves
performance on accuracy and learning speed compared with previous VI-based
algorithms on various scenarios.

  * •

We provide two training methods: end-to-end training and capability module
training alone on small maps. Both of two methods greatly improves learning
speed compared with the previous VI-based model.

  * •

Experiments on 2D, 3D maps and hexapod robot demonstrate that CIN outperforms
previous models on accuracy and learning speed, and can be utilized for path
planning with real robots.

##  II Related Work

###  II-A Path Planning

Path planning is one of the essential tasks in the automation process of a
system that moves in the environment while avoiding obstacles and respecting
various constraints[18]. It has been widely applied to lots of scenarios
including auto driving, autonomous underwater vehicle control and video
games[19, 20, 21, 22]. It has been widely studied by robot researchers, and
plenty of methods have been proposed.

$A^{*}$ algorithm developed by Hart et al. [3] is widely used in path planning
because of the following properties:(1)$A^{*}$ gives optimality when applied
with the visible graph. (2) the path given by $A^{*}$ is unique (3)
computation cost is small with a good heuristic function. Many various
variants of $A^{*}$ are proposed such as $D^{*}$ Lite[23] and any-angle
$A^{*}$ [21]. Artificial potential field(APF) introduced by Oussama Khatib[4]
is another important method which creates an artificial potential field to
attract the agent around the goal and repulse them around the obstacles[24].
APF is easy to be implemented with smooth paths in real-time. However, the
agent may stay in local minimum regions, which is the main drawback of APF
method. Probabilistic path planning algorithms such as Rapidly-exploring
random trees[5] and probabilistic roadmaps[25] are also effective methods,
which randomly select non-collision points in motion space and then connect
them to find the best path. These algorithms don’t require any environment
modeling, which outperform previous algorithms such as $A^{*}$ in terms of
computation cost[22]. Genetic algorithm is another effective path planning
problem with high robustness in various scenearios such as robot manipulators
and unmanned surface vehicle [26, 27].

Recently, with the development of reinforcement learning and deep learning,
deep reinforcement learning(DRL) has been applied to solve path planning
problems in many scenarios such as robotic manipulator[7, 8], drilling
robot[9] and planetary rover [10]. DRL based methods try to summarize patterns
through numerous attempts to generate a suitable path in a new environment,
which doesn’t require modeling of environment and robot before planning paths.
However, current DRL methods rely heavily on training data, which limits the
application scenarios.

###  II-B Reinforcement Learning

Reinforcement learning (RL) is a powerful paradigm to solve the sequential
decision making problem for Markov Decision Process (MDP), whose objective is
to find an optimal policy for the agent through interaction with the
environment and maximizing the rewards. RL methods can be categorized into
model-free and model-based methods depend on whether the policy has access to
the underlying model of the environment.

Model-free methods learn the policy directly from interactions with the
environment[12]. For example, Mnih et al. proposed Deep Q network [6] which
approximates Q function with a neural network to find the best policy;
Lillicrap et al. proposed Deep Deterministic Policy Gradient[28] which
utilizes policy networks to generate actions while using Q-networks to
criticize the policy. Model-free methods have achieved great success in many
application scenarios such as video games[29] and robot control[30]. However,
they are suffering from low sample efficiency for training data and random
seed sensitivity, which limits the real application of model free methods.

Model-based RL methods can simulate the transition probability of the
environment, resulting in increased sample efficiency, which is particularly
important in domains where each interaction with the environment is
expensive[12]. For example, PILCO[31] and its variants utilize Gaussian
process to learn the environment dynamics and achieve good performance in
nonlinear control problems[32]. World Model[33] utilizes MDN-RNN to simulate
the environment model with raw images as input and outperforms previous
methods in gym tasks. However, learning the environment model may introduces
extra complexities and the errors in the environment model may in turn
influence the policy accuracy. The state-conditioned transition probabilities
utilized in CIN is similar to the environment dynamics in model based RL. In
this paper, we combine the idea of the environment dynamics with the Value
Iteration based deep learning methods to improve the performance in the path
planning tasks.

###  II-C Value Iteration Based Methods

For a markov decision process (MDP) problem with known transition probability
$P_{ss^{\prime}}^{a}$, a dynamic programming method called Value Iteration[14]
can give an optimal policy $\pi^{*}$. The following equations (1) and (2) are
the basis of VI algorithm called the Bellman equations:

$Q(s,a)=R(s,a)+\gamma\sum_{s^{\prime}}P^{a}_{s,s^{\prime}}V_{(}s^{\prime})\\\
$ (1) $V^{\prime}(s)=\max_{a}Q(s,a)$ (2)

The optimal value map $V^{*}(s)$ and Q value map $Q^{*}(s,a)$ can be obtained
by calculating between $V(s)$ and $Q(s)$ iteratively until convergence, after
which the optimal policy $\pi^{*}$ can be obtained utilizing equation (3).

$\pi^{*}=\arg\max_{a}Q^{*}(s,a)$ (3)

Figure 2: Structure of the _Value Iteration Network_ (VIN)[13] VIN combines
the VI algorithm with CNN to implement a RL model with planning capability.

As is shown in Fig. 2, in order to implement planning computation in existing
RL models, Tamar et al.[13] proposed a VI-based RL model called _Value
Iteration Network_ (VIN). VIN combines value iteration with convolution and
max-pooling operators to plan paths with unknown environment dynamics. Niu et
al. proposed generalized value iteration network (GVIN)[34] which is able to
learn and plan paths on irregular spatial graphs using a novel graph
convolution operator. Lee et al. proposed Gated Path Planning
Network(GPPN)[35] which replaces the unconventional recurrent update in VIN
with a gated LSTM recurrent operator to alleviate optimization issues like
random seed sensitivity. Experiment results on 2D and 3D path planning
demonstrate that GPPN outperforms VIN on accuracy, learning speed, training
instability and sensitivity to random seeds. Pflueger et al. proposed Rover-
IRL [10] based on VIN and inverse reinforcement learning(IRL) and achieves
good performance on planetary rover path planning problems.

Currently, existing Value Iteration based methods implement planning
computation without access to the environment dynamics $P_{ss^{\prime}}^{a}$,
which is not consistent with the Bellman equations. In this paper, a new Value
Iteration based path planning model called CIN is proposed, which utilizes
$P_{ss^{\prime}}^{a}$ and the sparse reward to implement VI algorithm with
convolution neural networks.

##  III Capability Iteration Network

Figure 3:  Structure of _Capability Iteration Network_ (CIN). The input of CIN
is the original map, current state $s_{t}$ and the goal position, which is
used to generate a sparse reward map firstly. For each state $s_{i}$, the
algorithm extracts its neighborhood and generates the state-conditioned
transition probability $P_{ss^{\prime}}^{a}$ utilizing the capability module.
Then the algorithm computes the Q-value map $Q(s,a)$ and the value map $V(s)$
utilizing convolution and max-pooling layers iteratively combined with the
Bellman equations. After $K$ times iteration, the best action $a^{*}$ is
obtained by max-pooling on the Q-value map: $a^{*}=\arg\max_{a}Q(s_{t},a)$.

###  III-A Framework

The framework of CIN is illustrated in Fig. 3 which consists of two important
modules: the capability module and value iteration module. The input of CIN is
a grid map $\mathcal{M}(m\times m)$, current state $s_{t}$ and the goal
$s_{g}$. A sparse reward map $R(s)$ is generated with positive reward $r_{p}$
at goal state and small negative living cost $r_{n}$ otherwise, i.e.

$R(s)=\left\\{\begin{aligned} r_{p}&,&s=s_{g},\\\ r_{n}&,&s\neq
s_{g}.\end{aligned}\right.$

The core component of CIN is a capability module where environment states are
processed as state-conditioned transition probabilities utilizing a neural
network classifier. The value iteration module computes the Q-value map
$Q(s,a)$ and value map $V(s)$ iteratively utilizing convolution layers with
kernels given by the capability module and max-pooling layers respectively.
The optimal Q-value map $Q^{*}(s,a)$ obtained after $K$ iterations is utilized
to find the best action $a^{*}$ for current state $s_{t}$ by maximization
through the action space.

###  III-B Capability Module

Capability module is the key component of CIN, which utilizes neural network
classifiers and environment information to predict probability distribution of
the next state $s_{t+1}$, i.e. $P(s_{t+1}|s_{t},a)$, often denoted as
$P_{ss^{\prime}}^{a}$. Actually, it is similar to the idea of the environment
dynamics utilized in the model based RL methods.

The input of the capability module is an $F\times F$ local map, which contains
information of state $s_{i,j}$ at the center and states can be reached in one
step from $s_{i,j}$. $F$ is the kernel size and a hyper-parameter for CIN,
which is related to the action capability of the agent. In maze environment,
$F=3$ is utilized. The output is the probability distribution of the next
state from the $s_{i,j}$ with all the available actions, i.e.

$\displaystyle f(\mathcal{M}_{[i,j,F]})_{k,l,a}$
$\displaystyle=P(s_{i+k,j+l}|s_{i,j,a}),$ $\displaystyle-\frac{2F-1}{2}\leq
k,l\leq\frac{2F-1}{2},a\in A$

where $f(\cdot)$ is the capability module, $X_{[i,j,F]}$ denotes the image
patch centered at position $(i,j)$ with kernel size $F$. Capability module
learns and encodes the capability of the agent in the learned neural network,
which is the only component needed to train in CIN and determines the
performance of the whole model.

For path planning tasks like 2D and 3D navigation, a multilayer perceptron
(MLP) classifier is enough; more complex classifiers like CNN[36] and
xgboost[37] are necessary for more difficult tasks. More details of the
capability module are provided in the experiment section.

###  III-C Value Iteration Module

The value iteration module utilizes convolution and max-pooling operators to
simulate the process of Value Iteration algorithm, which computes the Q-value
map $Q(s,a)$ and value map $V(s)$ iteratively. The following equations are the
Bellman equations, the basis of VI algorithm:

$Q(s,a)=\gamma\sum_{s^{\prime}}P_{ss^{\prime}}^{a}V(s^{\prime})+R(s,a)$ (4)
$V(s)=\max_{a}{Q(s,a)}$ (5)

As is shown in Fig. 3, CIN utilizes convolution to simulate the equation (4).
The transition probabilities $P_{ss^{\prime}}^{a}$ provided by the capability
module are the convolution kernels, and the sparse reward $R(s)$ are the bias,
i.e.

$Q^{(k+1)}(s_{i,j},a)=\gamma f(\mathcal{M}_{[i,j,F]})_{i,j}\cdot
V^{(k)}(s_{[i,j,F]})+R(s_{i,j}),\quad 1\leq i,j\leq m$

Then the new value map $V(s)$ is obtained utilizing max-pooling illustrated in
equation (5) and the Q-value map $Q(s,a)$ given by the previous convolution
operator:

$\displaystyle V^{(k+1)}(s_{i,j})=$ $\displaystyle\max\limits_{a\in
A}Q^{(k)}(s_{i,j},a),$ $\displaystyle\quad 1\leq i,j\leq m$

The VI algorithm tell us that the value iteration module will output the
optimal Q-value map $Q^{*}(s,a)$ after convergence:

$\lim_{k\to\infty}Q^{(k)}(s,a)=Q^{*}(s,a)$

In practical, the number of iteration $K$ is a hyper-parameter, which needs to
be large enough to ensure convergence. Then the best action $a^{*}$ for
current state $s_{t}$ is obtained utilizing maximization on the Q value map
through the action space:

$a^{*}=\arg\max_{a\in A}Q^{*}(s_{t},a)$

Compared to VIN model described in Fig. 2, CIN utilizes state-conditioned
transition probabilities $P_{ss^{\prime}}^{a}$ in the value iteration module
instead of $P_{s^{\prime}}^{a}$ in VIN, which is completely consistent with
the Bellman equation. This makes what CIN learns closer to the nature of the
task compared with VIN. The sparse reward map used in CIN consists with
positive reward at the goal and small living cost otherwise, which is simple
to generate and avoids the accuracy of CIN being influenced by the defects in
reward maps generated by the neural network, especially on large maps.

###  III-D Training Methods

There is no parameter to train in the value iteration module and the
capability module is just a simple classifier to predict the probability
distribution of the next state in MDP. Thus, we can train capability module
alone by supervised learning on local maps. Besides, training CIN end-to-end
is also a good choice.

End-to-end: The whole network of CIN is differentiable, which means we can
train CIN end-to-end using any RL or Imitation learning (IL)[38] algorithms.
The Cross Entropy Loss is utilized for IL in our experiment:

$\mathcal{L}_{CE}=-\sum^{|A|}_{i=1}p_{i}\log{\hat{p_{i}}},\quad
p_{i}=\left\\{\begin{aligned} 1&,&a_{i}=a^{*},\\\ 0&,&a_{i}\neq
a^{*}.\end{aligned}\right.$

where $\hat{p_{i}}$ is the probability of the model to choose the action
$a_{i}$.

CIN utilizes sparse reward maps immediately, rather than complex reward maps
generated by CNN in previous models. This makes CIN need less computation cost
and has faster learning speed compared to VIN, especially on large maps. The
experiment results show that CIN can achieve better accuracy and faster
learning speed compared to existing VI based models and reactive policies.
More details are described in the experiment section.

Train Capability Module Alone: We can also train capability module alone by
supervised learning. Take 2D grid world task as an example: the agent is
controlled with a random policy to obtain random trajectories. For each state
$s_{t}$ on the trajectories, $3\times 3$ small local map with center $s_{t}$
is sampled, which is utilized as the training data combined with the action
$a_{t}$ taken by the agent. The corresponding position of the next state
$s_{t+1}$ is the training label. The capability module is seen as a classifier
and trained utilizing supervised learning with the Mean Squared Error(MSE)
loss:

$\displaystyle\mathcal{L}_{MSE}=\frac{1}{F^{2}}\sum^{F}_{k=1}\sum^{F}_{l=1}(f(\mathcal{M}_{[i,j,F]})_{k,l,a}-P_{i,j,k,l})^{2}$

, where
$P_{i,j,k,l}=P(s_{i+k-\frac{2F-1}{2},\;j+l-\frac{2F-1}{2}}|s_{i,j},a)$.

To speed up convergence of the capability module, curriculum learning[39] is
needed. Take 3D terrain map task as an example, our training goal is to teach
capability module the capability of the agent to overcome obstacles (i.e., the
max height difference $\Delta h^{*}$ the agent can overcome). Thus, we start
training capability module by “simple” local maps with small or big height
difference. Gradually, “harder” local maps with height difference approach to
$\Delta h^{*}$ are needed for training. Training with curriculum learning can
speed up learning and achieves better accuracy.

##  IV Experiment and Dicussion

In this section, we evaluate the performance of CIN in three types of
experiment environments(2D, 3D grid world and real hexapod robot) and make a
comparison with existing methods including VIN, GPPN and ResNet-based reactive
policy to demonstrate the effectiveness and improvement of CIN.

###  IV-A 2D Grid-World Domain

| VIN | GPPN | CIN | ResNet18
---|---|---|---|---
$m$ | %Opt | %Suc | %Opt | %Suc | %Opt | %Suc | %Opt | %Suc
8 | 98.9 | 99.3 | 98.6 | 99.2 | 99.9 | 99.9 | 99.7 | 99.9
15 | 89.1 | 90.5 | 97.6 | 98.4 | 99.7 | 99.8 | 97.4 | 98.2
28 | 80.9 | 81.0 | 96.3 | 96.7 | 99.5 | 99.5 | 86.5 | 91.8
45 | 63.3 | 72.1 | 92.3 | 95.2 | 99.0 | 99.1 | 54.0 | 66.3
TABLE I: Experiment result in 2D grid world with varying map size ($m\times
m$).

(a) $28\times 28$ grid world (b) Reward map of CIN

(c) Value map of CIN

(d) Reward map of VIN (e) Value map of VIN

(f) Value map of VI

Figure 4: Visualization of the 2D maze, reward maps and value maps.

Grid world path planning with obstacles randomly located in the map is one of
the most common experiment environments for path planning problem. In 2D grid-
world domain, the input map contains binary numbers 0 (obstacles) and 1 (free-
space) only. The 2D grid world dataset is created with a maze generation
algorithm which uses Depth-First Search and the Recursive Backtracker
algorithm[40] with random starting and ending points. The size of the training
dataset $N_{tr}=10K$; validation dataset size $N_{val}=1K$; test dataset size
$N_{te}=1K$.

CIN: In this experiment, we train the capability module alone with supervised
learning for simplicity. The sparse reward map consists of a high reward +10
at the goal state and living cost -0.5 at other states. A multilayer
perceptron(MLP) consisting of four letent layers, one softmax layer and ReLU
activation function is utilized for the capability module in this task. We
utilize Mean Square Error (MSE) loss with Adam optimizer[41] for optimization.

Two baselines are utilized in our experiments: VIN[13]: VIN is a classic VI-
based model combined with VI algorithm and CNN. VIN achieves good performance
in path planning problems because of its planning computation. ResNet18[42]:
The problem setting for this task is similar to the image segmentation, in
which problem each pixel in the given image needs to be assigned a label (the
best action choice in our case). Thus, image segmentation model ResNet18 can
be seen as a reactive path planning model and utilized as a baseline in our
experiments.

We conduct experiments on maps with four different sizes: $8\times 8$,
$15\times 15$, $28\times 28$ and $45\times 45$. Two metrics are utilized to
evaluate the performance of algorithms, including %Optimal(%Opt)—the
percentage of states whose predicted paths under the policy estimated has
optimal length, %Success(%Suc)— the percentage of states whose predicted paths
under the policy estimated reach the goal state.

Table I shows the performance of models on each metric. Two conclusions can be
concluded: (1) CIN outperforms other models in 2D domain on two metrics,
especially on maps with large sizes. (2) Models with planning computation has
better performance on big maps, compared to reactive policies (ResNet18 in
this experiment). The conclusion (2) has been proposed in [13].

###  IV-B Model Result Visualization

The visualization of the model results are depicted in Fig. 4. The reward map
used in CIN (Fig. 4(b)) is sparse, which is consistent with the value
iteration algorithm. Reward map of VIN (Fig. 4(d)) is generated by CNN, which
is more messy, sensitive to the map settings and may influence the accuracy of
VIN on unseen or large mazes. Compared to VIN (Fig. 4(e)), value map of CIN
(Fig. 4(c)) is sharp at the edges of obstacles with obvious tendency towards
goal point, and is generally consistent with optimal value map (Fig. 4(f))
generated by the value iteration algorithm.

Figure 5: The learning curve of each model on the $15\times 15$ 2D mazes with
different random seeds. All the models are trained end-to-end using imitation
learning for fairness. CIN outperforms other models on learning speed,
training stability and sensitivity to random seeds.

###  IV-C Learning Speed and Stability

The learning curve of each model on the $15\times 15$ 2D mazes is depicted in
Fig. 5. In order to compare the learning speed of each model fairly, CIN is
trained end-to-end utilizing Imitation Learning(IL) in this experiment. All
the models are trained utilizing at least 5 different random seeds with
training dataset size $N_{tr}=10K$ in 50 epochs.

The %Error in Fig. 5 is the percentage of states whose predicted action by the
model estimated is different with the action given by the expert. As shown in
the figure, CIN and ResNet18 have faster learning speed compared to GPPN and
VIN, because CIN utilizes sparse reward immediately while VIN and GPPN uses
CNN to generate reward maps with more parameters to train. However, ResNet18
appears overfitted because it belongs to reactive policy without planning
computation. Besides, GPPN utilizes gated LSTM recurrent operator which
improves the training stability and random seed sensitivity compared to VIN.
As shown in the figure, CIN has better random seed sensitivity and converges
more stably than other models.

###  IV-D 3D Terrain Map

Unlike maps in 2D domain only contain 0 and 1, matrix elements in 3D terrain
maps can be any float value, representing height in a specific position. The
agent can only walk to its neighborhoods if the height difference is smaller
than a certain value. Other experiment settings are same to the 2D grid world
experiment. Compared to 2D mazes, 3D terrain maps are more complex and
challenging for path planning models.

| VIN | GPPN | CIN | ResNet18
---|---|---|---|---
$m$ | %Opt | %Suc | %Opt | %Suc | %Opt | %Suc | %Opt | %Suc
8 | 92.0 | 93.7 | 93.6 | 94.0 | 95.2 | 96.5 | 94.0 | 94.2
15 | 83.7 | 86.8 | 89.3 | 89.7 | 93.7 | 94.3 | 79.1 | 81.3
28 | 74.1 | 76.5 | 85.3 | 87.5 | 89.6 | 89.9 | 56.2 | 60.0
45 | 60.9 | 63.1 | 80.3 | 81.5 | 85.3 | 87.6 | 47.3 | 51.3
TABLE II: Experiment result in 3D grid world with varying map size ($m\times
m$).

Two metrics used in this experiment are same to which utilized in the 2D
environment. CIN is trained by training capability module alone while other
models are trained end-to-end with imitation learning. As shown in the Table.
II, CIN still achieves better performance than VIN and GPPN on various metrics
and map settings, especially on big maps.

###  IV-E Real Hexapod Robot

Hexapod robot. The hexapod robot has the characteristics of good stability,
large carrying capability and flexibility, and can adapt to complex
environments and special terrains. The robot used in this experiment is a
parallel structure robot, also known as Parallel-Parallel (PP) structure six-
legged robot. Each leg of the robot consists of three connecting rods
[$l_{1},l_{2},l_{3}$]. The lifting, lowering and stepping of the legs are
realized by the expansion and contraction of the three rods, and each leg can
be taken in multiple directions.

(a) Stairs (b) Uneven surface (c) U turn (d) Path planned

Figure 6: Simulation in CoppeliaSim (V-rep) with a hexapod robot. CIN achieves
good performance for the path planning on stairs, uneven terrain and U turn.

(a) Uneven surface (b) Narrow path (c) Left turn (d) Goal point

Figure 7: Experiment on the real hexapod robot. The same hexapod robot and
environment adaptation methods are utilized as in the simulation. CIN trained
in the simulation still achieves good performance in the real world. More
information is provided in the attached video.

Adaption to hexapod robot. Currently, VI-based algorithms are only effective
for discrete state and action space. Thus, simplification and discretization
are needed to be utilized. Given a real environment, we construct a 3D map
through point cloud data of Lidar with SLAM algorithm. Then the 2D map with
small squares of $0.2m\times 0.2m$ is obtained by projecting the 3D map to 2D
planes. Each small square is given a constant float number as the average
height. The 2D coordinate of the square which the robot mass center lies in is
utilized to represent the current state of the robot.

To simplify the large continuous action space of the hexapod robot, we design
a triangle-gait for the robot. In triangle-gait, six legs are divided into two
group with three non-adjacent feet as a group. The robot take three feet first
and then the other ones for a walk cycle. The gait is defined as (step length,
step height, rising length) with eight walking directions. Available
combinations are $(0.1m,0.3m,0.3m)$, $(0.2m,0.2m,0.2m)$, $(0.3m,0.1m,0.1m)$.
For simplicity, We train capability module alone to reduce the time cost.

We simulate the robot in CoppeliaSim (V-rep) which is a versatile, scalable
and powerful general-purpose robot simulation framework[43]. As is shown in
Fig. 6, the robot can step up stairs, make U-turn on the corner and handle
uneven surface, utilizing the path planned by CIN. Experiment in the real
environment is also conducted. As is shown in Fig. 7, we apply CIN trained in
the simulation to the real hexapod robot which also achieves good performance.
More details are provided in the attached video.

##  V Conclusion

In this paper, we propose a new VI-based path planning model, _Capability
Iteration Network_ , which improves accuracy and learning speed compared with
existing VI-based algorithms, especially on large-scale maps. CIN utilizes
sparse reward maps and state-conditioned transition probabilities
$P_{ss^{\prime}}^{a}$ instead of reward maps generated by CNN and global
convolution kernel $P_{s^{\prime}}^{a}$ in the previous algorithms. This makes
CIN consistent with VI algorithm, only need to learn the capability of the
agent with higher accuracy and faster learning speed. Several experiments are
conducted, including path planning on 2D, 3D grid world and a real hexapod
robot. The results demonstrate that CIN outperforms existing algorithms and is
able to be utilized in path planning tasks for real robots.

## References

  * [1] T. T. Mac, C. Copot, D. T. Tran, and R. De Keyser, “Heuristic approaches in robot path planning: A survey,” _Robotics and Autonomous Systems_ , vol. 86, pp. 13–28, 2016.
  * [2] E. Plaku, L. E. Kavraki, and M. Y. Vardi, “Motion planning with dynamics by a synergistic combination of layers of planning,” _IEEE Transactions on Robotics_ , vol. 26, no. 3, pp. 469–482, 2010.
  * [3] P. E. Hart, N. J. Nilsson, and B. Raphael, “A formal basis for the heuristic determination of minimum cost paths,” _IEEE transactions on Systems Science and Cybernetics_ , vol. 4, no. 2, pp. 100–107, 1968.
  * [4] O. Khatib, “Real-time obstacle avoidance for manipulators and mobile robots,” in _Autonomous robot vehicles_. Springer, 1986, pp. 396–404.
  * [5] S. M. LaValle, “Rapidly-exploring random trees: A new tool for path planning,” 1998.
  * [6] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, _et al._ , “Human-level control through deep reinforcement learning,” _Nature_ , vol. 518, no. 7540, p. 529, 2015.
  * [7] M. Sadeghzadeh, D. Calvert, and H. A. Abdullah, “Autonomous visual servoing of a robot manipulator using reinforcement learning,” _International Journal of Robotics and Automation_ , vol. 31, no. 1, 2016.
  * [8] T. Yan, W. Zhang, S. X. Yang, and L. Yu, “Soft actor-critic reinforcement learning for robotic manipulator with hindsight experience replay,” _International Journal of Robotics & Automation_, vol. 34, no. 5, pp. 536–543, 2019.
  * [9] Y. Liu, M. Cong, H. Dong, and D. Liu, “Reinforcement learning and ega-based trajectory planning for dual robots,” _International Journal of Robotics & Automation_, vol. 33, no. 4, pp. 367–378, 2018.
  * [10] M. Pflueger, A. Agha, and G. S. Sukhatme, “Rover-irl: Inverse reinforcement learning with soft value iteration networks for planetary rover path planning,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 1387–1394, 2019.
  * [11] A. Faust, K. Oslund, O. Ramirez, A. Francis, L. Tapia, M. Fiser, and J. Davidson, “Prm-rl: Long-range robotic navigation tasks by combining reinforcement learning and sampling-based planning,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 5113–5120.
  * [12] K. Arulkumaran, M. P. Deisenroth, M. Brundage, and A. A. Bharath, “A brief survey of deep reinforcement learning,” _arXiv preprint arXiv:1708.05866_ , 2017.
  * [13] A. Tamar, Y. Wu, G. Thomas, S. Levine, and P. Abbeel, “Value iteration networks,” in _Advances in Neural Information Processing Systems_ , 2016, pp. 2154–2162.
  * [14] R. Bellman, “Dynamic programming,” _Science_ , vol. 153, no. 3731, pp. 34–37, 1966.
  * [15] M. Andrychowicz, F. Wolski, A. Ray, J. Schneider, R. Fong, P. Welinder, B. McGrew, J. Tobin, O. P. Abbeel, and W. Zaremba, “Hindsight experience replay,” in _Advances in neural information processing systems_ , 2017, pp. 5048–5058.
  * [16] M. Vecerik, T. Hester, J. Scholz, F. Wang, O. Pietquin, B. Piot, N. Heess, T. Rothörl, T. Lampe, and M. Riedmiller, “Leveraging demonstrations for deep reinforcement learning on robotics problems with sparse rewards,” _arXiv preprint arXiv:1707.08817_ , 2017.
  * [17] M. Riedmiller, R. Hafner, T. Lampe, M. Neunert, J. Degrave, T. Wiele, V. Mnih, N. Heess, and J. T. Springenberg, “Learning by playing solving sparse reward tasks from scratch,” in _International Conference on Machine Learning_ , 2018, pp. 4344–4353.
  * [18] O. Souissi, R. Benatitallah, D. Duvivier, A. Artiba, N. Belanger, and P. Feyzeau, “Path planning: A 2013 survey,” in _Proceedings of 2013 International Conference on Industrial Engineering and Systems Management (IESM)_. IEEE, 2013, pp. 1–8.
  * [19] B. Paden, M. Čáp, S. Z. Yong, D. Yershov, and E. Frazzoli, “A survey of motion planning and control techniques for self-driving urban vehicles,” _IEEE Transactions on intelligent vehicles_ , vol. 1, no. 1, pp. 33–55, 2016\.
  * [20] G. Che, L. Liu, and Z. Yu, “An improved ant colony optimization algorithm based on particle swarm optimization algorithm for path planning of autonomous underwater vehicle,” _Journal of Ambient Intelligence and Humanized Computing_ , vol. 11, no. 8, pp. 3349–3354, 2020.
  * [21] P. K. Y. Yap, N. Burch, R. C. Holte, and J. Schaeffer, “Any-angle path planning for computer games,” in _Seventh Artificial Intelligence and Interactive Digital Entertainment Conference_ , 2011.
  * [22] W. Gong, X. Xie, and Y.-J. Liu, “Human experience–inspired path planning for robots,” _International Journal of Advanced Robotic Systems_ , vol. 15, no. 1, p. 1729881418757046, 2018.
  * [23] S. Koenig and M. Likhachev, “D* lite,” in _Eighteenth national conference on Artificial intelligence_ , 2002, pp. 476–483.
  * [24] J. Wang, X.-j. Lin, H.-y. Zhang, G.-d. Lu, Q.-l. Pan, and H. Li, “Path planning of manipulator using potential field combined with sphere tree model,” _International Journal of Robotics and Automation_ , vol. 35, no. 2, 2020.
  * [25] L. E. Kavraki, P. Svestka, J.-C. Latombe, and M. H. Overmars, “Probabilistic roadmaps for path planning in high-dimensional configuration spaces,” _IEEE transactions on Robotics and Automation_ , vol. 12, no. 4, pp. 566–580, 1996.
  * [26] Y. Liu, M. Cong, H. Dong, D. Liu, and Y. Du, “Time-optimal motion planning for robot manipulators based on elitist genetic algorithm,” _International Journal of Robotics & Automation_, vol. 32, no. 4, pp. 396–405, 2017.
  * [27] J. Xin, J. Zhong, J. Sheng, P. Li, and Y. Cui, “Improved genetic algorithms based on data-driven operators for path planning of unmanned surface vehicle,” _International Journal of Robotics and Automation_ , vol. 34, no. 6, 2019.
  * [28] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” _arXiv preprint arXiv:1509.02971_ , 2015.
  * [29] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. Riedmiller, “Playing atari with deep reinforcement learning,” _arXiv preprint arXiv:1312.5602_ , 2013.
  * [30] Y. Duan, X. Chen, R. Houthooft, J. Schulman, and P. Abbeel, “Benchmarking deep reinforcement learning for continuous control,” in _International Conference on Machine Learning_ , 2016, pp. 1329–1338.
  * [31] M. Deisenroth and C. E. Rasmussen, “Pilco: A model-based and data-efficient approach to policy search,” in _Proceedings of the 28th International Conference on machine learning (ICML-11)_ , 2011, pp. 465–472.
  * [32] R. McAllister and C. E. Rasmussen, “Data-efficient reinforcement learning in continuous state-action gaussian-pomdps,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 2040–2049.
  * [33] D. Ha and J. Schmidhuber, “Recurrent world models facilitate policy evolution,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 2450–2462.
  * [34] S. Niu, S. Chen, H. Guo, C. Targonski, M. C. Smith, and J. Kovačević, “Generalized value iteration networks: Life beyond lattices,” in _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018.
  * [35] L. Lee, E. Parisotto, D. S. Chaplot, E. Xing, and R. Salakhutdinov, “Gated path planning networks,” in _International Conference on Machine Learning_ , 2018, pp. 2947–2955.
  * [36] Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, and L. D. Jackel, “Backpropagation applied to handwritten zip code recognition,” _Neural computation_ , vol. 1, no. 4, pp. 541–551, 1989.
  * [37] T. Chen and C. Guestrin, “Xgboost: A scalable tree boosting system,” in _Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining_ , 2016, pp. 785–794.
  * [38] P. Abbeel and A. Y. Ng, “Apprenticeship learning via inverse reinforcement learning,” in _Proceedings of the twenty-first international conference on Machine learning_ , 2004, p. 1.
  * [39] Y. Bengio, J. Louradour, R. Collobert, and J. Weston, “Curriculum learning,” in _Proceedings of the 26th annual international conference on machine learning_. ACM, 2009, pp. 41–48.
  * [40] Wikipedia, “Maze generation algorithm — Wikipedia, the free encyclopedia,” 2020, [Online; accessed 01-September-2020].
  * [41] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [42] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [43] E. Rohmer, S. P. Singh, and M. Freese, “V-rep: A versatile and scalable robot simulation framework,” in _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2013, pp. 1321–1326.


# Engineering bound states in the continuum at telecom wavelengths with non-
Bravais lattices

Shunsuke Murai  Department of Material Chemistry, Graduate School of
Engineering, Kyoto University, Nishikyo-ku, Kyoto 615-8510, Japan
<EMAIL_ADDRESS>Diego R. Abujetas  Physics Department,
Fribourg University, Chemin de Musée 3, 1700 Fribourg Switzerland  Libei Liu
Department of Material Chemistry, Graduate School of Engineering, Kyoto
University, Nishikyo-ku, Kyoto 615-8510, Japan  Gabriel W. Castellanos
Department of Applied Physics and Institute for Photonic Integration,
Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The
Netherlands  Vincenzo Giannini  Instituto de Estructura de la Materia (IEM-
CSIC), Consejo Superior de Investigaciones Científicas, Serrano 121, 28006
Madrid, Spain<EMAIL_ADDRESS>José A. Sánchez-Gil  Instituto de Estructura
de la Materia (IEM-CSIC), Consejo Superior de Investigaciones Científicas,
Serrano 121, 28006 Madrid, Spain<EMAIL_ADDRESS>Katsuhisa Tanaka
Department of Material Chemistry, Graduate School of Engineering, Kyoto
University, Nishikyo-ku, Kyoto 615-8510, Japan  Jaime Gómez Rivas  Department
of Applied Physics and Institute for Photonic Integration, Eindhoven
University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands
<EMAIL_ADDRESS>
###### Abstract

Various optical phenomena can be induced in periodic arrays of nanoparticles
by the radiative coupling of the local dipoles in each particle. Probably the
most impressive example is bound states in the continuum (BICs), which are
electromagnetic modes with a dispersion inside the light cone but infinite
lifetime, i.e., modes that cannot leak to the continuum. Symmetry-protected
BICs appear at highly symmetric points in the dispersion of periodic systems.
Although the addition of nonequivalent lattice points in a unit cell is an
easy and straightforward way of tuning the symmetry, BICs in such particle
lattice, i.e., non-Bravais lattice, are less explored among periodic systems.
Starting from a periodic square lattice of Si nanodisks, we have prepared
three non-Bravais lattices by detuning size and position of the second disk in
the unit cell. Diffraction-induced coupling excites magnetic/electric dipoles
in each nanodisk, producing two surface lattice resonances at the $\Gamma$
point with a band gap in between. The high/low energy branch becomes a BIC for
the size/position-detuned array, respectively, while both branches are bright
(or leaky) when both size and position are detuned simultaneously. The role of
magnetic and electric resonances in dielectric nanoparticles and the change of
BIC to bright character of the modes is explained by the two different origins
of BICs in the detuned arrays, which is further discussed with the aid of a
coupled electric and magnetic dipole model. This study gives a simple way of
tuning BICs at telecom wavelengths in non-Bravais lattices, including both
plasmonic and dielectric systems, thus scalable to a wide range of
frequencies.

###### keywords:

American Chemical Society, LaTeX

Department of Applied Physics and Institute for Photonic Integration,
Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The
Netherlands Centre of Excellence ENSEMBLE3 sp. z o.o., Wolczynska 133, Warsaw,
01-919, Poland Technology Innovation Institute, Masdar City 9639, Abu Dhabi,
United Arab Emirates IR,NMR,UV

##  1 Introduction

Periodic arrays of resonant particles represent an arena of intensive research
in modern optics. Localized resonances at each particle can be radiatively
coupled via in-plane diffraction, causing hybridized modes that are called
surface lattice resonances (SLRs)1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12. SLRs
are characterized by narrow and strong resonances that appear in a wide
spectrum range from terahertz13 to ultraviolet14, depending on the periodicity
of the lattice and the size of the particles. The arrays include one
dimensional chains3, 2, 15 and two dimensional lattices of particles. The
most-studied two dimensional lattices are square and hexagonal patterns, which
fall into two of the five Bravais lattices in two dimensions. In a Bravais
lattice, all the particles are crystallographically equivalent and can be
generated by a set of translation operations. Some works have explored non-
Bravais lattices, where a second set of basis is introduced to the initial
Bravais lattice16, 17, 18, 19, 20, 21, 22, 23, 24, 25. This provides the
nonequivalent lattice points in the unit cell, giving another channel for
particle interaction. In many cases the second basis lowers the lattice
symmetry, which changes the modes that can be excited in the system.

One related phenomenon that appears in periodic lattices is bound states in
the continuum (BICs), which are bounded (non-radiative or dark) eigenstates of
an optical system despite having a dispersion within the light cone of the
surrounding medium 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45. The quality ($Q$) factor of a BIC approaches infinity in
the absence of material losses, which is of particular interest for lasing and
sensing applications29, 46, 47, 48, 49. Symmetry-protected BICs often appear
in high-symmetry points of the first Brillouin zone, such as at the $\Gamma$
point. Defects and perturbations breaking the structural symmetries allow
coupling of BICs to the radiative modes, transforming them into leaky
resonances or quasi-BICs with a finite $Q$-factor.

Since the addition of a second basis lowers the lattice symmetry, non-Bravais
lattices should represent a rich field for investigating and controlling BICs.
BICs in a non-Bravais lattice have been recently reported in the terahertz
(THz) and visible regions using metallic systems. In the THz region, a pair of
gold resonant bars (dimer) sustain anti-parallel resonances that cancel each
other, forming a BIC at the $\Gamma$ point50. Such dimers have been also
studied at optical frequencies, where pairs of silver nanoparticles in a
lattice sustain a BIC mode17, 16, 20. These works have focused on the
interaction between pairs of electric dipoles, leaving unexplored the
interplay with magnetic dipoles and the relation between the BICs and the
symmetry of the unit cell.

In this study, we unveil BICs emerging in non-Bravais rectangular lattices
with different unit cell symmetries. The array consists of Si nanodisks which
support both (in-plane and out-of-plane) electric and magnetic dipoles at
telecom wavelengths. Therefore, the results presented here can be easily
adaptable to any other dipolar lattice in contrast to metallic systems where
only electric dipoles are excitable. The emergence of BICs is explained in an
intuitive way based on the symmetry matching between the lattice and the pair
of dipoles constituting the mode, which is further verified by a theoretical
model based on coupled electric and magnetic dipoles (CEMD)51, 52, 53.

##  2 non-Bravais lattice with size and position detunings

We design the non-Bravais lattices that support BICs in the telecom region
starting from a square array of Si nanodisks. The array consists of
polysilicon nanodisks (200 nm height, 224 nm diameter) arranged in a square
lattice with a period of $a_{\rm x}=a_{\rm y}$ = 425 nm (see Fig.1(a)). This
is the initial Bravais lattice structure that we call symmetric array. Note
that the period is too short to sustain SLRs in the targeted telecom region.
In the space group terminology, the symmetric array is classified as P4mm.
Then we lower the symmetry of the lattice by shifting the position of the
second disk by $\delta$x = 10 $\%$ from the original position (Fig.1(b)),
while the diameters of the nanodisks are the same (236 nm). The size of the
unit cell is doubled by this operation ($a_{x}$ = 850 nm) and the second disk
lies at 2$a_{x}$/5 = 340 nm in the unit cell. It is noted that the two disks
in the unit cell are not equivalent: i.e., the lattice is non-Bravais. We call
this array the position-detuned array (Fig.1(b)). The space group of the
position-detuned array is Pmm2, where a mirror plane lies at the center of the
unit cell along the $y$-direction as indicated by a dotted line in the
scanning electron microscopic (SEM) image. Another way of lowering the
symmetry is to detune the size of the disk. The diameter of the second
nanodisk is thus increased by $\delta D$ = 31% (256 nm) with respect to the
first disk (194 nm). This operation also makes $a_{x}$ = 850 nm and we call it
the size-detuned array (Fig.1(c)). This array also falls into Pmm2 in space
group, i.e., the position- and size-detuned arrays have the same elements in
terms of the symmetry. We also make a doubly-detuned array, where both the
position and the size of the second disk are detuned from the first disk by
$\delta$x = 10% and $\delta D$ = 31%, respectively. The sample is called
doubly-detuned array and the diameters of the nanodisks are 210 and 276 nm
(Fig.1(d)). The space group of the doubly-detuned array is Pm; i.e., the
mirror planes along the $y$-direction are lost in this array.

##  3 Results and Discussion

###  3.1 Mode dispersion at telecom frequencies

We measure the zeroth-order (forward) transmission spectrum by illuminating
the sample with a collimated beam from a supercontuniuum laser. To realize a
homogeneous medium surrounding the nanodisks that provides out-of-plane
symmetry, we place a quartz coverslip on top of the nanodisk array with
refractive index matching oil ($\varepsilon=2.1$) in between. We evaluate the
optical extinction ($E$), defined as $E=1-T/T_{0}$, where $T$ is the
transmission through the sample and $T_{0}$ is the transmission through a
reference consisting of a substrate, index-matching oil and superstrate. The
extinction spectra are shown for TE-polarization in Figs.1(e)-(h) as a
function of the angle of incidence, $\theta_{\rm in}$, where the incident
plane is defined in the $z-x$ plane.

For the symmetric array (Fig.1(e)), no extinction features appear in the range
of measurement ($\lambda$ between 1100 and 1500 nm) because the lattice
resonances occur in the visible frequency range. On the other hand, for the
position-detuned array (Fig.1(f)), highly dispersive features appear following
the in-plane diffraction orders (Rayleigh anomalies), which correspond to the
SLRs. The Rayleigh anomalies satisfy

$\vec{k}_{0}=\vec{k}_{x}(\theta_{\rm{in}})\pm\vec{G}(m_{1},m_{2}),$ (1)

where $\vec{k}_{0}$ is the wave vector of the in-plane diffracted orders,
$\vec{k_{x}}$ is the in-plane component of the wave vector of the incident
beam, and
$\vec{G}(m_{1},m_{2})=m_{1}(\frac{2\pi}{a_{x}})\vec{x}+m_{2}(\frac{2\pi}{a_{y}})\vec{y}$
is the reciprocal lattice vector of the array with $m_{1}$ and $m_{2}$ the
orders of diffraction in the $x-$ and $y-$ directions, respectively. At the
$\Gamma$ point, the SLRs bend with a splitting energy of $\sim$ 41 meV between
the upper and lower modes. The longer wavelength or lower-energy SLR becomes
narrower as $\theta_{\rm in}$ approaches 0∘, leading to a BIC at $\theta_{\rm
in}$ = 0∘ as seen in Fig.1(f). SLRs also appear for the size-detuned array
(Fig.1(g)), and interestingly enough, the BIC mode corresponds to the higher-
energy SLR. For the doubly-detuned array (Fig. 1(h)), both SLRs are bright.
The cuts to the extinction dispersion measurements further show the evolution
of the SLRs around the $\Gamma$ point (Figs. 1(i-l)).

The emergence of these BICs at the $\Gamma$ point can be understood in an
intuitive way upon considering the modes supported by the symmetric array and
how they evolve as the symmetry of the lattice is changed, as depicted in
Figs. 2(a-b). In the frequency range of interest, the symmetric array supports
two guided modes given by the out-of-phase (anti-parallel) oscillation of in-
plane electric dipoles (EDy) and out-of-plane magnetic dipoles (MDz). We will
refer to them as the electric and magnetic guided modes, respectively. Since
guided modes cannot be excited with propagating plane waves, the extinction
spectrum shows no features (see Figs. 1(e) and (i)). In contrast, for the
other arrays, the perturbation of the unit cell doubles their size along the
$x$ axis and the dispersion relations of the guided modes fold down to the
continuum of radiation. Therefore, the initially guided modes of the symmetric
array evolve to leaky (bright) resonances or even to BICs.

For the position-detuned array, the displacement of the second disk in each
(new) unit cell allows the coupling between EDy and MDz, forming hybrid modes.
At the $\Gamma$ point, the eigenstate stemming from the magnetic (electric)
guided mode is given by anti-parallel oscillations of MDz (EDy) mixed with
parallel oscillations of EDy (MDz), as depicted in Fig. 2(a) (2(b)). Note that
both, EDy and MDz, oscillating parallel or anti-parallel are not self-
consistent dipole configurations in the position-detuned array. Since the MDz
does not radiate in any case at the $\Gamma$ point54 and the anti-parallel EDy
cancels with each other, the electric guided mode evolves to a BIC (Fig.
2(b)). In contrast, the magnetic guided mode evolves to a bright mode, due to
the parallel oscillations of the EDy associated to this eigenmode (Fig. 2(a)).
On the other hand, for the size-detuned array, EDy and MDz do not mix, but the
relative strength between the dipoles in each cell is different. Therefore,
the magnetic guided mode evolves to a BIC (MDz does not radiate at the
$\Gamma$ point), while the electric guided mode evolves to a bright resonance,
since the unequal amplitudes of the EDy break the symmetry at normal incidence
as seen in Figs. 2(a) and (b), respectively. Finally, for the doubly-detuned
array, both resonances can couple to a plane wave at the $\Gamma$ point, so
that both become bright. This is understood by combining the phenomenology
exhibited for the position- and size-detuned arrays. The absence of symmetry-
protected BICs is in accordance with the fact that the doubly-detuned array
shows the lowest symmetry among the arrays in the present study.

The cuts to extinction spectra at normal incidence are shown in Figs. 2(c-f),
together with the schematic of the modes supported by each lattice. The
phenomenology described in terms of the eigenmodes of the systems adequately
suits with the experimental results: The modes related to the magnetic guided
mode are always at higher energies than those derived from the electric guided
mode. Also, the mode frequency is only slightly altered upon a lattice
perturbation (note that BIC frequencies can be estimated from the extinction
dispersion spectra), suggesting that the framework depicted in Figs. 2(a) and
(b) is correct.

For the sake of completeness, the extinction for TM light is shown in
Supporting Fig. S1. The two SLRs also appear for TM polarization, but the
splitting is smaller ($\sim$ 21 meV). For this polarization, in-plane MDs
instead of EDs couple to the diffraction order55, 56, 57, and a smaller
splitting indicates a smaller in-plane magnetic polarizability of the Si
nanodisk compared to the electric polarizability. For the position-detuned
array, the BIC from anti-parallel MDy mixed with parallel EDz is supported,
while anti-parallel EDz creates the BIC for the size-detuned array.

###  3.2 CEMD calculation

To confirm our intuitive explanation of the underlying physical mechanism, we
have investigated theoretically the extinction of the arrays using a CEMD
model. 58, 52 We extract numerically the polarizability of electric and
magnetic dipoles from the single disks, and calculate analytically $1-T$
(equivalent to the experimental extinction, since $T$ is the specular
transmittance) when the dipoles are arranged in an infinite lattice. The
calculations shown in Figs. 3(a-d) reproduce fairly well the experimental
results in Figs. 1(e-h). No modes appear for the symmetric array as in the
experiment. The flip of energy of the bright and BIC modes in position- and
size-detuned arrays is clearly seen in Figs. 3(b) and (c), and the two bright
modes at the $\Gamma$ point are obvious for the doubly-detuned array (Fig.
3(d)).

We further verify the intuitive assignment of the modes in Fig.2 by projecting
the total calculated extinction into the in-plane and out-of-plane ED and MD
components. The bottom panels in Fig. 3 show the cuts to the extinction
spectra around $\theta_{\rm in}$ = 0∘. The color of the lines for all three
detuned arrays in Figs. 3(b)-(d) reflects the origin of the modes, in
accordance with the assignments in Fig. 2 for the lower and higher energy
bands: red represents the projection over the low-energy state with anti-
parallel EDy mixed with parallel MDz, while blue denotes the projection over
the high-energy state with anti-parallel MDz mixed with parallel EDy. This
projection is done because these states are the theoretical eigenstates of the
system at normal incidence. The brightness of the mode is determined by the
position and magnitude of the dipoles, which are governed by the position and
size of the disks, respectively, as explained in Fig. 2. The projections of
the theoretical results support the above-mentioned nature of the bands and
the emergence of BICs at the $\Gamma$ point in Figs. 3(f) and (g), only
plausible with combinations of anti-parallel (identical) in-plane ED and/or
arbitrary vertical MDs. Away from the $\Gamma$ point the SLRs bands may
exhibit more complex combinations of dipolar states, as revealed in Figs.
3(f-h) by the mixture of red/blue colors. We also plot separately the
contributions from the EDy and MDz components to the total extinction as a
function of frequency and angle of incidence in Supporting Fig. S2; note that
the combination of both in-plane ED and out-of-plane MD contributions are
needed to reproduce the discussed SLR bands and resulting BICs.

The direct evidence of the transition from quasi-BIC to symmetry-protected BIC
is the resonance narrowing as it approaches $\theta_{\rm in}$ = 0 ∘. We
extract the $Q$-factors from the two quasi-BICs appearing for the position-
and size-detuned arrays in Fig. 4. For the position-detuned array in Fig.
4(a), the $Q$-factor is so large that the value is above the resolution limit
of the spectrometer (indicated by the dotted line in the plot) for all the
measured angles. We also have plotted the $Q$-values extracted from the CEMD
calculation by solid lines, showing a diverging behaviour as $\theta_{\rm in}$
approaches 0∘. For the size-detuned array (Fig. 4(b)), the $Q$-factor
increases as $\theta_{\rm in}$ approaches 0 ∘ and reaches the detection limit
imposed by the resolution of the spectrometer. The agreement between the
experiment and calculation is excellent. This result might be surprising given
the rough surface of the nanoparticles as observed in the SEM images shown in
Fig. 1. One possible reason for the excellent agreement between experiments
and calculations is the working wavelength in the present work: Telecom
wavelengths are relatively long compared to the surface imperfections and
scattering can be neglected. Also the absorption of Si is negligible in this
wavelength range, allowing the disks to be regarded as pure magnetic and
electric dipoles.

###  3.3 Spectral tuning of the BIC position

In order to demonstrate the universality of the BICs appearing in rectangular
non-Bravais lattice, we change the parameters in the design to shift the
features. When the size of nanodisks is enlarged while the lattice constants
remain the same, the two bands redshift with a larger splitting ($\sim$ 83 meV
for TE light) due to the larger polarizability of the disks. Details are found
in Supporting Sec. 3.

We further shift the spectral features to the visible range by scaling the
structure. The prepared size-detuned array consists of polysilicon nanodisks
(86 and 107 nm diameters) arranged with $a_{x}$ = 400 nm and $a_{y}$ = 200 nm
(see the inset in Fig. 5(c) for the SEM image). The extinction shows a
localized dispersion-less band at $\lambda$ $\leq$ 500 nm and the two SLRs
(Fig. 5). At $\theta_{\rm in}$=0∘, the high and low-energy bands appear with
an energy gap of $\sim$ 105 meV for TE light (Fig. 5(a)), and the higher-
energy SLR becomes a BIC. For TM polarization (Fig. 5(b)), the BIC mode is
flipped to the lower energy band although less obvious due to a smaller
magnetic polarizability. The appearance of BIC is evidenced by the $Q$-factor
as a function of $\theta_{\rm in}$, where $Q$ diverges as $\theta_{\rm in}$
approaches 0∘ in Fig. 5(c).

##  4 Conclusions

We have demonstrated experimentally the emergence of symmetry-protected BICs
at telecom wavelengths by detuning square arrays of polysilicon nanodisks into
rectangular non-Bravais lattices. The initial square array supports two guided
modes, comprising anti-parallel EDs and MDs, respectively, that turn into SLRs
upon doubling the unit cell length in one direction by modifying the
position/size of every other disk. For the position-detuned array, the TE (TM)
mode polarized along the short (long) axis of the lattice, consisting of
horizontal anti-parallel EDs (MDs) mixes with vertical parallel MDs (EDs),
becoming a BIC; for the size-detuned array, conversely, vertical anti-parallel
MDs (EDs) lead to a BIC for the TE (TM) mode. When the symmetry of the lattice
is further lowered by detuning size and position of the second nanodisks
simultaneously, the symmetry-protected BICs disappear. The measurements of the
$Q$-factor of the BIC is limited by the resolution of the spectrometer, being
this $Q$-factor larger than 400. This work provides a robust design strategy
of symmetry-protected BICs via introduction of a second meta-atom (nanodisk in
this case) per unit cell to Bravais lattices, which is readily applicable to
lasing, sensing, and light-matter coupling researches.

##  5 Methods

Fabrication. Polycrystalline Si thin films with a thickness of 80 nm were
grown on a synthetic silica glass substrate by low-pressure chemical vapor
deposition using SiH4 gas as a source of Si. A resist (NEB22A2, Sumitomo) was
cast on the Si film and exposed to electron-beam lithography, followed by
development to make nanoparticle arrays of resist on the Si film. The Si film
was vertically etched using a selective dry etching (Bosch process) with SF6
and C4H8 gases, and the resist residue was etched away by oxygen dry etching.
The fabricated area of the array covers 3 $\times$ 3 mm2.

Extinction measurements. The nanodisk arrays were in an optically homogeneous
environment by placing an upperstrate of synthetic silica glass with index-
matching oil. The samples were placed in a rotational stage, and illuminated
by a collimated and polarized beam. The source of illumination was a broadband
supercontinuum laser (FIANIUM WHITELASE MICRO). The beam diameter was $0.1$
mm, smaller than the area of the array. The optical axis corresponds to the
$z-$axis and the sample was rotated around the $y-$axis (see Fig. 1a). The
polarization of the incident light was fixed either along the $y$-axis (TE
polarization) or the $x$-axis (TM polarization). As the sample was rotated
along the $y$-axis by an angle $\theta_{\rm in}$, light was incident with a
wave vector component that was parallel to the array plane along the
$x$-direction, $k_{x}=\frac{2\pi n}{\lambda}\sin{\theta_{\rm in}}$, where
$\lambda$ is the free-space wavelength and $n$ is refractive index of the
glass. We measured the spectra with a fiber-coupled spectrometer (NIR Quest,
Ocean Optics) with 3.1 nm spectral resolution. We calculated the extinction as
$E=1-T/T_{0}$, where $T$ and $T_{0}$ are the zeroth-order transmission
measured on and off the array, respectively.

CEMD theoretical calculations. We make use of our coupled electric/magnetic
dipole theory (CEMD) for infinite 2D planar arrays developed in Refs. 51, 52,
53. In the case of disk meta-atoms with axes perpendicular to the metasurface,
as explained in Ref. 59, we need the polarizabilities of the in-plane and out-
of-plane, electric and magnetic resonances to apply the CEMD theory. These
polarizabilities are numerically obtained through SCUFF 60, 61, an open-source
software package for analysis of electromagnetic scattering problems using the
method of moments, the dielectric function of polycrystalline silicon
$\varepsilon=3.5$,62, 63, 64 and assuming a medium surrounding the nanodisks
with $\varepsilon=2.1$. The four calculated polarizibilites are introduced in
our CEMD model, i.e., in- and out-of-plane MD and ED, which accounts for all
lattice-induced interactions.

This work was partly supported by the Nanotechnology Hub, Kyoto University and
Kitakyusyu FAIS in the “Nanotechnology Platform Project,” sponsored by MEXT,
Japan. We gratefully acknowledge the financial support from Kakenhi (17KK0133,
19H02434) and the Netherlands Organisation for Scientific Research (NWO)
through Gravitation grant “Research Centre for Integrated Nanophotonics” and
Innovational Research Activities Scheme (Vici project SCOPE no. 680-47-628).
J.A.S.-G. and D.R.A. acknowledge support from the Spanish Ministerio de
Ciencia e Innovación through grants MELODIA PGC2018-095777-B-C21 and FPU PhD
Fellowship FPU15/03566 (MCIU/AEI/FEDER, UE).

Supporting Information is available: Extinction for the TM-polarized light;
Contributions of in- and out-of-plane components to the extinction; BICs in
the size-detuned array with larger disk diameters.

## References

  * Zou et al. 2004 Zou, S.; Janel, N.; Schatz, G. C. Silver nanoparticle array structures that produce remarkably narrow plasmon lineshapes. _J. Chem. Phys._ 2004, _120_ , 10871
  * Hicks et al. 2005 Hicks, E. M.; Zou, S.; Schatz, G. C.; Spears, K. G.; Van Duyne, R. P.; Gunnarsson, L.; Rindzevicius, T.; Kasemo, B.; Käll, M. Controlling Plasmon Line Shapes Through Diffractive Coupling in Linear Arrays of Cylindrical Nanoparticles Fabricated by Electron Beam Lithography. _Nano Lett._ 2005, _5_ , 1065
  * Zou and Schatz 2006 Zou, S.; Schatz, G. C. Theoretical Studies of Plasmon Resonances in One-Dimensional Nanoparticle Chains: Narrow Lineshapes with Tunable Widths. _Nanotechnology_ 2006, _17_ , 2813
  * Chu et al. 2008 Chu, Y. Z.; Schonbrun, E.; Yang, T.; Crozier, K. B. Experimental observation of narrow surface plasmon resonances in gold nanoparticle arrays. _Applied Physics Letters_ 2008, _93_ , 181108
  * Auguie and Barnes 2008 Auguie, B.; Barnes, W. L. Collective resonances in gold nanoparticle arrays. _Physical Review Letters_ 2008, _101_ , 143902
  * Kravets et al. 2008 Kravets, V. G.; Schedin, F.; Grigorenko, A. N. Extremely Narrow Plasmon Resonances Based on Diffraction Coupling of Localized Plasmons in Arrays of Metallic Nanoparticles. _Phys. Rev. Lett._ 2008, _101_ , 087403
  * Vecchi et al. 2009 Vecchi, G.; Giannini, V.; Gómez Rivas, J. Shaping the Fluorescent Emission by Lattice Resonances in Plasmonic Crystals of Nanoantennas. _Phys. Rev. Lett._ 2009, _102_ , 146807
  * Vecchi et al. 2009 Vecchi, G.; Giannini, V.; Gómez Rivas, J. Surface Modes in Plasmonic Crystals Induced by Diffractive Coupling of Nanoantennas. _Phys. Rev. B: Condens. Matter Mater. Phys._ 2009, _80_ , 201401
  * Zhou and Odom 2011 Zhou, W.; Odom, T. W. Tunable Subradiant Lattice Plasmons by Out-of-plane Dipolar Interactions. _Nat. Nanotechnol._ 2011, _6_ , 423
  * Collin 2014 Collin, S. Nanostructure arrays in free-space: optical properties and applications. _Reports on Progress in Physics_ 2014, _77_ , 126402
  * Wang et al. 2018 Wang, W.; Ramezani, M.; Väkeväinen, A. I.; Törmä, P.; Rivas, J. G.; Odom, T. W. The rich photonic world of plasmonic nanoparticle arrays. _Mater. Today_ 2018, _21_ , 303–314
  * Kravets et al. 2018 Kravets, V. G.; Kabashin, A. V.; Barnes, W. L.; Grigorenko, A. N. Plasmonic surface lattice resonances: A review of properties and applications. _Chemical Reviews_ 2018, _118_ , 5912–5951
  * Schaafsma et al. 2016 Schaafsma, M. C.; Bhattacharya, A.; Rivas, J. G. Diffraction Enhanced Transparency and Slow THz Light in Periodic Arrays of Detuned and Displaced Dipoles. _ACS Photonics_ 2016, _3_ , 1596–1603
  * Kawachiya et al. 2018 Kawachiya, Y.; Murai, S.; Saito, M.; Sakamoto, H.; Fujita, K.; Tanaka, K. Collective plasmonic modes excited in Al nanocylinder arrays in the UV spectral region. _Optics Express_ 2018, _26_ , 5970–5982
  * Pocock et al. 2018 Pocock, S. R.; Xiao, X.; Huidobro, P. A.; Giannini, V. Topological Plasmonic Chain with Retardation and Radiative Effects. _ACS Photonics_ 2018, _5_ , 2271–2279
  * Humphrey et al. 2016 Humphrey, A. D.; Meinzer, N.; Starkey, T. A.; Barnes, W. L. Surface Lattice Resonances in Plasmonic Arrays of Asymmetric Disc Dimers. _ACS Photonics_ 2016, _3_ , 634–639
  * Humphrey and Barnes 2016 Humphrey, A.; Barnes, W. Plasmonic Surface Lattice Resonances in Arrays of Metallic Nanoparticle Dimers. _J. Opt._ 2016, _18_ , 035005
  * Guo et al. 2017 Guo, R.; Hakala, T. K.; Törmä, P. Geometry Dependence of Surface Lattice Resonances in Plasmonic Nanoparticle Arrays. _Phys. Rev. B: Condens. Matter Mater. Phys._ 2017, _95_ , 155423
  * Hong and Reinhard 2014 Hong, Y.; Reinhard, B. M. Collective photonic-plasmonic resonances in noble metal - dielectric nanoparticle hybrid arrays. _Optical Materials Express_ 2014, _4_ , 2409–2422
  * Cuartero-González et al. 2020 Cuartero-González, A.; Sanders, S.; Zundel, L.; Fernández-Domínguez, A. I.; Manjavacas, A. Super- and Subradiant Lattice Resonances in Bipartite Nanoparticle Arrays. _ACS Nano_ 2020, _14_ , 11876–11887
  * Esposito et al. 2019 Esposito, M.; Todisco, F.; Bakhti, S.; Passaseo, A.; Tarantini, I.; Cuscunà, M.; Destouches, N.; Tasco, V. Symmetry Breaking in Oligomer Surface Plasmon Lattice Resonances. _Nano Letters_ 2019, _19_ , 1922–1930
  * Baur et al. 2018 Baur, S.; Sanders, S.; Manjavacas, A. Hybridization of Lattice Resonances. _ACS Nano_ 2018, _12_ , 1618–1629
  * Todisco et al. 2016 Todisco, F.; Esposito, M.; Panaro, S.; De Giorgi, M.; Dominici, L.; Ballarini, D.; Fernández-Domínguez, A. I.; Tasco, V.; Cuscunà, M.; Passaseo, A.; Ciracì, C.; Gigli, G.; Sanvitto, D. Toward Cavity Quantum Electrodynamics with Hybrid Photon Gap-Plasmon States. _ACS Nano_ 2016, _10_ , 11360–11368
  * Mahi et al. 2017 Mahi, N.; Lévêque, G.; Saison, O.; Marae-Djouda, J.; Caputo, R.; Gontier, A.; Maurer, T.; Adam, P.-M.; Bouhafs, B.; Akjouj, A. In Depth Investigation of Lattice Plasmon Modes in Substrate-Supported Gratings of Metal Monomers and Dimers. _The Journal of Physical Chemistry C_ 2017, _121_ , 2388–2401
  * Guo et al. 2017 Guo, K.; Du, M.; Osorio, C. I.; Koenderink, A. F. Broadband light scattering and photoluminescence enhancement from plasmonic Vogel’s golden spirals. _Laser & Photonics Reviews_ 2017, _11_ , 1600235
  * Marinica et al. 2008 Marinica, D. C.; Borisov, A. G.; Shabanov, S. V. Bound States in the Continuum in Photonics. _Physical Review Letters_ 2008, _100_ , 183902
  * Plotnik et al. 2011 Plotnik, Y.; Peleg, O.; Dreisow, F.; Heinrich, M.; Nolte, S.; Szameit, A.; Segev, M. Experimental Observation of Optical Bound States in the Continuum. _Physical Review Letters_ 2011, _107_ , 183901
  * Rybin et al. 2017 Rybin, M. V.; Koshelev, K. L.; Sadrieva, Z. F.; Samusev, K. B.; Bogdanov, A. A.; Limonov, M. F.; Kivshar, Y. S. High-$Q$ Supercavity Modes in Subwavelength Dielectric Resonators. _Physical Review Letters_ 2017, _119_ , 243901
  * Hsu et al. 2016 Hsu, C. W.; Zhen, B.; Stone, A. D.; Joannopoulos, J. D.; Soljačić, M. Bound states in the continuum. _Nature Reviews Materials_ 2016, _1_ , 16048
  * Han et al. 2020 Han, S.; Rybin, M. V.; Pitchappa, P.; Srivastava, Y. K.; Kivshar, Y. S.; Singh, R. Guided-Mode Resonances in All-Dielectric Terahertz Metasurfaces. _Advanced Optical Materials_ 2020, _8_ , 1900959
  * Kirill et al. 2019 Kirill, K.; Gael, F.; Andrey, B.; Yuri, K.; Andrea, F. Nonradiating photonics with resonant dielectric nanostructures. _Nanophotonics_ 2019, _8_ , 725–745
  * Fan et al. 2019 Fan, K.; Shadrivov, I. V.; Padilla, W. J. Dynamic bound states in the continuum. _Optica_ 2019, _6_ , 169–173
  * Kim et al. 2019 Kim, S.; Kim, K.-H.; Cahoon, J. Optical Bound States in the Continuum with Nanowire Geometric Superlattices. _Physical Review Letters_ 2019, _122_ , 187402
  * Koshelev et al. 2018 Koshelev, K.; Lepeshov, S.; Liu, M.; Bogdanov, A.; Kivshar, Y. Asymmetric Metasurfaces with High-$Q$ Resonances Governed by Bound States in the Continuum. _Physical Review Letters_ 2018, _121_ , 193903
  * Hsu et al. 2013 Hsu, C. W.; Zhen, B.; Lee, J.; Chua, S.-L.; Johnson, S. G.; Joannopoulos, J. D.; Soljačić, M. Observation of trapped light within the radiation continuum. _Nature_ 2013, _499_ , 188–191
  * Liu and Choi 2018 Liu, M.; Choi, D.-Y. Extreme Huygens’ Metasurfaces Based on Quasi-Bound States in the Continuum. _Nano Letters_ 2018, _18_ , 8062–8069
  * Yoon et al. 2015 Yoon, J. W.; Song, S. H.; Magnusson, R. Critical field enhancement of asymptotic optical bound states in the continuum. _Scientific Reports_ 2015, _5_ , 18301
  * Karl et al. 2019 Karl, N.; Vabishchevich, P. P.; Liu, S.; Sinclair, M. B.; Keeler, G. A.; Peake, G. M.; Brener, I. All-optical tuning of symmetry protected quasi bound states in the continuum. _Applied Physics Letters_ 2019, _115_ , 141103
  * Abujetas et al. 2019 Abujetas, D. R.; van Hoof, N.; ter Huurne, S.; Gómez Rivas, J.; Sánchez-Gil, J. A. Spectral and temporal evidence of robust photonic bound states in the continuum on terahertz metasurfaces. _Optica_ 2019, _6_ , 996
  * Abujetas et al. 2019 Abujetas, D. R.; Barreda, Á.; Moreno, F.; Sáenz, J. J.; Litman, A.; Geffrin, J.-M.; Sánchez-Gil, J. A. Brewster quasi bound states in the continuum in all-dielectric metasurfaces from single magnetic-dipole resonance meta-atoms. _Sci. Rep._ 2019, _9_ , 16048
  * Koshelev et al. 2018 Koshelev, K. L.; Sychev, S. K.; Sadrieva, Z. F.; Bogdanov, A. A.; Iorsh, I. V. Strong coupling between excitons in transition metal dichalcogenides and optical bound states in the continuum. _Physical Review B_ 2018, _98_ , 161113
  * Nguyen et al. 2019 Nguyen, T. G.; Ren, G.; Schoenhardt, S.; Knoerzer, M.; Boes, A.; Mitchell, A. Ridge Resonance in Silicon Photonics Harnessing Bound States in the Continuum. _Laser & Photonics Reviews_ 2019, _13_ , 1900035
  * Han et al. 2019 Han, S.; Cong, L.; Srivastava, Y. K.; Qiang, B.; Rybin, M. V.; Kumar, A.; Jain, R.; Lim, W. X.; Achanta, V. G.; Prabhu, S. S.; Wang, Q. J.; Kivshar, Y. S.; Singh, R. All-Dielectric Active Terahertz Photonics Driven by Bound States in the Continuum. _Advanced Materials_ 2019, _31_ , 1901921
  * Sadrieva et al. 2017 Sadrieva, Z. F.; Sinev, I. S.; Koshelev, K. L.; Samusev, A.; Iorsh, I. V.; Takayama, O.; Malureanu, R.; Bogdanov, A. A.; Lavrinenko, A. V. Transition from Optical Bound States in the Continuum to Leaky Resonances: Role of Substrate and Roughness. _ACS Photonics_ 2017, _4_ , 723–727
  * Algorri et al. 2021 Algorri, J. F.; Dell’Olio, F.; Roldán-Varona, P.; Rodríguez-Cobo, L.; López-Higuera, J. M.; Sánchez-Pena, J. M.; Zografopoulos, D. C. Strongly resonant silicon slot metasurfaces with symmetry-protected bound states in the continuum. _Optics Express_ 2021, _29_ , 10374–10385
  * Kodigala et al. 2017 Kodigala, A.; Lepetit, T.; Gu, Q.; Bahari, B.; Fainman, Y.; Kanté, B. Lasing action from photonic bound states in continuum. _Nature_ 2017, _541_ , 196–199
  * Ha et al. 2018 Ha, S. T.; Fu, Y. H.; Emani, N. K.; Pan, Z.; Bakker, R. M.; Paniagua-Domínguez, R.; Kuznetsov, A. I. Directional lasing in resonant semiconductor nanoantenna arrays. _Nat. Nanotechnol._ 2018, _13_ , 1042–1047
  * Yesilkoy et al. 2019 Yesilkoy, F.; Arvelo, E. R.; Jahani, Y.; Liu, M.; Tittl, A.; Cevher, V.; Kivshar, Y.; Altug, H. Ultrasensitive hyperspectral imaging and biodetection enabled by dielectric metasurfaces. _Nat. Photonics_ 2019, _13_ , 390–396
  * Wu et al. 2020 Wu, M.; Ha, S. T.; Shendre, S.; Durmusoglu, E. G.; Koh, W. K.; Abujetas, D. R.; Sánchez-Gil, J. A.; Paniagua-Domínguez, R.; Demir, H. V.; Kuznetsov, A. I. Room-Temperature Lasing in Colloidal Nanoplatelets via Mie-Resonant Bound States in the Continuum. _Nano Lett._ 2020, _20_ , 6005–6011
  * Abujetas et al. 2019 Abujetas, D. R.; van Hoof, N.; ter Huurne, S.; Gómez Rivas, J.; Sánchez-Gil, J. A. Spectral and temporal evidence of robust photonic bound states in the continuum on terahertz metasurfaces. _Optica_ 2019, _6_ , 996–1001
  * Abujetas et al. 2018 Abujetas, D. R.; Sánchez-Gil, J. A.; Sáenz, J. J. Generalized Brewster effect in high-refractive-index nanorod-based metasurfaces. _Opt. Express_ 2018, _26_ , 31523
  * Abujetas et al. 2020 Abujetas, D. R.; Olmos-Trigo, J.; Sáenz, J. J.; Sánchez-Gil, J. A. Coupled electric and magnetic dipole formulation for planar arrays of particles: Resonances and bound states in the continuum for all-dielectric metasurfaces. _Phys. Rev. B_ 2020, _102_ , 125411
  * Abujetas and Sánchez-Gil 2021 Abujetas, D. R.; Sánchez-Gil, J. A. Near-Field Excitation of Bound States in the Continuum in All-Dielectric Metasurfaces through a Coupled Electric/Magnetic Dipole Model. _Nanomaterials_ 2021, _11_ , 998
  * Murai et al. 2020 Murai, S.; Abujetas, D. R.; Castellanos, G. W.; Sánchez-Gil, J. A.; Zhang, F.; Rivas, J. G. Bound States in the Continuum in the Visible Emerging from out-of-Plane Magnetic Dipoles. _ACS Photonics_ 2020, _7_ , 2204–2210
  * Murai et al. 2020 Murai, S.; Gabriel W. Castellanos, G. W.; Raziman, T.; Curto, A. G.; Gómez Rivas, J. Enhanced light emission by magnetic and electric resonances in dielectric metasurfaces. _Adv. Opt. Mater._ 2020, _8_ , 1902024
  * Castellanos et al. 2020 Castellanos, G. W.; Murai, S.; Raziman, T. V.; Wang, S.; Ramezani, M.; Curto, A. G.; Gómez Rivas, J. Exciton-Polaritons with Magnetic and Electric Character in All-Dielectric Metasurfaces. _ACS Photonics_ 2020, _7_ , 1226–1234
  * Wang et al. 2020 Wang, S.; Raziman, T. V.; Murai, S.; Castellanos, G. W.; Bai, P.; Berghuis, A. M.; Godiksen, R. H.; Curto, A. G.; Gómez Rivas, J. Collective Mie Exciton-Polaritons in an Atomically Thin Semiconductor. _The Journal of Physical Chemistry C_ 2020, _124_ , 19196–19203
  * Abujetas et al. 2018 Abujetas, D. R.; Paniagua-Domínguez, R.; Sánchez-Gil, J. A. Impedance-Matched, Double-Zero Optical Metamaterials Based on Weakly Resonant Metal Oxide Nanowires. _Photonics_ 2018, _5_ , 7
  * Murai et al. 2020 Murai, S.; Abujetas, D. R.; Castellanos, G. W.; Sánchez-Gil, J. A.; Zhang, F.; Rivas, J. G. Bound States in the Continuum in the Visible Emerging from out-of-Plane Magnetic Dipoles. _ACS Photonics_ 2020, _7_ , 2204–2210
  * Reid and Johnson 2013 Reid, M. H.; Johnson, S. G. Efficient Computation of Power, Force, and Torque in BEM Scattering Calculations. _arXiv:1307.2966_ 2013,
  * 61 Reid, M. T. H. SCUFF-EM. _http://homerreid.github.io/scuff-em-documentation/_
  * Aspnes and Studna 1983 Aspnes, D. E.; Studna, A. A. Dielectric functions and optical parameters of Si, Ge, GaP, GaAs, GaSb, InP, InAs, and InSb from 1.5 to 6.0 eV. _Physical Review B_ 1983, _27_ , 985–1009
  * Schinke et al. 2015 Schinke, C.; Peest, P. C.; Schmidt, J.; Brendel, R.; Bothe, K.; Vogt, M. R.; Kröger, I.; Winter, S.; Schirmacher, A.; Lim, S.; Nguyen, H. T.; MacDonald, D. Uncertainty analysis for the coefficient of band-to-band absorption of crystalline silicon. _AIP Advances_ 2015, _5_ , 067168
  * Green 2008 Green, M. A. Self-consistent optical parameters of intrinsic silicon at 300K including temperature coefficients. _Solar Energy Materials and Solar Cells_ 2008, _92_ , 1305–1310

Figure 1: Structure and extinction of the nanoparticle arrays. SEM images of
(a) the symmetric array, where Si nanodisks (200 nm height, 224 $\pm$ 7.2 nm
diameter) are arranged in the $x$\- and $y$-directions with the period of 425
nm, (b) the position-detuned array, where the pair of Si nanodisks (200 nm
height, 236 $\pm$ 6.6 nm diameter, center-to-center distance 340 nm) are
arranged in the $x$\- and $y$-directions with the period of 850 and 425 nm,
respectively, and (c) the size-detuned array, where the pair of Si nanodisks
(200 nm height, 194 $\pm$ 3.6 and 256 $\pm$ 6.6 nm diameters, center-to-center
distance 425 nm) are arranged in the $x$\- and $y$-directions with the period
of 850 and 425 nm, respectively, and (d) doubly-detuned array, where the pair
of Si nanodisks (200 nm height, 210 $\pm$ 7.0 and 276 $\pm$ 2.6 nm diameters,
center-to-center distance 340 nm) are arranged in the $x$\- and $y$-directions
with the period of 850 and 425 nm, respectively. The unit cell is denoted as
white dashed boxes. The sketches of the disk arrangement along the
$x$-direction are on the top. Extinction spectra (1-zeroth-order
transmittance) for TE light as a function of angle of incidence $\theta_{\rm
in}$ for the (e) symmetric array, (f) position-detuned array, (g) size-detuned
array, and (h) doubly-detuned array. The in-plane diffraction condition is
denoted as dotted lines. Extinction spectra (TE) around the normal incidence
($\theta_{\rm in}$ between 0 and 5.14 ∘ in steps of 0.514 ∘ step) and the mode
assignment: (i) the symmetric array, (j) the size-detuned array, (k) the
position-detuned array, and (l) the doubly-detuned array. The in-plane
diffraction condition is denoted by vertical ticks.

Figure 2: Schematic of mode evolution under lattice modification and
correlation with extinction spectra (TE) at normal incidence ($\theta_{\rm
in}$ = 0∘). (a) Evolution of the out-of-plane magnetic dipoles and (b) of the
in-plane electric dipoles eigenmodes from the guided modes in symmetric arrays
to bright and BIC modes in other lattice symmetries. GM: guided mode, B:
bright mode. (c-f) Extinction spectra (TE) at normal incidence taken from Fig.
1 (i-l). (c) Symmetric, (d) position-detuned, (e) size-detuned and (f) double-
detuned arrays. The insets show the eigenmodes supported by each lattice
symmetry and the vertical dashed lines indicate the in-plane diffraction
condition.

Figure 3: Extinction spectra ($1-T$, $T$ being the specular transmittance) of
the nanoparticle arrays calculated using the CEMD method for: (a) a symmetric
array, where Si nanodisks (200 nm height, 225 nm diameter) are arranged in a
square lattice with a period of 425 nm; (b) a position-detuned array, where Si
nanodisk dimers (200 nm height, 225 nm diameter, center-to-center distance 340
nm) are arranged in a rectangular lattice with periods of 850 and 425 nm; (c)
a size-detuned array, where Si nanodisk dimers (200 nm height, 200 and 250 nm
diameters, center-to-center distance 425 nm) are arranged in a rectangular
lattice with periods of 850 and 425 nm; and (d) a doubly-detuned array, where
Si nanodisk dimers (200 nm height, 210 and 270 nm diameters, center-to-center
distance 340 nm) are arranged in a rectangular lattice with periods of 850 and
425 nm, respectively. Upper panels show the extinction spectra (TE) as a
function of $\theta_{\rm in}$. (e-h) The bottom panels are the spectra shown
in (a-d) between $\theta_{\rm in}$= 0∘ and 5.14∘ in steps of 0.514∘. The color
of the curves denotes the projection over the collective modes shown in the
insets in Figs. 2(b-d): (red) low-energy mode with anti-parallel pairs of EDy
and parallel MDz; and (blue) high-energy mode with parallel EDy and anti-
parallel MDz. The non-diffraction region is delimited by the vertical ticks.

Figure 4: $Q$-factors of quasi-BIC mode extracted from the extinction as a
function of $\theta_{\rm in}$ for TE polarized light. (a) the position-detuned
array, and (b) the size-detuned array. The dotted line indicate the upper
limit of the $Q$-factor imposed by the resolution of the spectrometer. The
$Q$-factors extracted from the CEMD calculation in Figure 3 are shown as solid
lines. The upper and lower insets are the schematic representation of the
array and the optical extinction around normal incidence.

Figure 5: Optical properties of the size-detuned array resonating in the
visible. The array consists of the pair of Si nanodisks (200 nm height, 86 and
107 nm diameters, center-to-center distance 200 nm) are arranged in the $x$\-
and $y$-directions with the period of 400 and 200 nm, respectively. (a) and
(b): Extinction spectra as a function of $\theta_{\rm in}$ for (a) TE and (b)
TM polarized light. (c) $Q$-factors of quasi-BIC mode extracted from the
extinction as a function of $\theta_{\rm in}$ for TE-polarized light. The top
inset is the extinction spectra around normal incidence ($\theta_{\rm in}$
between 0 to 5.14∘ every 0.68 ∘ step). The spectra are shifted vertically for
the sake of clarity. The bottom inset is a top-view SEM image of the array.
The dotted rectangle denotes the unit cell (400 nm $\times$ 200 nm) of the
array.