rosshort = ROS, long = Robot Operating System slamshort = SLAM, long =
simultaneous localisation and mapping lidarshort = LiDAR, long = light
detection and ranging imushort = IMU, long = inertial measurement unit
gprshort = GPR, long = Gaussian process regression gpshort = GP, long =
Gaussian process sarshort = SAR, long = search and rescue tcpshort = TCP, long
= tool center point nbcshort = NBC, long = nuclear biological chemical
enrichshort = EnRicH, long = European Robotics Hackathon ugvshort = UGV, long
= uncrewed ground vehicle uavshort = UAV, long = uncrewed areal vehicle
dofshort = DOF, long = degree of freedom

# Semi-Autonomous Mobile Search and Rescue Robot
for Radiation Disaster Scenarios

Simon Schwaiger∗1, Lucas Muster∗1,2, Georg Novotny1, Michael Schebek1,
Wilfried Wöber1, Stefan Thalhammer1 and Christoph Böhm1 *Equal
contributionThis work was supported by the Austrian Research Promotion Agency
(FFG) through the research project UGV-ABC-Probe (FFG project Call 2020) and
the Austrian Armed Forces.1Simon Schwaiger, Lucas Muster, Georg Novotny,
Michael Schebek, Wilfried Wöber, Stefan Thalhammer and Christoph Böhm are with
the Faculty of Industrial Engineering, University of Applied Sciences
Technikum Wien, 1200 Vienna, Austria<EMAIL_ADDRESS><EMAIL_ADDRESS>Muster is with the University of Natural
Resources and Life Sciences, Department of Biotechnology, Institute for
Computational Biology, Muthgasse 18, 1190 Vienna, Austria

###### Abstract

This paper describes a novel semi-autonomous mobile robot system designed to
assist search and rescue (SAR) first responders in disaster scenarios. While
robots offer significant potential in SAR missions, current solutions are
limited in their ability to handle a diverse range of tasks. This gap is
addressed by presenting a system capable of (1) autonomous navigation and
mapping, allowing the robot to autonomously explore and map areas affected by
catastrophic events, (2) radiation mapping, enabling the system to triangulate
a radiation map from discrete radiation measurements to aid in identifying
hazardous areas, (3) semi-autonomous substance sampling, allowing the robot to
collect samples of suspicious substances and analyze them onboard with
immediate classification, and (4) valve manipulation, enabling teleoperated
closing of valves that control hazardous material flow. This semi-autonomous
approach balances human control over critical tasks like substance sampling
with efficient robot navigation in low-risk areas. The system is evaluated
during three trials that simulate possible disaster scenarios, two of which
have been recorded during the enrich (enrich). Furthermore, we provide
recorded sensor data as well as the implemented software system as
supplemental material through a GitHub repository: https://github.com/TW-
Robotics/search-and-rescue-robot-IROS2024.

##  I INTRODUCTION

In the context of sar, robots are deployed to aid emergency teams in
responding faster and more efficiently in case of a disaster scenario [1, 2,
3, 4]. Robots are used in dangerous environments to collect important
information about the area of operation [5, 6, 7] to keep human first
responders away from a potentially hazardous environment and minimise their
exposure to harmful substances. These robotic systems need to handle diverse
tasks, such as exploration and mapping of disaster zones [8], localizing
hazardous material [9], and solving manipulation tasks [10]. However, sar
robots with a varying skill-set are not a standard piece of equipment of SAR
response teams, therefore, their potential is not fully realized [11].

Recent advances have been made in the implementation of mapping [12],
hazardous material localization [13, 14], and manipulation tasks [3]. Those
systems address and solve individual tasks. However, for efficient deployment,
their capabilities need to be combined on a single platform.

Figure 1: Robot capabilities. The robot system solves four tasks on one
platform to assist first responders in disaster scenarios: Autonomous
navigation and mapping, radiation mapping, semi-autonomous substance sampling
and valve manipulation.

Our presented system performs a diverse set of tasks in areas affected by
catastrophes based on requirements set by real world problems encountered by
the Austrian Armed Forces. Our robot, as depicted in Fig. 1, fulfills these
requirements by being capable of

  * •

autonomously navigating and exploring the vicinity and creating a map of the
robot’s surrounding geometry,

  * •

measuring radiation and triangulating a radiation map,

  * •

taking samples of substances on the floor in a semi-autonomous fashion with
subsequent on-board analysis, and

  * •

closing valves that control the flow of hazardous material.

This paper advances the current state of the art in sar robotics by designing
and implementing a semi-autonomous mobile robot system. Our experiments show,
that the robot provides these capabilities in a unified platform. Therefore,
our robot potentially improves support for rescue teams during disasters.

The remainder of this paper is structured as follows. Sec. II presents recent
advances in autonomous sar systems. Sec. III provides a methodological
overview of the sar robot and the fulfilled tasks, while Sec. IV discusses
robot hardware and software. Sec. V describes the robot’s evaluation in three
scenarios and Sec. VI summarizes our research and provides an outlook on
future work.

##  II RELATED WORK

The increased interest in mobile robots aiding in sar missions motivates
researchers and companies to create solutions for the tasks encountered in
these scenarios. These tasks range from teleoperated exploration or inspection
to semi-autonomous mission execution. Challenges stemming from given real-
world problems allow the evaluation of mobile robotic systems in safe and
controlled environments. The result of all those efforts are fast to deploy
and easy to use systems for real-world disaster sites.

###  II-A Competitions and Challenges

Past and current competitions and challenges emulate a single task or a whole
series of tasks (mission) found in real-world problems. These tasks can vary
in scale, from enclosed rooms to whole sites, and their complexity, from
exploring and inspection to fully autonomous interactions with the
environment.

Conferences often hold competitions on a smaller scale with single task
requirements, e.g., ICRA 2023 BARN
Challenge111https://www.icra2023.org/competitions. Bigger challenges that
contain whole missions for mobile robots such as the DARPA Subterranean (SubT)
Challenge222https://www.darpa.mil/program/darpa-subterranean-challenge, Total
Autonomous Robot for Gas and Oil Sites (ARGOS), or The European Land Robot
Trial (ELROB)333https://www.elrob.org/ often originate from the need of
practical sar solutions. Hence, they provide a valuable proving ground for sar
robots. Our robotic system was present at the enrich
(enrich)444https://enrich.european-robotics.eu/ in 2023 and data was collected
to evaluate the effectiveness of the system. During the challenge, a hazardous
material incident response operation was staged.

###  II-B Mobile Robotic Systems for Search and Rescue

One of the first uses of mobile robots in sar missions at the World Trade
Center (WTC) disaster site 2001 is documented in [1]. Tethered tracked robot
vehicles with cameras that are teleoperated helped rescuers on site. The team
conducted multiple runs with the following tasks: (1) searching for victims,
(2) searching for faster excavation paths, (3) structural inspection, and (4)
detection of hazardous materials. These robots were able to enter spaces too
small or dangerous for humans or rescue dogs. Another example of mobile
robotic systems aiding humans was showcased 2011 in Japan during the Fukushima
nuclear accident. Researchers deployed the teleoperated mobile robotic
platform TALON equipped with a radiation sensor [2]. The heavy radiation
prohibited humans from entering the area, therefore, a robot was used to
localize radiation sources, using lidar depth information and radiation
measurements.

Over the last two decades, researchers came up with many solutions to mobile
robots in sar missions. One example of such a robotic system is the work by
the authors of [10] which proposed a remotely operated mobile robot for
radiation mapping and semi-autonomous sampling of radioactive sources. Many
robots [15, 9, 16, 14] solely focusing on radiation mapping were presented by
the research community. These systems use discrete maps of the contaminated
area and navigate through user predefined way-points or coverage areas. The
robots measure radiation levels during the operation and generate a radiation
map. This map can be used by the users themselves or navigation packages use
them to avoid sources. [17] is an example of autonomous exploration (coverage
planning) and environmental measuring based on behavior trees for task
execution. This robot participated at the EnRicH in 2019. Almost all solutions
rely on the user selecting points of interest, while the robotic system
presented in this work uses exploration to autonomously determine interesting
spots to complete the radiation map.

An alternative to ground vehicles in the context of sar robotics are uav. [18,
19] automate radiation mapping using flying autonomous systems, allowing for
coverage of bigger areas. However, uavs are limited by their maximum take-off
weight and flight time.

Most solutions in radiation mapping lack the capability of manipulation, a
needed skill in sar missions. This is emphasized by the application shown in
research [20, 21, 3]. All robots are able to interact with their environment,
e.g., removing potential obstacles with a robotic manipulator or evacuate
victims from the disaster site by carrying them away. These robots are
teleoperated, therefore, control and manipulation happens based on direct user
input or a list of way-points.

Our approach incorporates semi-autonomous manipulation tasks, e.g., sampling
and analysis of potentially hazardous materials, to aid disaster response.

In the current state of the art, research gaps are identified with regard to
robots’ abilities to perform diverse tasks consisting of reconnaissance as
well as manipulation. Furthermore, systems incorporating manipulation act in a
teleoperated manner and lack autonomy. We address these gaps by implementing
reconnaissance in the form of geometry and radiation mapping as well as
autonomous exploration. Additionally, our system is capable of performing two
manipulation tasks, i.e., substance sampling and valve manipulation.

Figure 2: Functional system overview. The proposed system aims to provide 2D
and 3D mapping, radiation mapping, autonomous navigation and manipulation
capabilities (blue rounded border) for sar missions. The figure depicts the
data flow between functional units, sensors and the multi-purpose end-
effector.

##  III METHOD

The presented sar robot’s capabilities are depicted in the form of a high-
level overview in Fig. 2. The used tracked ground vehicle [22] has a five dof
arm mounted onto the chassis for manipulation tasks. The mobile platform is
equipped with sensors and an end-effector that enable it to complete required
tasks in a semi-autonomous fashion. Two lidar sensors and an imu allow for 2D
and 3D slam using lidar-odometry and mapping in 3D and a graph-based slam in
2D. The generated 2D map serves as the basis for autonomous navigation. Robot
movement is planned using a finite state machine with recovery behaviors. A
frontier-based exploration algorithm proposes goal poses, that allow the robot
to extend its discovered map. The generated 2D map is fused with measurements
from a radiation sensor in order to perform radiation mapping. The setup is
complemented by a multi-purpose end-effector that enables sampling of
substances encountered during a mission using a Raman spectrometer.
Furthermore, the end-effector is capable of manipulating valves that control
the flow of hazardous material.

###  III-A 2D and 3D Mapping

When implementing the system, we adhered to the conventions set by the ros
[23] and its reference robots [24] to enable compatibility with existing
software implementations. The mapping setup consists of a separate 3D as well
as a more efficient 2D slam system running on a industrial PC. This approach
was chosen to enable autonomous navigation without relying on the
computationally expensive 3D mapper. The utilised 2D mapper was chosen based
on its wide adoption as the by default supported mapper in ros 2 [25].

The 3D mapper is used for visualization with the intent of providing the
robot’s non-technical users with a more realistic representation of the
robot’s environment compared to a 2D map. 3D mapping is performed on data from
the 3D lidar and imu [26]. The combination of lidar and imu is more robust
than lidar-only mapping, especially when handling aggressive motion and sudden
robot rotation [27]. Since a separate 2D slam is the basis of the navigation
system, the robot behavior is unaffected by disabling online 3D mapping in
order to save performance.

###  III-B Autonomous Navigation Stack

Autonomous navigation is achieved by planning a trajectory from the current
pose to a goal pose based on the map and robot pose estimated by the 2D
mapper. Layered costmaps [28] have been established as the standard approach
to mobile robot trajectory planning in ROS [29]. However, to deal with track
slippage and make planning more robust, costmap-based trajectory planning is
performed using hierarchical finite state machines that are complemented by
custom recovery behaviors that are activated if the robot does not
successfully navigate to a desired pose. The implemented recovery behaviors
allow the robot to autonomously (1) move away in case of a collision and (2)
clear the local costmap if no path can be planned in the robot’s direct
vicinity.

The goal pose is either provided by a human operator or the frontier-based
exploration module, which proposes goal poses for the robot to navigate to in
order to expand the mapped area. Therefore, autonomous and teleoperated modes
are intelligently blended to achieve semi-autonomous behavior. Control signals
given by the human operator are favored over control signals given by the
exploration or trajectory planning modules. This approach allows the robot’s
users to interrupt autonomous behavior by providing direct control signals or
setting goal poses using the graphical user interface.

The robot’s standard behavior is also interrupted if it loses connection to
the operator PC. If the operator can not be communicated with, for example,
due to a wireless communication breakdown, a routine is executed that preempts
exploration and navigates the robot back towards its starting pose.

###  III-C Radiation Mapping

The mapping of nuclear radiation is based on count data gathered by a Geiger
counter. While the mapping of such radiation was previously done by estimating
the posterior predictive distribution described by a Poisson distribution
[30], autonomous mapping using robots is a challenging task due to the
presence of both, process and sensor noise. A radiation mapping pipeline
optimized for autonomous mobile robots is proposed in [31]. The authors apply
non-parametric gpr [32] using adapted kernels to map the radiation based on
Geiger counts. In this context, a gp models the spatial distribution of
radiation and interpolates between sampling points.

Our proposed approach for mapping is inspired by [10], however, we extend
their work by applying gp to mapping, as seen in [31]. The developed radiation
mapping pipeline combines the mobile robot’s current pose (see Sec. IV-B) with
count data gathered by the digital Geiger counter. At each time step $t$, the
gathered radiation count data $c_{t}$ is synchronized with the robot’s pose
$x_{t}$ estimated by the 2D mapper [25]. The training data $\mathcal{D}$
describes the gathered poses and count data.

In a similar manner to [31], the map is estimated by defining the posterior
predictive distribution
$p(c_{new}|x_{new},x_{1:t},c_{1:t})=p(c_{new}|x_{new},\mathcal{D})$ relying on
gpr with radial base function kernels [32]. Hence, for each location $(i,j)$
in the map $m$, we can estimate the radiation using Eq. (1), where the
functions $GP_{\mu}(.)$ and $GP_{\Sigma}(.)$ describe the prediction functions
of the gp.

$\forall{i,j}\in m,p(c_{i,j}|x_{i,j},\mathcal{D})\approx\\\
\mathcal{N}\left(GP_{\mu}(x_{i,j},\mathcal{D}),GP_{\Sigma}(x_{i,j},\mathcal{D})\right)$
(1)

The mapping procedure described in Eq. (1) is applied in an offline manner,
i.g., after the robot has explored unknown environments.

###  III-D Manipulation System

The manipulation system consists of an industrial robotic arm with five dof
that is complemented by a custom three-finger gripper as end-effector. The
end-effector is designed to perform two diverse tasks, namely closing valves
as well as grasping the probe used for sampling. The gripper’s tcp is located
along the robot arm’s fifth axis in order to turn valves and correctly orient
samples into the analysis tray (see Fig. 3).

Figure 3: Custom gripper for valve manipulation. The left part of the figure
shows a real image of the custom gripper and the right part shows the modelled
end-effector with highlighted individual parts.

The end-effector grasps a probe used to absorb potentially harmful substances
for on-board analysis. This probe consists of an interface to the gripper in
cylindrical form and a glass wool swab. The construction allows the swab to be
correctly rotated towards the sensor for substance analysis.

Up to two probes can be stored in a storage tray on the sar robot. One
additional probe slot is present at the back of the robot, where samples of
substances are examined by a Raman spectrometer.

##  IV IMPLEMENTATION DETAIL

As presented in [33], our implementation utilizes a synergistic integration of
Docker containers555https://www.docker.com/ to encapsulate software packages,
systemd service management666https://systemd.io/ to start, stop and restart
groups of containers, and the Cockpit monitoring tool777https://cockpit-
project.org/ for monitoring software packages running within each systemd
service.

Docker containers encapsulate and isolate our robot’s individual software
components. This approach allows each software package’s diverse dependencies
to be managed individually. It enables us to set up software in a portable and
hardware-agnostic manner, while selectively allowing communication between
software components when needed.

Systemd, a service management system, automates and manages the lifecycle of
Docker containers in groups that represent the main stages (sensing,
perception, mapping, exploration, navigation) of the sar robot’s data
processing pipeline. Systemd ensures that these groups are reliably started,
stopped, and restarted as needed based on whether or not the robot’s core
drivers are running. The setup utilizes systemd template unit files that start
Docker compose files, allowing the concurrent management of multiple container
instances.

The Cockpit monitoring tool provides a comprehensive view of system status,
including the health and activity of both Docker containers and systemd
services, facilitating efficient monitoring by human operators. This approach
allows for immediate identification and troubleshooting of potential problems,
aiding in the sar robot’s continuous availability. Furthermore, Cockpit allows
Docker containers and systemd services to be started, stopped and restarted
manually using the graphical user interface. In combination with having
software components encapsulated, the user interface allows for parts of the
software stack to be restarted during a mission without shutting down the
whole system. This functionality is not natively supported by ros when using
the robot control pipeline using a launchfile.

The remainder of this section details the specific components of the robotic
system, including utilized hardware, the software implementing autonomous
navigation and the deployed endeffector.

###  IV-A Hardware

The sar robot’s built-in industrial PC features an Intel i7 7700T CPU, Nvidia
GTX 1050 Ti GPU and 32 GB of system memory. An Nvidia Jetson Xavier NX 16GB
single board computer interfaces with the front cameras. A Realsense D455 is
mounted to the robot’s arm for teleoperated manipulation, while a Realsense
L515 serves as the front lidar used to fill the main 3D lidar’s scan. The main
3D lidar is a Velodyne VLP-16, while the utilized imu is a LORD
3DMGX5-GNSS/INS. Radiation is measured using a Flir Identifinder R300 sensor,
while the sample analysis sensor has the model name Thermo Scientific
FirstDefender RMX Handheld Chemical Identification.

###  IV-B Navigation and Mapping Software

The system leverages ros 1 (version Noetic) on the industrial PC, while the
single board PC runs on ros 2 (version Humble). Slam Toolbox [25] serves as
the system’s 2D mapper, while trajectory planning is performed by Move Base
Flex [29]. LIO-SAM [26] performs 3D mapping based on 3D lidar and imu
measurements. Explore Lite888https://github.com/hrnr/m-explore generates goal
poses, where the map depicting the discovered area can be expanded.

Autonomous and teleoperated operation modes are automatically switched between
using priority-based multiplexers. When switching from a signal source with a
higher priority to a lower priority (e.g., switching from a teleoperation-
based command velocity to a velocity generated by the trajectory planner)
switching occurs with a predefined cooldown to prevent the robot base from
oscillating between conflicting control signals.

Two multiplexers are deployed across the pipeline. The first multiplexer
switches between command velocities generated by the operator (teleoperation)
and the autonomous navigation stack. The second multiplexer switches the goal
pose between the following sources in descending priority: (1) Go home
routine, (2) planning preemption, (3) user interface and (4) explorer. Go home
routine depicts the stored home pose that the robot navigates to in case the
connection to the operator is interrupted, while planning preemption depicts a
node that preempts trajectory planning if the operator performs teleoperation.
In this context, user interface refers to goal poses that the operator can
send by clicking on the 2D map in the user interface and explorer refers to
goalposes proposed by the map exploration module.

###  IV-C End-effector

The sar robot is equipped with a custom end-effector that is capable of both,
closing valves as well as grasping the probe used for sampling. Manipulated
objects are automatically centered through the gripper’s kinematics and can
therefore be rotated using the robot arm’s fifth axis that is aligned with the
end-effector’s tcp. After accounting for the reduction in the planetary
gearbox, the resulting gripping force of 80 Newtons on the tcp is generated by
a Robotis XM430-W350 servo motor. Assuming centered placement of the tcp along
a valve’s rotation axis, the end-effector can manipulate valves of up to 160
millimeters in diameter by grasping them from the outside. Additionally,
larger valves with gaps between the center of rotation and handle can be
grasped inside the valve’s handle.

Figure 4: Autonomous Mapping results. The 3D maps resulting from online
mapping are visualized as pointclouds. a), b) and c) correspond to scenarios
1, 2 and 3. The legend indicating map height using the colormap applies to all
maps. In c), improperly mapped geometry is highlighted by two black circles.
Figure 5: Radiation mapping results. The 2D map (black) is overlaid with the
robot’s trajectory (green) and the estimated radiation map (yellow shading)
during scenario 1, depicted in a) and scenario 3, depicted in b). Red borders
in b) indicate map truncation to zoom into relevant parts of the map.

##  V EXPERIMENTS

Our robot was evaluated during trials in the Tritolwerk nbc and disaster
relief training area as well as the enrich European Robotics Hackathon. The
Tritolwerk nbc and disaster relief training area, situated in Eggendorf,
Austria, serves as a training ground for the Austrian Armed Forces and various
emergency response organizations. It includes ruins and bunker facilities,
which serve as our setting for simulating diverse disaster scenarios. enrich
is hosted at the decommissioned nuclear power plant in Zwentendorf, Austria
and aims to simulate catastrophic events and evaluate robots’ ability to
navigate hazardous environments and locate survivors.

###  V-A Experimental Setup

We test the robot’s ability to perform the four tasks outlined in Sec. I,
which revolve around mapping the geometry and radiation in the robot’s
environment, taking samples of substances for the Raman spectrometer to
analyze and closing valves in a semi-autonomous manner. The following chapters
present experiments conducted to demonstrate the robot’s capability to perform
these tasks in the following scenarios:

  * •

Scenario 1: Tritolwerk main production hall
This test location is two stories tall and features pillars and scaffolding
within the room. During the trial, a radiation source is placed near one of
the hall’s exits. The radiation source is partly occluded using bricks to test
the radiation mapping’s ability to estimate the spread of radioactive
emissions from the discrete radiation measurements.

  * •

Scenario 2: Indoor navigation in nuclear power plant
(enrich training area)
This test consists of navigation through hallways and the power plant’s
turbine rooms. The area features narrow hallways connecting multiple rooms,
resulting in more rooms compared to scenario 1. However, no radiation source
is present, therefore, no radiation map is plotted in Fig. 5.

  * •

Scenario 3: Radiation source pattern
During this scenario, the robot is teleoperated to move towards an
unobstructed radiation source in an open space. The robot’s trajectory follows
a flower-shaped pattern to imitate how human disaster responders triangulate
radiation sources. The scenario is located outside near the nuclear power
plant Zwentendorf.

###  V-B Autonomous Mapping

Fig. 4 plots the 3D maps generated from all three scenarios. In scenario 1
(Fig. 4a), details, such as the scaffolding in the middle of the room, belts
used to secure the scaffolding as well as the barrier at the top entrance of
the hall are properly captured and visualized by the system. The neighboring
rooms are incomplete since the mission only involves the main hall.

Scenario 2’s area (see Fig. 4b) consists of narrow corridors, which results in
a denser pointcloud compared to scenario 1. The hallways form a loop, where
loop closure is correctly applied by the 3D mapper. The turbine rooms on the
bottom of the plot are not fully shown in the map, because the robot is too
wide to navigate into these rooms.

Scenario 3 (see Fig. 4c) shows an open outside area. Contrary to Fig. 4a) and
Fig. 4b), the resulting map shows deficiencies highlighted using black
circles. In these areas, the same geometry has been mapped multiple times due
to a drift in odometry. This behavior is most likely caused by the direct
sunlight experienced during the trials, which introduces significant noise
into the lidar measurements. Similar to scenarios 2 and 3, [8] also provide
maps of the nuclear power plant Zwentendorf recorded during enrich, however,
they omit the section depicted in Fig. 4c).

###  V-C Radiation Mapping

Fig. 5 presents the 2D map used for navigation (black) and the robot
trajectory (green) overlaid onto the radiation map (yellow shading) for
scenario 1 (Fig. 5a) and scenario 3 (Fig. 5b). The resulting heatmap
represents a triangulation of each radiation source based on discrete
radiation measurements within the 2D map. Fig. 5a) shows the dampening of
radioactive emissions to the top right of the radiation source. This dampening
corresponds to the locations where bricks were placed to partially occlude the
radiation source. The trajectory shown in Fig. 5b) shows the robot circling
the radiation source at different distances. This teleoperated trajectory aims
to imitate how human first responders manually triangulate the position of
radiation sources by circling the estimated position of the radiation source
in a pattern resembling a flower. Scenario 2 is omitted from Fig. 5 because
there is no radiation source present during the trial.

###  V-D Semi-Autonomous Sampling

The goal of the semi-autonomous sampling procedure is to collect solid powder-
like as well as fluid substances using one of the stored probes and to analyze
them on-board. Fig. 6 presents a series of autonomous routines that are called
sequentially by the operator to sample a powdered substance in front of the
robot. The process includes a) preparation to get a probe from storage, b)
probe grasping, c) swiping of the probe through the sampled substance, d)
preparation to place the sample, and e) probe placement in the analysis tray.
The robot arm trajectory between the poses presented in Fig. 6 is planned in a
manner that prevents direct fly-over of the contaminated probe. This approach
prevents the sampled substance from trickling from the probe onto the robot
chassis while it is moved from the sampling pose to the analysis tray.

a)

b)

c)

d)

e)

Figure 6: Semi-autonomous sampling sequence. The sampling process consists of
the following poses: a) ready to get probe from storage, b) grasp probe, c)
swipe substance, d) ready to place sample, and e) sample placed in analysis
tray.

###  V-E Teleoperated Valve Manipulation

Fig. 7 presents the robot performing the valve closing task in scenario 1. The
sar robot manipulates a tank’s 100 millimeter output valve by gripping the
valve from the outside.

Figure 7: Valve manipulation. The sar mobile robot manipulates a tank’s output
valve at Tritolwerk, Eggenburg, Austria. The gripper is closed partially to
match the valve’s 100 mm outer diameter ensuring proper grip of the handle.

##  VI CONCLUSIONS

In this paper, we present a semi-autonomous sar (sar) mobile robot system,
aimed at assisting first responders in disaster scenarios. It fulfils four
main tasks in one combined solution: (1) autonomous exploration and mapping of
unknown areas, (2) triangulation of a map from discrete radiation
measurements, (3) sampling and analysis of encountered substances in a semi-
autonomous manner, and (4) teleoperated manipulation of valves that control
the flow of hazardous substances. The system is based on a tracked robot base
with a five dof (dof) robot arm for manipulation and uses 2D, 3D and radiation
mapping pipelines for autonomous navigation as well as visualization of the
robot environment to the human operators. Semi-autonomous robot navigation is
achieved using priority-based multiplexing of control signals, while the
sampling and valve manipulation tasks are solved using a combination of
prerecorded routines and teleoperation. The system was evaluated at a disaster
response training site as well as at the enrich (enrich) to show its ability
to perform the required tasks.

Future work on the presented system intends to improve perception
capabilities, especially in direct sunlight by tightly coupling the camera’s
image data to the lidar (lidar) odometry and mapping pipeline. Furthermore,
the level of autonomy during valve manipulation can be increased by
incorporating methods for vision-guided arm control, such as visual servoing,
into the robot control pipeline.

## References

  * [1] R. R. Murphy, “Trial by fire [rescue robots],” _IEEE Robotics & Automation Magazine_, vol. 11, no. 3, pp. 50–61, 2004.
  * [2] K. Ohno, S. Kawatsuma, T. Okada, E. Takeuchi, K. Higashi, and S. Tadokoro, “Robotic control vehicle for measuring radiation in fukushima daiichi nuclear power plant,” in _2011 IEEE International Symposium on Safety, Security, and Rescue Robotics_. IEEE, 2011, pp. 38–43.
  * [3] K. Berns, A. Nezhadfard, M. Tosa, H. Balta, and G. D. Cubber, “Unmanned ground robots for rescue tasks,” in _Search and Rescue Robotics_. Rijeka: IntechOpen, 2017, ch. 4.
  * [4] H. Chitikena, F. Sanfilippo, and S. Ma, “Robotics in search and rescue (sar) operations: An ethical and design perspective framework for response phase,” _Applied Sciences_ , vol. 13, no. 3, 2023.
  * [5] G. Lunghi, R. Marin, M. Di Castro, A. Masi, and P. J. Sanz, “Multimodal human-robot interface for accessible remote robotic interventions in hazardous environments,” _IEEE Access_ , vol. 7, pp. 127 290–127 319, 2019.
  * [6] M. P. Manuel, M. Faied, and M. Krishnan, “A novel lora lpwan-based communication architecture for search & rescue missions,” _IEEE Access_ , vol. 10, pp. 57 596–57 607, 2022.
  * [7] D. H. Stolfi, M. R. Brust, G. Danoy, and P. Bouvry, “Uav-ugv-umv multi-swarms for cooperative surveillance,” _Frontiers in Robotics and AI_ , vol. 8, 2021.
  * [8] J. Bedkowski, “Open source, open hardware hand-held mobile mapping system for large scale surveys,” _SoftwareX_ , vol. 25, p. 101618, 2024.
  * [9] K. Groves, E. Hernandez, A. West, T. Wright, and B. Lennox, “Robotic exploration of an unknown nuclear environment using radiation informed autonomous navigation,” _Robotics_ , vol. 10, no. 2, 2021.
  * [10] R. Guzman, R. Navarro, J. Ferre, and M. Moreno, “Rescuer: Development of a modular chemical, biological, radiological, and nuclear robot for intervention, sampling, and situation awareness,” _Journal of Field Robotics_ , vol. 1–15 (2015), 06 2015.
  * [11] J. P. Queralta, J. Taipalmaa, B. Can Pullinen, V. K. Sarker, T. Nguyen Gia, H. Tenhunen, M. Gabbouj, J. Raitoharju, and T. Westerlund, “Collaborative multi-robot search and rescue: Planning, coordination, perception, and active vision,” _IEEE Access_ , vol. 8, pp. 191 617–191 643, 2020.
  * [12] F. G. Rosas, F. Hoeller, and F. E. Schneider, “Autonomous robot exploration and measuring in disaster scenarios using behavior trees,” in _2020 IEEE 10th International Conference on Intelligent Systems (IS)_ , 2020, pp. 469–474.
  * [13] N. A. Abd Rahman, K. S. M. Sahari, N. A. Hamid, and Y. C. Hou, “A coverage path planning approach for autonomous radiation mapping with a mobile robot,” _International Journal of Advanced Robotic Systems_ , vol. 19, no. 4, 2022.
  * [14] B. Nouri Rahmat Abadi, A. West, M. Nancekievill, C. Ballard, B. Lennox, O. Marjanovic, and K. Groves, “Carma ii: A ground vehicle for autonomous surveying of alpha, beta and gamma radiation,” _Frontiers in Robotics and AI_ , vol. 10, p. 1137750, 2023.
  * [15] N. A. Abd Rahman, K. S. M. Sahari, M. F. A. Jalal, A. A. Rahman, M. I. Abd Adziz, and M. Z. Hassan, “Mobile robot for radiation mapping in indoor environment,” in _IOP Conference Series: Materials Science and Engineering_ , vol. 785, no. 1. IOP Publishing, 2020, p. 012021.
  * [16] A. West, T. Wright, I. Tsitsimpelis, K. Groves, M. J. Joyce, and B. Lennox, “Real-time avoidance of ionising radiation using layered costmaps for mobile robots,” _Frontiers in Robotics and AI_ , vol. 9, p. 862067, 2022.
  * [17] F. G. Rosas, F. Hoeller, and F. E. Schneider, “Autonomous robot exploration and measuring in disaster scenarios using behavior trees,” in _2020 IEEE 10th International Conference on Intelligent Systems (IS)_. IEEE, 2020, pp. 469–474.
  * [18] G. Christie, A. Shoemaker, K. Kochersberger, P. Tokekar, L. McLean, and A. Leonessa, “Radiation search operations using scene understanding with autonomous uav and ugv,” _Journal of Field Robotics_ , vol. 34, no. 8, pp. 1450–1468, 2017.
  * [19] D. Connor, P. G. Martin, and T. Scott, “Airborne radiation mapping: overview and application of current and future aerial systems,” _International Journal of Remote Sensing_ , vol. 37, pp. 5953–5987, 12 2016.
  * [20] J. Zhao, J. Gao, F. Zhao, and Y. Liu, “A search-and-rescue robot system for remotely sensing the underground coal mine environment,” _Sensors_ , vol. 17, no. 10, p. 2426, 2017.
  * [21] A. Denker and M. C. İşeri, “Design and implementation of a semi-autonomous mobile search and rescue robot: Salvor,” in _2017 International Artificial Intelligence and Data Processing Symposium (IDAP)_. IEEE, 2017, pp. 1–6.
  * [22] G. A. Novotny and W. Kubinger, “Design and implementation of a mobile search and rescue robot,” in _Proceedings of the Joint Austrian Computer Vision and Robotics Workshop_ , 2020, pp. 21–26.
  * [23] M. Quigley, K. Conley, B. P. Gerkey, J. Faust, T. Foote, J. Leibs, R. Wheeler, and A. Y. Ng, “Ros: an open-source robot operating system,” in _ICRA Workshop on Open Source Software_ , 2009.
  * [24] U. Jahn, D. Heß, M. Stampa, A. Sutorma, C. Röhrig, P. Schulz, and C. Wolff, “A taxonomy for mobile robots: Types, applications, capabilities, implementations, requirements, and challenges,” _Robotics_ , vol. 9, no. 4, 2020.
  * [25] S. Macenski and I. Jambrecic, “Slam toolbox: Slam for the dynamic world,” _Journal of Open Source Software_ , vol. 6, no. 61, p. 2783, 2021.
  * [26] T. Shan, B. Englot, D. Meyers, W. Wang, C. Ratti, and R. Daniela, “Lio-sam: Tightly-coupled lidar inertial odometry via smoothing and mapping,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 5135–5142.
  * [27] D. Lee, M. Jung, W. Yang, and A. Kim, “Lidar odometry survey: recent advancements and remaining challenges,” _Intelligent Service Robotics_ , Feb 2024.
  * [28] D. V. Lu, D. Hershberger, and W. D. Smart, “Layered costmaps for context-sensitive navigation,” in _2014 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2014, pp. 709–715.
  * [29] S. Pütz, J. Santos Simón, and J. Hertzberg, “Move base flex a highly flexible navigation framework for mobile robots,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. Madrid, Spain: IEEE, 01.-05. October 2018, pp. 3416–3421.
  * [30] A. Abdelhakim, “Machine learning for localization of radioactive sources via a distributed sensor network,” _Soft Comput._ , vol. 27, no. 15, p. 10493–10508, 2023.
  * [31] A. West, I. Tsitsimpelis, M. Licata, A. Jazbec, L. Snoj, J. Malcolm, and B. Lennox, “Use of gaussian process regression for radiation mapping of a nuclear reactor with a mobile robot,” _Scientific Reports_ , vol. 11, no. 1, pp. 1–11, 2021.
  * [32] C. E. Rasmussen and C. K. I. Williams, _Gaussian processes for machine learning._ MIT Press, 2006.
  * [33] G. Novotny, S. Schwaiger, L. Muster, M. Aburaia, and W. Wöber, “On the applicability of docker containers and systemd services for search and rescue applications,” in _Proceedings of the Austrian Robotics Workshop_. Johannes Kepler University, 4 2023, pp. 19–24.


# Thermo-optic hysteresis with bound states in the continuum

D. N. Maksimov1,2,3, A. S. Kostyukov1, A. E. Ershov1,3, M. S. Molokeev1,2, E.
N. Bulgakov2,3 and V. S. Gerasimov1,3 1IRC SQC, Siberian Federal University,
660041, Krasnoyarsk, Russia  2Kirensky Institute of Physics, Federal Research
Center KSC SB RAS, 660036, Krasnoyarsk, Russia  3Institute of Computational
Modelling SB RAS, Krasnoyarsk, 660036, Russia

###### Abstract

We consider thermo-optic hysteresis in a silicon structure supporting bound
state in the continuum. Taking into account radiative heat transfer as a major
cooling mechanism we constructed a non-linear model describing the optical
response. It is shown that the thermo-optic hysteresis can be obtained with
low intensities of incident light $I_{0}\approx 1~{}\rm{W/m^{2}}$ at the red
edge of the visible under the critical coupling condition.

##  I Introduction

Recently, we have seen a surge of interest to bound states in the continuum
(BICs) Hsu  _et al._ (2016); Koshelev  _et al._ (2019a); Sadreev (2021) that
have grown to an important tool in nanophotonics paving a way to optical
devices with enhanced light-matter interaction. BICs do not couple the
incident light, however, if the symmetry of the system is broken the BICs are
observed as narrow Fano resonances in the scattering spectrum Kim  _et al._
(1999); Shipman and Venakides (2005); Sadreev  _et al._ (2006); Blanchard  _et
al._ (2016). In more detail, the BICs are spectrally surrounded by a leaky
band of high-quality resonances (quasi-BICs) which can be excited from the
far-field Yuan and Lu (2017). The excitation of the strong resonances results
in critical field enhancement Yoon  _et al._ (2015); Mocella and Romano (2015)
with the near-field amplitude controlled by the frequency and the angle of
incidence of the incoming monochromatic wave. The critical field enhancement
triggers nonlinear optical effects even with a low intensity of the incident
light. This resonant enhancement of nonlinear effects can lead to symmetry
breaking Bulgakov  _et al._ (2011), channel dropping Bulgakov  _et al._
(2013), excitation of non-linear standing waves Yuan and Lu (2016) as well as
self-adaptive robust Bulgakov and Sadreev (2014) and tunable Fabry-Perot
Bulgakov and Sadreev (2010) BICs.

In the field of nonlinear optics the BICs have been applied for second
harmonic generation Ndangali and Shabanov (2013); Koshelev  _et al._ (2019b);
Wang and Zhang (2018). Quasi-BICs in subwavelength dielectric resonators Rybin
_et al._ (2017); Bogdanov  _et al._ (2019) are also shown to demonstrate to
enhance second harmonic generation Carletti  _et al._ (2018); Koshelev  _et
al._ (2020). In Yuan and Lu (2020) it was found that, otherwise decoupled,
BICs can be excited via second harmonic generation by illuminating the
structure from the far field. At the same time it was shown theoretically that
quasi-BICs allow for optical bistability Yuan and Lu (2017); Krasikov  _et
al._ (2018); Bulgakov and Maksimov (2019); Maksimov  _et al._ (2020) due to
the Kerr effect.

More recently, the research focus shifted towards BICs in lossy structures. In
the presence of material absorption the BIC are shown to acquire finite-life,
albeit remain localized and decoupled from the outgoing channels Hu  _et al._
(2020). Quasi-BIC in lossy periodic structures are found to be instrumental
for enhancement of light absorption Saadabad  _et al._ (2021); Tan  _et al._
(2022) in the critical coupling regime even in low loss dielectrics. This
opens novel opportunities for highly efficient light absorbers Zhang and Zhang
(2015); Wang  _et al._ (2020); Sang  _et al._ (2021); Xiao  _et al._ (2021);
Cai  _et al._ (2022).

Lately, it has been suggested that thermo-optical effects can be the
dominating nonlinear effects in BIC supporting structures due to heating by
absorbed radiationZograf  _et al._ (2021). In this work we investigate
resonantly enhanced thermo-optical bistability Sun and Reano (2010); Khandekar
and Rodriguez (2017); Gao  _et al._ (2017); Pottier and Bellon (2021); Ryabov
_et al._ (2022) in a system supporting an optical BICs. Taking into account
radiative heat transfer as a major cooling mechanism we shall construct a non-
linear model based on the temporal coupled mode theory (TCMT) Fan  _et al._
(2003) and theoretically demonstrate thermo-optic hysteresis.

##  II Scattering theory

We consider an array of identical dielectric rods of radius $R_{0}$ linearly
arranged with period $L$ in vacuum. The rods are made of amorphous silicon.
The axes of the rods are collinear and aligned with the $z$-axis as shown in
Fig. 1 (a). Such a system is known to support an abundance of BICs as
demonstrated in Bulgakov and Sadreev (2015); Bulgakov and Maksimov (2017a). In
this work we performed numerical simulations with application of FDTD
Lumerical to examine the properties of the BIC induced optical response taking
into account both temperature and frequency dependence of the refractive
index. In Fig. 1 (b) we compare the eigenmode profile of an optical BIC with
vacuum wavelength $\lambda_{\scriptscriptstyle{\rm BIC}}=782~{}\rm{nm}$
against the scattering solution obtained under illumination by a TE plane wave
with $\lambda=785~{}\rm{nm}$ at the incidence angle
$\theta=8.97~{}\mathrm{deg}$, and the incident wave vector on the $x0y$-plane
as shown in Fig. 1 (a). The results are obtained with the numerical values of
the refractive index from Jellison and Modine (1994). In Fig. 1 (b) one can
see a striking similarity between the two field profiles.

By definition a BIC can not couple to incident light. However, each BIC is a
singular point on the dispersion sheet of a leaky band where the $Q$-factor
diverges to infinity. The BIC is, thus, spectrally surrounded by a family of
high-$Q$ leaky modes (quasi-BICs) with the resonant $Q$-factor controlled by
variation of the angle of incidence. In our case the BIC occurs in the
$\Gamma$-point being symmetrically mismatched from the single radiation
channel of the zeroth diffraction order. Therefore, the high-$Q$ resonant
response is triggered by setting off the angle of incidence from zero. Thus,
the similarity between the BIC and the scattering solution is explained by
both BIC and quasi-BIC sitting in the same dispersion band.

Our goal in this section is to set up the optimal regime for enhanced light
absorption leading to the most significant change of the refractive index by
heating. The problem of enhanced absorption by quasi-BICs has been previously
considered in the literature Saadabad  _et al._ (2021). The central result for
upside-down symmetric structures, such as the one in Fig. 1 (a), is that the
maximal absorption of $50\%$ can be achieved in the critical coupling point
where the radiation and absorption loss rates are equal to one another.

Figure 1: (a) Set-up of the array of dielectric rods. (b) Comparison between
the eigenmode profile for a BIC with $\lambda_{\scriptscriptstyle{\rm
BIC}}=782~{}\rm{nm}$ and scattering solution under illumination by a plane
wave with $\lambda=785~{}\rm{nm}$. The solutions are visualized as the
$z$-components of the electric field. The geometric parameters:
$R_{0}=128~{}{\rm{nm}}$, $L=428~{}\rm{nm}$.

The BIC induced optical response can be understood in the framework of
temporal coupled mode theory (TCMT) Fan  _et al._ (2003). For the reader’s
convenience we list the most important TCMT formulas below. Let us consider
two-channel scattering for TE-polarized light. The $S$-matrix is implicitly
defined through the following equation

$\left(\begin{array}[]{c}s_{1}^{\scriptscriptstyle(\\!-\\!)}\\\
s_{2}^{\scriptscriptstyle(\\!-\\!)}\\\
\end{array}\right)=\widehat{S}\left(\begin{array}[]{c}s_{1}^{\scriptscriptstyle(\\!+\\!)}\\\
s_{2}^{\scriptscriptstyle(\\!+\\!)}\\\ \end{array}\right),$ (1)

where $s_{m}^{\scriptscriptstyle(\\!\pm\\!)}$ are the amplitudes of the plane
waves in the far-field with subscript $m=1,2$ corresponding to the upper and
lower half-spaces while superscripts
${}^{\scriptscriptstyle(\\!+\\!)},^{\scriptscriptstyle(\\!-\\!)}$ stand for
incident and outgoing waves respectively. We assume that the system is
illuminated by a monochromatic wave of frequency $\omega$. In what follows we
introduce the vectors of incident and outgoing amplitudes
$|s^{\scriptscriptstyle(\\!\pm\\!)}(t)\rangle$ which oscillate in time with
the harmonic factor $e^{-i\omega t}$. The TCMT equations Fan  _et al._ (2003)
take the following form

$\displaystyle\frac{da(t)}{dt}=-(i\omega_{0}+\gamma+\gamma_{0})a(t)+\langle
d^{*}|s^{\scriptscriptstyle(\\!+\\!)}\rangle,$
$\displaystyle|s^{\scriptscriptstyle(\\!-\\!)}\rangle=\widehat{C}|s^{\scriptscriptstyle(\\!+\\!)}\rangle+a(t)|d\rangle,$
(2)

where $\widehat{C}$ is the matrix of direct (non-resonant) process,
$\omega_{0}$ is the resonance center frequency, $\gamma$ is the radiation
decay rate, $\gamma_{0}$ is material loss decay rate, $a$ is the amplitude of
the resonant eigenmode, and $|d\rangle$ is the $2\times 1$ vector of coupling
constants satisfying

$\langle d|d\rangle=2\gamma.$ (3)

The solution for the $S$-matrix reads

$\widehat{S}=\widehat{C}+\frac{|d\rangle\langle
d^{*}|}{i(\omega_{0}-\omega)+\gamma+\gamma_{0}}.$ (4)

In the case of the center-plane mirror symmetry we have

$\widehat{C}=e^{i\phi}\left(\begin{array}[]{cc}\rho&i\tau\\\
i\tau&\rho\end{array}\right),$ (5)

where $\phi,\rho,\tau$ are real valued parameters such as

$\rho^{2}+\tau^{2}=1.$ (6)

The following equation also holds true Fan  _et al._ (2003)

$\widehat{C}|d^{*}\rangle=-|d\rangle.$ (7)

Using Eq. (3) together with Eq. (7) one finds

$|d\rangle=e^{i\frac{\phi}{2}}\sqrt{\frac{\gamma}{2(1+\rho)}}\left(\begin{array}[]{l}\tau-i(1+\rho)\\\
\tau-i(1+\rho)\end{array}\right),$ (8)

for symmetric modes, and

$|d\rangle=e^{i\frac{\phi}{2}}\sqrt{\frac{\gamma}{2(1+\rho)}}\left(\begin{array}[]{l}\tau+i(1+\rho)\\\
-\tau-i(1+\rho)\end{array}\right),$ (9)

for anti-symmetric ones. The reflectance can be found from Eq. (4) as

$R=\frac{\rho^{2}(\omega-\omega_{0})^{2}\mp
2\rho\tau\gamma(\omega-\omega_{0})+\tau^{2}\gamma^{2}+\rho^{2}\gamma_{0}^{2}}{(\omega-\omega_{0})^{2}+(\gamma+\gamma_{0})^{2}},$
(10)

while the transmittance is given by

$T=\frac{\tau^{2}(\omega-\omega_{0})^{2}\pm
2\rho\tau\gamma(\omega-\omega_{0})+\rho^{2}\gamma^{2}+\tau^{2}\gamma_{0}^{2}}{(\omega-\omega_{0})^{2}+(\gamma+\gamma_{0})^{2}},$
(11)

where the top sign is used for symmetric modes. Finally, the absorbance is
given by

$A=\frac{2\gamma\gamma_{0}}{(\omega-\omega_{0})^{2}+(\gamma+\gamma_{0})^{2}}.$
(12)

The maximal absorbance is obtained at the critical coupling condition
$\gamma=\gamma_{0}$ and $\omega=\omega_{0}$. Substituting the above into Eq.
(12) one finds

$A_{\mathrm{max}}=\frac{1}{2}.$ (13)

Figure 2: Scattering spectra for the system shown in Fig. 1 (a),
$R_{0}=128~{}{\rm{nm}}$, $L=428~{}\rm{nm}$. (a) The reflectance,
transmittance, and absorbance spectra found from FDTD simulations at the
critical coupling angle $\theta=8.97~{}\mathrm{deg}$. (b) The TCMT
approximation of the data in subplot (c); $L\omega/c=3.425$,
$L\gamma/c=L\gamma_{0}/c=0.0056$, $\rho=0$.

Since the system supports its antisymmetric BIC in the $\Gamma$-point, we
expect that the resonant frequency $\omega_{0}$ and the radiation decay rate
$\gamma$ are given by the following equations Bulgakov and Maksimov (2017b)

$\displaystyle\omega_{0}=\omega_{\scriptscriptstyle\rm{BIC}}+\kappa_{\omega}\theta^{2}+\mathcal{O}(\theta^{4}),$
$\displaystyle\gamma=\kappa_{\gamma}\theta^{2}+\mathcal{O}(\theta^{4}),$ (14)

where $\kappa_{\omega,\gamma}$ are the coefficients of the Taylor expansion.
At the same time the absorption decay rate can be assessed in the perturbative
manner Vainstein (1988) as

$\gamma_{0}=\epsilon^{\prime\prime}\mkern-1.2mu\int\limits_{S_{R}}dS\frac{{\bf
E}^{\dagger}{\bf E}+{\bf H}^{\dagger}{\bf H}}{4},$ (15)

where $\epsilon^{\prime\prime}\mkern-1.2mu$ and ${\bf E}$ and ${\bf H}$ are
the electric and magnetic vectors of the BIC eigenfield satisfying the
normalization condition

$\int\limits_{S}dS\frac{\epsilon(x,y){\bf E}^{\dagger}{\bf E}+{\bf
H}^{\dagger}{\bf H}}{4}=1,$ (16)

with $S$ as the area of the elementary cell, and $S_{R}$ as the area of the
rod cross section. By a close examination of Eq. (II) and Eq. (15) one
concludes that the critical coupling condition can always be achieved by
increasing the angle of incidence. In our simulations we found that the
critical coupling condition is fulfilled at $\theta=8.97~{}\mathrm{deg}$ which
is the case shown in Fig. 1 (b). In Fig. 2 (a) we show transmittance,
reflectance, and absorbance at the critical angle found via direct numerical
simulations. In Fig. 2 (b) we plotted the TCMT approximation with Eqs. (10,
11, 12). The non-resonant reflection coefficient is taken $\rho=0$ according
to Bulgakov and Maksimov (2019). The other approximation parameters are listed
in the caption to Fig. 2. By comparing Fig. 2 (a) against Fig. 2 (b) one can
see a certain deviation from the TCMT approach which is due to the
perturbative nature of the solution given by Eqs. (10, 11, 12). In particular,
the absorption in the direct process is not accounted for by Eq. (5). A better
accuracy of the perturbative solution can be achieved in the near infrared
where the imaginary part of the refractive index drops by two orders of
magnitude Schinke  _et al._ (2015). In Fig. 3 we compare the numerical data
for dielectric array with $L=580~{}\mathrm{nm},R_{0}=174~{}{\rm{nm}}$
supporting a BIC at wavelength
$\lambda_{\scriptscriptstyle\rm{BIC}}=1027~{}\rm{nm}$. One can see an
excellent coincidence between the TCMT and numerical simulations. Also notice
that the position of the critical coupling spot in the absorption spectrum is
shifted towards smaller angles of incidence due to a smaller $\gamma_{0}$.
Further on in this work we stay in the visible relying on approximation shown
in Fig. 2 (b), and bearing in mind that the accuracy of the TCMT can be
improved by going to the near infrared.

Figure 3: Absorbance, reflectance, and transmittance in the spectral vicinity
of the BIC in the near infrared,
$\lambda_{\scriptscriptstyle\rm{BIC}}=1027~{}\rm{nm}$,
$L=580~{}\mathrm{nm},R_{0}=174~{}{\rm{nm}}$. Column (a) shows the data
obtained by direct FDTD simulations. Column (b) shows the TCMT approximation,
$L\omega_{\scriptscriptstyle\mathrm{BIC}}/c=3.549$, $\rho=0$,
$L\kappa_{\omega}/c=-0.00016$, $L\kappa_{\gamma}/c=0.000072$.

##  III Heat transfer

Since our system consists of dielectric rods in vacuum, the only mechanism of
heat transfer is radiation from the rods. The total energy lost per unit
surface area across all wavelengths per unit time is given by the Stefan-
Boltzmann law

$j=\alpha\sigma[(\Theta_{0}+\Delta\Theta)^{4}-\Theta_{0}^{4}],$ (17)

where $\alpha$ is the emissivity, $\sigma=5.670\times
10^{-8}~{}\mathrm{W}\cdot\mathrm{m}^{-2}\cdot\mathrm{K}^{-4}$ is the Stefan-
Boltzmann constant, $\Theta_{0}$ is the temperature of the surrounding medium
and $\Delta\Theta$ is the temperature increase due to absorption of light. The
emissivity, $\alpha$, is a quantity of key importance, that describes the
difference in the radiative properties between the dielectric rod and the
black body. In this work we assessed the emissivity by computing the
absorption coefficient of a single dielectric rod in the temperature range
$25^{\circ}\mathrm{C}\leq\Theta\leq 300^{\circ}\mathrm{C}$. The absorption
coefficient was found as the average over all possible direction of the
incident wave vector for both polarizations and, then, averaged with the
Planck’s distribution

Figure 4: Emissivity of a single dielectric rod as a function of temperature;
$L=428$ nm, $R_{0}=128$ nm. The average emissivity, $\tilde{\alpha}$, is shown
by the horizontal line.

$p(\omega)=\frac{\hbar\omega^{3}}{4\pi^{3}c^{2}}\frac{1}{e^{\hbar\omega/(k_{\mathrm{B}}\Theta)}-1}$
(18)

peaked at about $10~{}\mu{\mathrm{m}}$. The numerical values of the refractive
index were taken from Chandler-Horowitz and Amirtharaj (2005). The temperature
dependance of the refractive index was neglected in the mid-infrared. The
result is demonstrated in Fig. 4. One can see that the emissivity is dependant
on temperature of the rods due to the shift of the peak of Planck’s
distribution. For further convenience we computed the temperature average
$\tilde{\alpha}=2.8\times 10^{-5}$ which will be applied in our simulations.
The average emissivity is shown in Fig. 4 as a horizontal line.

The regime of heating can be determined by the value of the Biot number
Bergman  _et al._ (2011)

$\mathrm{Bi}=\frac{h_{S}}{k}L_{0},$ (19)

where $L_{0}$ is the characteristic length, $h_{S}$ is the surface heat
transfer coefficient, and $k$ is the thermal conductivity of the body. Values
of the Biot number smaller than $0.1$ imply that the heat conduction inside
the body is much faster than the heat lost at the interface between the body
and the surrounding medium. In this regime we can assume that the temperature
is constant across the body. By using Eq. (17) we find

$h_{S}=4\tilde{\alpha}\sigma\Theta_{0}^{3}\Delta\Theta.$ (20)

Taking $L_{0}$ as the perimeter of the rod, and
$k=1.8~{}\mathrm{W}/(\mathrm{m}\cdot\mathrm{K})$ according to Wada and Kamijoh
(1996) we find at the room temperature $\mathrm{Bi}=6.5\cdot 10^{-11}$. This
allows us to assume that the temperature is uniform in the cross-section of
the rod and avoid simulating heat transfer within the rods. Now our TCMT and
heat transfer models can be linked by equating the energy absorption rate in a
unit length of the rod

$\frac{dE}{dt}=\gamma_{0}(\Theta)|a|^{2}$ (21)

to the radiation rate from the surface per unit length of the same rod

$\frac{dE}{dt}=2\pi
R_{0}\tilde{\alpha}\sigma[(\Theta_{0}+\Delta\Theta)^{4}-\Theta_{0}^{4}].$ (22)

Figure 5: Real $n^{\prime}$ and imaginary $n^{\prime\prime}$ parts of the
refractive index of the amorphous silicon as a function of temperature at
$\lambda=782~{}\mathrm{nm}$. Solid lines show linear approximation.

##  IV Non-linear TCMT equation and optical bistability

Figure 6: Thermo-optic hysteresis in the system shown in Fig. 1 (a) with the
same values of geometric parameters at incident angle $\theta=8.97^{\circ}$.
(a) Temperature increment against frequency of incident light,
$I_{0}=0.54\rm{W/m^{2}}$. (c) Scattering spectra, $I_{0}=0.54\rm{W/m^{2}}$.
(b) Temperature increment against the intensity of incident light with vacuum
wavelength $\lambda=795.5\rm{nm}$. (d) Absorbance, reflectance, and
transmittance against the intensity, $\lambda=795.5\rm{nm}$.

Frequency domain coupled mode equation reads

$[i\omega_{0}(\Theta)-i\omega+\gamma(\Theta)+\gamma_{0}(\Theta)]a=d_{1}(\Theta)\sqrt{LI_{0}},$
(23)

where $I_{0}$ is the intensity of the incident monochromatic light. Notice
that according to Eq. (16) the BIC eigenfield is normalized to store a unit
energy in the elementary cell per unit length of the rod. Thus, the intensity
in Eq. (23) is multiplied by the period under the basic TCMT assumption that
the scattering channel function are normalized to supply (retract) a unit
energy in a unit of time Fan  _et al._ (2003). All parameters including the
resonant eigenfrequency, absorption and radiation decay rates are dependent on
temperature. Hence we write

$\displaystyle\omega_{0}(\Theta)=\omega_{0}^{\scriptscriptstyle{(0)}}+\omega_{0}^{\scriptscriptstyle{(1)}}\Delta\Theta,$
$\displaystyle\gamma_{0}(\Theta)=\gamma_{0}^{\scriptscriptstyle{(0)}}+\gamma_{0}^{\scriptscriptstyle{(1)}}\Delta\Theta,$
$\displaystyle\gamma(\Theta)_{=}\gamma^{\scriptscriptstyle{(0)}}+\gamma^{\scriptscriptstyle{(1)}}\Delta\Theta.$
(24)

As it was mentioned the temperature dependence is due to the temperature
dependance of the refractive index. By using the experimental data from
Jellison and Modine (1994) we plotted the temperature dependance of both real
$n^{\prime}$ and imaginary $n^{\prime\prime}$ parts of the refractive index of
the amorphous silicon at $\lambda=782~{}\mathrm{nm}$ in Fig. 5. One can see in
Fig. 5 that both quantities can be linearly approximated in the temperature
range of interest.

By numerically computing the scattering spectra with
$\theta=8.97~{}\mathrm{deg}$ at elevated temperatures and performing a TCMT
fit we found the numerical values of the expansion coefficients in Eq. (IV) as
follows

$\displaystyle L\omega_{0}^{\scriptscriptstyle{(1)}}/c=-2.80\times
10^{-4},~{}\mathrm{K}^{-1},$ $\displaystyle
L\gamma_{0}^{\scriptscriptstyle{(1)}}/c=2.67\times
10^{-5},~{}\mathrm{K}^{-1},$ $\displaystyle
L\gamma^{\scriptscriptstyle{(1)}}/c=6.67\times 10^{-6},~{}\mathrm{K}^{-1}.$
(25)

One can see from the above data that the dominating effect is a shift of the
resonant frequency. The resonant frequency is at least one order of magnitude
more sensitive to heating than the other TCMT parameters. Also notice that the
absorption loss rate is much more sensitive to heating than the radiative loss
rate. This can be understood from the explicit linear dependance of
$\gamma_{0}$ on the imaginary part of the refractive index given by Eq. (15).
In what follows we neglect the temperature dependance of $\gamma$ and,
consequently, of $|\bf{d}\rangle$. From (23) and (3) we find

$|a|^{2}=\frac{\gamma
I_{0}L}{[\omega-\omega_{0}(\Theta)]^{2}+[\gamma+\gamma_{0}(\Theta)]^{2}},$
(26)

which in combination with Eq. (21) and Eq. (22) yields

$\displaystyle\frac{(\gamma_{0}^{\scriptscriptstyle{(0)}}+\gamma_{0}^{\scriptscriptstyle{(1)}}\Delta\Theta)\gamma
I_{0}L}{(\omega-\omega_{0}^{\scriptscriptstyle{(0)}}-\omega_{0}^{\scriptscriptstyle{(1)}}\Delta\Theta)^{2}+(\gamma+\gamma_{0}^{\scriptscriptstyle{(0)}}+\gamma_{0}^{\scriptscriptstyle{(1)}}\Delta\Theta)^{2}}=$
$\displaystyle 2\pi
R_{0}\tilde{\alpha}\sigma[(\Theta_{0}+\Delta\Theta)^{4}-\Theta_{0}^{4}].$ (27)

The above formula constitutes a non-linear equation which has to be solved for
$\Delta\Theta$.

By solving Eq. (IV) numerically in the frequency domain we find that a
bistabily window of width $\Delta\lambda=1.6\rm{nm}$ exists at

$\frac{2\pi R_{0}\tilde{\alpha}\sigma}{I_{0}L}=6.15\cdot 10^{-12}$ (28)

at the incident wavelength $\lambda=795.9~{}\rm{nm}$. This in combination with
our assessment for the emissivity $\tilde{\alpha}=2.80\cdot 10^{-5}$ yields
the incident intensity $I_{0}=0.54~{}\rm{W/m^{2}}$. Our findings are
illustrated in Fig. 6. Figure 6 (a), and Fig. 6 (b) demonstrate thermo-optical
hysteresis in the frequency domain. One can see that at the peak of the
resonance the system is heated up to $\Theta=187~{}^{\circ}\rm{C}$. Figure 6
(c), and Fig. 6 (d) demonstrate thermo-optical hysteresis in the intensity
domain with the width of bistability window $\Delta I_{0}=1.6~{}\rm{W/m^{2}}$.

##  V Summary and Conclusions

We considered thermo-optic hysteresis in a system supporting a symmetry
protected bound state in the continuum. It is found that bistability occurs at
surprisingly low intensities of incident light, about
$I_{0}=0.54~{}\rm{W/m^{2}}$. Two aspects are of key importance for obtaining
low bistability thresholds. (i) Critical coupling which is obtained by
breaking the symmetry of the system by incident light. With the increase of
the angle of incidence the radiation decay rate grows quadratically until it
becomes equal to the absorption decay rate, $\gamma=\gamma_{0}$. In this
situation, no matter how small the absorption decay rate is, $50\%$ of the
incident light is absorbed and transformed to heat. This approach paves a way
for engineering effective absorbers with low loss dielectrics. (ii) Low
emissivity which is due to a small optical path within the absorbing
dielectric at the peak wavelength of Planck’s distribution at about
$10~{}\rm{\mu m}$. The dielectric array is assembled of rods with a
subwavelength cross section in the visible. Thus, the optical path in the mid-
infrared is at least an order of magnitude less than the wavelength of the
thermal radiation. Given that the imaginary part of the refractive index of
silicon is vanishingly small in the mid-infrared Chandler-Horowitz and
Amirtharaj (2005), $n^{\prime\prime}\approx 10^{-4}$, this results in a small
value of the emissivity $\tilde{\alpha}=2.8\cdot 10^{-5}$.

The bistability occurs with intensities eleven orders of magnitude smaller
than those needed for hysteresis due to Kerr nonlinearity in silicon Maksimov
_et al._ (2020). Thus, that thermo-optic effects should dominate in non-linear
response of silicon nano-structures supporting bound states in the continuum.
In this work we have used a rough estimate for the emissivity. Yet, it is
worth mentioning that even with $\tilde{\alpha}=1$, i.e. with a mid-infrared
black body, the bistability thresholds would drop by less than five orders of
magnitude and still be dominating in the non-linear response. Here we
neglected diffusive and convective heat transfer. In particular, the system
under consideration has no substrate, which is necessary for realistic
dielectric metasurfaces. One would expect that the presence of substrate would
further decrease the bistability thresholds as the heat would be leaking into
the substrate. This issue will be addressed in future studies.

## Acknowledgments

This work received financial support through the grant of Russian Science
Foundation and Krasnoyarsk Regional Fund of Science No 22-22-20056,
https://rscf.ru/project/22-22-20056/, the authors acknowledge discussions with
A.A. Bogdanov.

## References

  * Hsu  _et al._ (2016) C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, Nature Reviews Materials 1, 16048 (2016).
  * Koshelev  _et al._ (2019a) K. Koshelev, G. Favraud, A. Bogdanov, Y. Kivshar, and A. Fratalocchi, Nanophotonics 8, 725 (2019a).
  * Sadreev (2021) A. F. Sadreev, Reports on Progress in Physics (2021).
  * Kim  _et al._ (1999) C. S. Kim, A. M. Satanin, Y. S. Joe, and R. M. Cosby, Physical Review B 60, 10962 (1999).
  * Shipman and Venakides (2005) S. P. Shipman and S. Venakides, Physical Review E 71, 026611 (2005).
  * Sadreev  _et al._ (2006) A. F. Sadreev, E. N. Bulgakov, and I. Rotter, Physical Review B 73, 235342 (2006).
  * Blanchard  _et al._ (2016) C. Blanchard, J.-P. Hugonin, and C. Sauvan, Physical Review B 94, 155303 (2016).
  * Yuan and Lu (2017) L. Yuan and Y. Y. Lu, Physical Review A 95, 023834 (2017).
  * Yoon  _et al._ (2015) J. W. Yoon, S. H. Song, and R. Magnusson, Scientific Reports 5, 18301 (2015).
  * Mocella and Romano (2015) V. Mocella and S. Romano, Physical Review B 92, 155117 (2015).
  * Bulgakov  _et al._ (2011) E. Bulgakov, K. Pichugin, and A. Sadreev, Physical Review B 83, 045109 (2011).
  * Bulgakov  _et al._ (2013) E. Bulgakov, K. Pichugin, and A. Sadreev, Journal of Physics: Condensed Matter 25, 395304 (2013).
  * Yuan and Lu (2016) L. Yuan and Y. Y. Lu, Physical Review A 94, 013852 (2016).
  * Bulgakov and Sadreev (2014) E. Bulgakov and A. Sadreev, Optics letters 39, 5212 (2014).
  * Bulgakov and Sadreev (2010) E. N. Bulgakov and A. F. Sadreev, Physical Review B 81, 115128 (2010).
  * Ndangali and Shabanov (2013) F. R. Ndangali and S. V. Shabanov, in  _Active Photonic Materials V_ , Vol. 8808 (International Society for Optics and Photonics, 2013) p. 88081F.
  * Koshelev  _et al._ (2019b) K. Koshelev, A. Bogdanov, and Y. Kivshar, Science Bulletin 64, 836 (2019b).
  * Wang and Zhang (2018) T. Wang and S. Zhang, Optics express 26, 322 (2018).
  * Rybin  _et al._ (2017) M. V. Rybin, K. L. Koshelev, Z. F. Sadrieva, K. B. Samusev, A. A. Bogdanov, M. F. Limonov, and Y. S. Kivshar, Physical review letters 119, 243901 (2017).
  * Bogdanov  _et al._ (2019) A. A. Bogdanov, K. L. Koshelev, P. V. Kapitanova, M. V. Rybin, S. A. Gladyshev, Z. F. Sadrieva, K. B. Samusev, Y. S. Kivshar, and M. F. Limonov, Advanced Photonics 1, 016001 (2019).
  * Carletti  _et al._ (2018) L. Carletti, K. Koshelev, C. De Angelis, and Y. Kivshar, Physical review letters 121, 033903 (2018).
  * Koshelev  _et al._ (2020) K. Koshelev, S. Kruk, E. Melik-Gaykazyan, J.-H. Choi, A. Bogdanov, H.-G. Park, and Y. Kivshar, Science 367, 288 (2020).
  * Yuan and Lu (2020) L. Yuan and Y. Y. Lu, SIAM Journal on Applied Mathematics 80, 864 (2020).
  * Krasikov  _et al._ (2018) S. D. Krasikov, A. A. Bogdanov, and I. V. Iorsh, Physical Review B 97, 224309 (2018).
  * Bulgakov and Maksimov (2019) E. N. Bulgakov and D. N. Maksimov, Scientific Reports 9, 7153 (2019).
  * Maksimov  _et al._ (2020) D. N. Maksimov, A. A. Bogdanov, and E. N. Bulgakov, Physical Review A 102, 033511 (2020).
  * Hu  _et al._ (2020) Z. Hu, L. Yuan, and Y. Y. Lu, Physical Review A 101, 013806 (2020).
  * Saadabad  _et al._ (2021) R. M. Saadabad, L. Huang, and A. E. Miroshnichenko, Physical Review B 104, 235405 (2021).
  * Tan  _et al._ (2022) L. Tan, L. Yuan, and Y. Y. Lu, Journal of the Optical Society of America B 39, 611 (2022).
  * Zhang and Zhang (2015) M. Zhang and X. Zhang, Scientific reports 5, 1 (2015).
  * Wang  _et al._ (2020) X. Wang, J. Duan, W. Chen, C. Zhou, T. Liu, and S. Xiao, Physical Review B 102, 155432 (2020).
  * Sang  _et al._ (2021) T. Sang, S. A. Dereshgi, W. Hadibrata, I. Tanriover, and K. Aydin, Nanomaterials 11, 484 (2021).
  * Xiao  _et al._ (2021) S. Xiao, X. Wang, J. Duan, T. Liu, and T. Yu, JOSA B 38, 1325 (2021).
  * Cai  _et al._ (2022) Y. Cai, X. Liu, K. Zhu, H. Wu, and Y. Huang, Journal of Quantitative Spectroscopy and Radiative Transfer 283, 108150 (2022).
  * Zograf  _et al._ (2021) G. P. Zograf, M. I. Petrov, S. V. Makarov, and Y. S. Kivshar, Advances in Optics and Photonics 13, 643 (2021).
  * Sun and Reano (2010) P. Sun and R. M. Reano, Optics letters 35, 1124 (2010).
  * Khandekar and Rodriguez (2017) C. Khandekar and A. W. Rodriguez, Applied Physics Letters 111, 083104 (2017).
  * Gao  _et al._ (2017) Y. Gao, W. Zhou, X. Sun, H. K. Tsang, and C. Shu, Optics Letters 42, 1950 (2017).
  * Pottier and Bellon (2021) B. Pottier and L. Bellon, SciPost Physics 10, 120 (2021).
  * Ryabov  _et al._ (2022) D. Ryabov, O. Pashina, G. Zograf, S. Makarov, and M. Petrov, Nanophotonics 11, 3981 (2022).
  * Fan  _et al._ (2003) S. Fan, W. Suh, and J. D. Joannopoulos, Journal of the Optical Society of America A 20, 569 (2003).
  * Bulgakov and Sadreev (2015) E. N. Bulgakov and A. F. Sadreev, Physical Review A 92, 023816 (2015).
  * Bulgakov and Maksimov (2017a) E. N. Bulgakov and D. N. Maksimov, Physical Review A 96, 063833 (2017a).
  * Jellison and Modine (1994) G. E. Jellison and F. A. Modine, Journal of Applied Physics 76, 3758 (1994).
  * Bulgakov and Maksimov (2017b) E. N. Bulgakov and D. N. Maksimov, Optics Express 25, 14134 (2017b).
  * Vainstein (1988) L. A. Vainstein,  _Electromagnetic waves_ (Moscow Izdatel Radio Sviaz, 1988).
  * Schinke  _et al._ (2015) C. Schinke, P. C. Peest, J. Schmidt, R. Brendel, K. Bothe, M. R. Vogt, I. Kröger, S. Winter, A. Schirmacher, S. Lim, H. T. Nguyen, and D. MacDonald, AIP Advances 5, 067168 (2015).
  * Chandler-Horowitz and Amirtharaj (2005) D. Chandler-Horowitz and P. M. Amirtharaj, Journal of Applied Physics 97, 123526 (2005).
  * Bergman  _et al._ (2011) T. L. Bergman, F. P. Incropera, D. P. Dewitt, and A. S. Lavine,  _Fundamentals of heat and mass transfer_ (John Wiley & Sons, 2011).
  * Wada and Kamijoh (1996) H. Wada and T. Kamijoh, Japanese Journal of Applied Physics 35, L648 (1996).