# Constrained Visual-Inertial Localization With Application And Benchmark in
Laparoscopic Surgery

Regine Hartwig1,3, Daniel Ostler1, Jean-Claude Rosenthal2
Hubertus Feußner1, Dirk Wilhelm1, Dirk Wollherr3 1Authors are with Research
Group MITI, Surgical Department, Technical University of Munich (TUM),
Germany2Authors are with Fraunhofer HHI, Berlin, Germany3Authors are with
Chair of Automatic Control Engineering, Department of Electrical and Computer
Engineering, Technical University of Munich (TUM), Germany
<EMAIL_ADDRESS>
###### Abstract

We propose a novel method to tackle the visual-inertial localization problem
for constrained camera movements. We use residuals from the different
modalities to jointly optimize a global cost function. The residuals emerge
from IMU measurements, stereoscopic feature points, and constraints on
possible solutions in $\operatorname{SE}(3)$. In settings where dynamic
disturbances are frequent, the residuals reduce the complexity of the problem
and make localization feasible. We verify the advantages of our method in a
suitable medical use case and produce a dataset capturing a minimally invasive
surgery in the abdomen. Our novel clinical dataset MITI is comparable to
state-of-the-art evaluation datasets, contains calibration and synchronization
and is available at [1].

##  I INTRODUCTION

Many researchers have studied simultaneous localization and mapping (SLAM) in
the computer vision, and robotics community [2], inspiring many innovative
products. Applications such as virtual reality, augmented reality [3],
unmanned vehicles, and autonomous robots [4] rely heavily on stable pose
estimates based on SLAM systems.

Challenging situations include a dynamic environment which prevents the SLAM
algorithms from being applicable in many use-cases. The works of [5, 6, 7]
show that it is possible to add residuals specially designed for a use-case
and jointly optimize a cost function. Motion priors, reprojection from a
calibrated multi-camera setup, or additional measurements reduce the degrees
of freedom on the solution space.

MIS has proven to have many medical advantages for the patient. Stable
virtual/ augmented reality and assisting roboter arms could lead to enhanced
hand-eye coordination in MIS and subsequently to a whole new way of conducting
surgeries. Furthermore, SLAM is a game-changer for applications such as
mapping detected metastasizing carcinoma for cancer staging, safety routines
(e.g., collision avoidance), autonomous systems, image-guided surgery, input
generation for artificial intelligence, and anatomic scene graphs. The huge
disadvantages (magnetic distortion, line-of-sight, room consumption)
disqualify available tracking equipment for most interventions. Data from this
domain was acquired, and implementation shows the applicability of the
proposed constrained visual-inertial SLAM method in this challenging dynamic
surrounding.

The paper is structured as follows: Section IV presents the visual-inertial
odometry algorithm with movement constraints. Section V discusses our novel
dataset [1] and provides useful information for users. In Section VI we will
discuss the implementation of the proposed algorithm, discuss properties of
the calibrated data from our dataset, and evaluate the performance by
experiments.

##  II PROBLEM STATEMENT AND CONTRIBUTION

#### Constrained Visual-Inertial Localization

[3], and further developments have become popular state-of-the-art methods
[8]. These algorithms focus on feature point matching based on binary
descriptors, even for small viewpoint changes, where descriptor matching is
made real-time applicable by speeding up the process with priorly trained or
online learned [9] bag-of-words vocabularies [10, 11].

However, we focus on domains where feature point descriptor matching and
tracking perform worse due to challenging image properties [12]. Sparse or
semi-dense optical flow [13, 5, 14] are alternatives that perform well forhigh
video framerate and can speed up trajectory generation for the tracking part
of the SLAM system.

Additional information from inertial measurement units (IMU) [15, 16, 17]
improves visual odometry methods in difficult situations, where tracking is
lost, and does not need an external measurement equipment, i.e., avoiding
line-of-sight problems.

In contrast to all previously mentioned strategies, we additionally exploit
naturally given motion constraints [19], which is also famous in the medical
domain as remote center-of-motion (RCM) [20][21], reducing the solution space
to tackle the more challenging image domains, which contain large movements in
the image, occlusions, fog, deformable surfaces or difficult illuminations
conditions [21].

Subsequently, we outline our approach in more detail: We formulate residuals
for IMU, stereo camera, and movement constraints and show that the solution
space reduces by pivoting the camera around a point from 6 to 4 degrees of
freedom. Afterward, we propose to use optical flow in combination with
gyroscope readings as an odometry alternative to descriptor-based feature
tracking for dynamic image data and high video frequency. Finally, for the
refinement of this odometry, we formulate the global cost function and
optimization strategies. Our method’s main advantage and contribution is the
ability to relocalize the camera in unseen environments by minimizing the
global cost function, which incorporates all residual terms emerging from
sensors and constraints. The gravity and magnetic field provide the necessary
absolute references. In addition, the Lucas-Kanade tracking of the stereo
reconstructed locally rigid parts in the environment provides information
about the distance from the pivot point.

#### Dataset

Publicly available datasets play a vital role in the evolution of new
algorithms, and technologies [22, 23, 24, 25]. Researchers focusing on the
minimally invasive application have already created endoscopic datasets [27,
21, 25], which are suitable for investigating computer tomography
registration, tissue deformation/ reconstruction [28, 29, 30], and instrument
segmentation. Since they are not complete or are not capturing a handheld
intervention, those are unsuitable for testing SLAM applications. In this
paper, we introduce a novel dataset suitable for testing SLAM in the MIS
domain. To the best of our knowledge, our dataset is the first publicly
available surgical dataset comparable to [22, 23, 24]. Our MITI Dataset [1],
available under the Creative Commons Attribution license (CC-BY 4.0), provides
all necessary data by a complete recording of a handheld surgical intervention
at Research Hospital Rechts der Isar of TUM. It contains multimodal sensor
information from IMU, stereoscopic video, and infrared (IR) tracking as ground
truth for evaluation. Furthermore, calibration for the stereoscope,
accelerometer, magnetometer, the rigid transformations in the sensor setup,
and time-offsets are available. We wisely chose a suitable intervention,
namely diagnostic laparoscopy, that contains very few cutting and tissue
deformation and shows a full scan of the abdomen with a handheld camera such
that it is ideal for testing SLAM algorithms. It incorporates 667 seconds with
images (1080p60 RGB) for each camera, readings of two IMU sensors at 220Hz,
and IR transformations to two different targets at 20Hz.

Intending to promote the progress of visual-inertial algorithms designed for
MIS application, we hope that our clinical training dataset helps and enables
researchers to enhance algorithms.

##  III NOTATION

The following notions are adapted from [31, 32, 33] using the rotation groups
$\operatorname{SO}(3)$ and the rigid motion group $\operatorname{SE}(3)$,
which fulfill the properties of a Lie group. We use the notation
$\exp:\mathfrak{m}\to\mathcal{M}$ to denote the mapping from tangent space,
i.e. Lie algebra $\mathfrak{m}$ to Lie group $\mathcal{M}$. The Lie algebra
$\mathfrak{so}(3){=}\operatorname{skew}(3)$ of $\operatorname{SO}(3)$ has
$d{=}3$ degrees of freedom and $\mathfrak{se}(3)$, tangent space of
$\operatorname{SE}(3)$, has $d{=}6$ degrees of freedom. They are isometric
isomorph to $\mathbb{R}^{d}$, and we use the function $\operatorname{Exp}$ as
the composition of the hat operator and the matrix exponential $\exp$. The
function can explicitly be given in closed form by the Rodriguez formula. The
inverse mapping is denoted by $\operatorname{Log}$, which is the mapping from
the Lie-Group to $\mathbb{R}^{d}$. We further use the
$\ominus:\mathcal{M}\times\mathcal{M}\to\mathbb{R}^{d}$ as

$\mathbf{T}_{1}\ominus\mathbf{T}_{2}:=\operatorname{Log}\left(\mathbf{T}_{1}^{-1}\mathbf{T}_{2}\right),$
(1)

which computes the difference between two transformation matrices mapped to
the tangent spaces’ isomorph space $\mathbb{R}^{d}$.

##  IV VISUAL-INERTIAL SENSOR FUSION

Our method splits into two parts: Firstly, estimating the pose starting from
previous keyframes by movement constraints, gyroscope readings, and Lucas-
Kanade tracks of stereoscopic matches. Secondly, the algorithm conducts a
windowed optimization over gathered keyframes, changing their location by
jointly minimizing the residuals from all modalities, eliminating any drift
that might accumulate.

###  IV-A Residuals

This subsection describes the residuals gathered from different modalities.
The residuals are functions of the localization pose
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}\in\operatorname{SE}(3)$
with the property that they are zero at the actual position and orientation.

#### Motion Constraints

The transformation
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}\in\operatorname{SE}(3)$
transforms a vector in homogeneous coordinates from T to C by multiplication
from the right. Since the motion of C is constrained to move around pivot
point T but can also move forward and backward in the x-direction, we can
reduce the 6 degrees of freedom to 4 degrees of freedom

$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}=\begin{bmatrix}\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{R}&\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{t}\\\
\boldsymbol{0}^{T}&1\end{bmatrix}=\begin{bmatrix}I&\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{t}\\\
\boldsymbol{0}^{T}&1\end{bmatrix}\begin{bmatrix}\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{R}&\boldsymbol{0}\\\
\boldsymbol{0}^{T}&1\end{bmatrix},$ (2)

with
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{t}=t_{x}\boldsymbol{e}_{x}$.
This means a point in T is first rotated in the pivot point and then
translated along the x-axis. The corresponding residual that enforces the
result to have only 4 degrees of freedom is

$\boldsymbol{r}_{pivot}=t_{y}\boldsymbol{e}_{y}+t_{z}\boldsymbol{e}_{z}.$ (3)

This adds a penalty to the minimization problem if the solution is not in the
reduced space.

#### Accelerometer

The accelerometer readings constrain the current camera pose by one
corresponding bearing vector in world coordinates. We assume constant velocity
as also proposed by [34]. Rotating the gravity
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{g}$
bearing vector from T to C and comparing it with the measurement
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{g}$
gives the residual

$\boldsymbol{r}_{a}=\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{g}-\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{R}\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{g}.$
(4)

#### Magnetometer

The magnetometer provides similar to the accelerometer a correspondence in
world coordinates

$\boldsymbol{r}_{m}=\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{m}-\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{R}\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{m},$
(5)

with the magnetic field
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{m}$
in T and measurement
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{m}$
in C.

#### Global Mappoints

The distance between a mappoint
${}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{x}$
in global coordinates and a matched 2D position $\boldsymbol{u}$ in image
plane of current camera coordinates is

$\boldsymbol{r}_{2D}=\Pi(\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}\,\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{x})-\boldsymbol{u}$
(6)

with $\Pi(\cdot)$ being the projection into distorted image plane and
$\boldsymbol{u}{:=}\boldsymbol{u}_{d}$ the distorted image coordinates. In
(6), the map points
${}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{x}$
are connected rigidly to the world coordinate system T at initialization by
using the current camera pose estimate. We need to rewrite (6), attaching the
mappoints to the camera coordinate system. Otherwise, since we optimize only
the poses
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}$,
the minimization of other residuals leads to the increase of (6). For this
sake, we store the landmark positions relative to the camera coordinate system
C, in which the stereo matcher initialized them. This adds another computation
step and one more parameter
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\mathbf{T}_{t^{\prime}}$
at detection time $t^{\prime}$ to the residual in (6), namely

$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{x}=\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\mathbf{T}_{t^{\prime}}\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{x}_{t^{\prime}}.$
(7)

Thus we have

$\boldsymbol{r}_{2D}=\Pi(\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\mathbf{T}_{t^{\prime}}\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{x}_{t^{\prime}})-\boldsymbol{u},$
(8)

which connects the observation frame with the initialization frame of a
landmark.

#### Gyroscope

The gyroscope measurements in tangent space of
$\boldsymbol{\omega}_{t^{\prime}}\in\mathbb{R}^{3}\equiv\mathfrak{so}(3)$ are
integrated by gyroscope sample time $\Delta t$ and the non-abelsh product

$\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\Delta\mathbf{R}_{t_{1},t_{2}}:=\prod_{t^{\prime}\in\\{t_{1},...,t_{2}-\Delta
t\\}}\operatorname{Exp}(\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{\omega}_{t^{\prime}}\Delta
t).$ (9)

This gives the update rule

$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\mathbf{R}_{t_{2}}=\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\mathbf{R}_{t_{1}}\,\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\Delta\mathbf{R}_{t_{1},t_{2}}$
(10)

and the residual term

$\boldsymbol{r}_{g}=\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\mathbf{R}_{t_{2}}\ominus\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\mathbf{R}_{t_{1}}\,\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\Delta\mathbf{R}_{t_{1},t_{2}}.$
(11)

which constraints consecutive transformations to have a certain deviation in
orientation, namely the integrated measurements of the angular velocity.

###  IV-B Visual-Inertial Odometry Tracking

This subsection describes our approach for the odometry tracking thread. It
updates from the pose of the previous keyframe, generating good
initializations for the optimization using gyroscope data and Lucas Kanade
tracks of stereo correspondences. Algorithm 1 sketches the procedure, and in
the following, we will explain line by line the computation steps the
algorithm conducts.

After entering the loop, it first updates the pose based on the gyroscope
measurements in Line 2. Although (11) usually only provides information about
the orientation, a position update is possible due to the camera movement
having only 4 degrees of freedom

$\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{t_{2}}=\begin{bmatrix}I&\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{t}_{t_{1}}\\\\[4.0pt]
\boldsymbol{0}^{T}&1\end{bmatrix}\begin{bmatrix}\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\Delta\mathbf{R}_{t_{1},t_{2}}^{T}\,\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{R}_{t_{1}}&\boldsymbol{0}\\\\[4.0pt]
\boldsymbol{0}^{T}&1\end{bmatrix}$ (12)

with
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\boldsymbol{t}_{t}=t_{x,t}\boldsymbol{e}_{x}$.
The new pose is on a sphere with radius $t_{x,t}$, around the trocar entry
point, leaving only the radius unknown but estimating all other degrees of
freedom.

During the first iteration, the Lines 3 to 4 have no effect. The first frame
is a keyframe, and the algorithm computes stereo correspondences by descriptor
matching using ORB descriptors [35] in Line 6. ORB descriptors are faster to
compute and compare than SIFT or SURF, whose benefits, e.g., scale and
rotation invariance, are less essential for stereo-matching. For faster
matching, it uses the projections provided by the camera calibration
parameters as initialization for finding matches in the second image [8]. The
epipolar error threshold [36] is adjusted to filter outliers

$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C0}}}\boldsymbol{x}^{T}\mathbf{E}\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C1}}}\boldsymbol{x}>d_{th},$
(13)

with epipolar matrix $\mathbf{E}=\hat{\boldsymbol{t}}\mathbf{R}$ and
$\hat{\cdot}:\mathbb{R}^{3}\to\text{skew}(3)$, defining the transformation
between left C0 and right C1 lens. The matched points initialize new landmarks
in Line 7.

Afterwards, in Line 3 for every incoming frame it updates the landmark
projections $\boldsymbol{u}\in\mathbb{R}^{2}$ in the 2D image plane, producing
residuals of the form (6). For that, it uses the variational Lucas-Kanade-
Method [37] taking two consecutive frames for every camera. This variational
method performs well (small baseline, photometric consistency) if video
frequency is high. For large displacements, the drift becomes large, and for
fast movements, we lose track. Nevertheless, optical flow is suitable for
local tracking.

With initialization computed in Line 2 it determines the camera pose
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{t}$
based on the 3D-2D-correspondences in Line 4 by a perspective-n-point
algorithm [38] inside a random sample consensus (RANSAC) [39] scheme with
refinement provided by OpenGV [40]. A sanity check on the rotation part using
the IMU measurements and on the translation part using the constraints ensures
locally rigid features.

Algorithm 1 Visual-Inertial Odometry with Optimization

0: Initialization
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{t_{0}}$,
Images $\\{I_{0,t_{i}},I_{1,t_{i}}\\}_{i\in I}$; IMU Data
$\\{\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\Delta\mathbf{R}_{t_{i},t_{i+1}},\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{a}_{t_{i}},\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{m}_{t_{i}}\\}_{i\in
I}$

0:
$\\{\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{t_{i}}\\}_{i\in
I}$

1: while True do

2: Update
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{t_{i-1}}$
based on
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\Delta\mathbf{R}_{t_{i-1}}$
and movement constraints

3: Lukas-Kanade tracking of landmarks

4: Update
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{t_{i}}$
based on 2D-landmark projections

5: if Decision for keyframe then

6: Stereo matching

7: Landmark insertion

8: end if

9: $i=i+1$

10: end while

###  IV-C Optimization

The optimization presented in this subsection jointly minimizes all previously
formulated residuals over a window of keyframes. As shown later the camera can
globally relocalize itself by solving the proposed optimization problem solely
by small tracks of locally rigid feature points and IMU measurements. In order
to use all the information from sensor readings and images we build a cost
function by summing up the normalized, weighted and robustified residuals that
have been collected, i.e. (3, 4, 5, 8, 11) which are available at each
timestep. We denote them as residuals $\boldsymbol{r}_{k}$, $k\in K$, where
$k$ indicates time and residual type.

The overall optimization problem is

$\operatorname*{arg\,min}_{\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{n:m}}\quad\sum_{k}\rho_{k}\left(\left\|\alpha_{k}\boldsymbol{r}_{k}\left(\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{n:m}\right)\right\|\right){:=}\sum_{k}f_{k}(\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{n:m}),$
(14)

optimizing over the window of keyframe localizations
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{n:m}:=\left(\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{t_{n}},...,\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{t_{m}}\right)$.
In above formulation, $\alpha_{k}$ stands for the scaling factor, and $\rho$
for the robustification function. Each residual gets an individual scaling
factor $\alpha_{k}$. The evaluation part of this paper gives a more detailed
overview of the chosen weight factors, optimization window size, and trigger
of the optimization problem.

For robustification of residuals that are prone to outliers, we choose the
Huber loss

$\rho_{k}(x)=\begin{cases}\frac{1}{2}{x^{2}}&\text{for }|x|\leq\delta,\\\
\delta(|x|-\frac{1}{2}\delta),&\text{otherwise.}\end{cases}$ (15)

The function input $x$ in our application is the normalized and weighted
residual $\|\alpha_{k}r_{k}\|$. For residuals with a dense error distribution
and no outliers, $\rho_{k}$ is simply the identity.

Each of the residuals
$\boldsymbol{r}_{k}\left(\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{n:m}\right)$
depends only on a maximum of two parameters
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}\in\operatorname{SE}(3)$,
and further, many parameters do not have a residual in common. The residual
terms give the optimization problem the typical graph structure of SLAM
problems. This sparse dependence on parameters makes the Hesse-matrix sparse,
and the optimization problem suitable for a sparse Schur solver [41].

##  V DATASET

In this section, we introduce our dataset and calibration procedures. In Fig.
1 we show the sensor setup during the data acquisition process. The sensor
data consists firstly of a stereoscopic video (Tipcam1, Karl Storz) from the
view inside the abdomen. Secondly, two 9DOF IMU (MetamotionR, Metawear) sensor
boards attached to the laparoscope handle sending: angular velocity,
acceleration, and magnetic field strength. Thirdly, an IR tracking camera
(Polaris Vega, NDI) detects two targets, each consisting of 3-4 reflecting
passive spheres. The anticipated precision of the IR tracking is 0.25mm [42],
making it suitable as ground truth to evaluate the SLAM algorithm’s
performance. We calibrated the camera and the IMU, determining sensor-specific
parameters and calculated time offsets using the state-of-the-art methods
defined in subsequent subsections. Parameters can be found in the path
calibration/*.csv of the dataset [1].

(a)

(b)

Figure 1:  Data Acquisition (a) IMU sensorboards and IR passive targets
attached to stereoscopic laparocopic camera (b) IR sensor and 3D monitor in
the operating room

###  V-A Camera Calibration

The pinhole projection [43] $\Pi_{0}:\mathbb{R}^{3}\to\mathbb{R}^{2}$
$\boldsymbol{x}\mapsto\boldsymbol{u}_{u}$ from a point
$\boldsymbol{x}\in\mathbb{R}^{3}$ in coordinate system of left C0 or right C1
lens to the undistorted image point $\boldsymbol{u}_{u}$ in the image plane is
parameterized by focal lengths $f_{x},f_{y}$ and optical center $c_{x},c_{y}$
both measured in $\mathrm{pixel}$.

Radial image distortion is compensated by a low order polynomial model
$d:\mathbb{R}^{2}\to\mathbb{R}^{2}$,
$\boldsymbol{u}_{u}\mapsto\boldsymbol{u}_{d}$. Additionally, the extrinsic
parameters are calibrated, which determine the rigid body transformation
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C0}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C1}}}T\in\operatorname{SE}(3)$.
The function $\Pi(\cdot)$ in residual (6) is a composition of rotation and
translation from C to C0 or C1 and projection to distorted 2D image
coordinates by composition of $\Pi_{0}$ and $d$. (LABEL:calib_cam_radial).

###  V-B IMU Calibration

Using a standard approach [34, 44], the accelerometer measurements
$\tilde{\boldsymbol{a}}(t)$ are modeled as

$\tilde{\boldsymbol{a}}(t)=(\boldsymbol{a}(t)-\boldsymbol{g}(t))d_{a}^{-1}+\boldsymbol{n}_{a}(t)+\boldsymbol{b}_{a}$
(16)

with acceleration $\boldsymbol{a}(t)\in\mathbb{R}^{3}$, gravity
$\boldsymbol{g}(t)\in\mathbb{R}^{3}$, errors
$\boldsymbol{b}_{a}\in\mathbb{R}^{3}$, $d_{a}\in\mathbb{R}$ and
$\boldsymbol{n}_{a}(t)\sim\mathcal{N}(0,\sigma_{a}^{2}\mathbf{I})$.
Calibration assumes constant velocity and gravity norm
$g=\|\boldsymbol{g}(t)\|=9.81\frac{m}{s^{2}}$. Furthermore, we model the error
of magnetometer readings $\tilde{\boldsymbol{m}}(t)$ [34, 44] by

$\tilde{\boldsymbol{m}}(t)=\boldsymbol{m}(t)d_{m}^{-1}+\boldsymbol{n}_{m}(t)+\boldsymbol{b}_{m},$
(17)

with magnetic field strength $\boldsymbol{m}(t)\in\mathbb{R}^{3}$, errors
$\boldsymbol{b}_{m}\in\mathbb{R}^{3}$, $d_{m}\in\mathbb{R}$ and
$\boldsymbol{n}_{m}(t)\sim\mathcal{N}(0,\sigma_{m}^{2}\mathbf{I})$. The
corrected measurements lie on a sphere with radius equal to the strength of
the earths magnetic field $m=\|\boldsymbol{m}\|=48.6\mu T$ in Munich, Germany.
This information is used in the calibration procedure. Fig. 2 shows calibrated
sensor readings.

(a)

(b)

Figure 2: Calibrated IMU measurements, i.e., red in front, blue behind sphere
(a)
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{a}-\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{g}$
in $9.81\frac{m}{s^{2}}$ and unit sphere (b)
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{m}$
in $\mu T$ and sphere with radius $48.6$

###  V-C Time Offset Calibration

The sensor data has been synchronized during the acquisition process by
configuring a network time protocol server/client at the data receiving
devices.

Figure 3: Time-calibrated angular velocity part
$\|\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{\omega}^{(k)}_{t_{i}}\|$

For synchronization of the video the displayed timestamp has been captured at
the start and end of the intervention. The sensor interfaces for IR/IMU are
ethernet/bluetooth respectively, which causes a remaining timeoffset in the
collected data. This offset is minimized by evaluating the tangent
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{\xi}_{t_{i}}\in\mathbb{R}^{6}\cong\mathfrak{se}(3)$
for each sensor at coordinate system C at time $t_{i}$ We then calibrated the
timeoffset $dt$ by comparing the velocities acquired from different sensors

$j,k\in\\{\text{"Visual Odometry", "IR", "IMU0", "IMU1"}\\}.$

The sum of differences at different timeoffsets is

$\min_{dt}\sum_{i}\left\|\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{\xi}^{(j)}_{t_{i}+dt}-\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{\xi}^{(k)}_{t_{i}}\right\|.$
(18)

For that, we determine the camera-based angular velocity using the visual
odometry algorithm based on Lucas Kanade feature tracks without IMU and IR.
The computed timeoffsets are 100ms/166ms for the IR/IMU readings. Time
synchronization is vital for short-term odometry, unlike windowed
optimization, where keyframes are more distant.

##  VI EVALUATION AND EXPERIMENTS

We conduct experiments on an implementation of the approach described in
Section IV. In Subsection VI-C we examine the influence of weight distribution
on residuals, and in Subsection VI-B we discuss the advantages of utilizing
robust loss functions. Subsection VI-A and VI-D are devoted to the evaluation
of our method with our MITI dataset. In doing that, we measure the deviation
of the algorithms output to ground-truth poses obtained by IR tracking.

We used several libraries: OpenCV for keypoint detection, feature tracking,
and reprojection; OpenGV [40] for triangulation, RANSAC, and nonlinear
refinement; as well as Ceres [45] and Sophus [46] for solving the optimization
problem. For the experiment, we took the image sequence starting from
timestamp 637273ms, which contains a fast movement from left to right with
dynamic occlusions and enables us to demonstrate the algorithm’s performance
in a scenario where standard SLAM approaches lose tracking and can not
relocalize in the unseen environment.

###  VI-A Evaluation of Visual Odometry Tracking

In Fig. 4 one can inspect the debug view of the output and Fig. 4(b) and Fig.
4(c) show huge deviations for different settings in the algorithm.

(a)

(b)

(c)

Figure 4: Debug environment (a) Stereo view with ORB-descriptor-based
stereomatches in green, Lucas-Kanade-tracks in teal, (b) and (c) Camera poses
connected to lines, that should go through pivot point and 3D map points of
left and right diaphragm, (b) Camera-only based tracking, (c) Incorporating
IMU sensordata and movement constraint

(a)

(b)

Figure 5: Algorithm 1 without optimization (a) V1 relies on IMU data and can
reduce rotational error in the long run, whereas V2 based on vision suffers
from drift, (b) error of V1 is unacceptable when movements do not stay on the
sphere, error of pure visual odometry V2 also suffers from drift

The isolated bundle of lines in Fig. 4(b) shows: without movement constraints
and IMU measurements, the optimization is unable to converge to the correct
poses since it loses tracking of 2D features in between keyframes. We computed
the Log distance between IR tracked and constrained visual-inertial odometry
(VO) estimated camera poses

$\boldsymbol{\xi}=[\boldsymbol{\omega},\boldsymbol{\upsilon}]^{T}=\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}^{(\text{IR})}\ominus\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}^{(\text{VO})}\in\mathbb{R}^{6}$
(19)

for variants of the computation steps in Algorithm 1, that are listed in Table
I. V1 integrates IMU sensor-readings in Line 2 into pose estimates as given in
(12). Therefore, it cannot detect changes in the x-translation component,
resulting in a large translational error as soon as the camera moves inside.
V2 only uses the camera-based pose estimation of Lines 3 and 4 but without
taking the initialization provided by the IMU in Line 2. The result of this
variant can additionally be inspected in the debug view of Fig. 4(b). V3 uses
both camera-based and IMU-based tracking, which already leads to a significant
decrease in both error residuals.

Table I: Settings of Alg. 1 for results shown in Fig. 5 Variant | Pivot point online | Line2 | Line 3 and 4 | Optimization
---|---|---|---|---
V1 | $\times$ | ✓ | $\times$ | $\times$
V2 | $\times$ | $\times$ | ✓ | $\times$
V3 | $\times$ | ✓ | ✓ | $\times$
V4 | ✓ | ✓ | ✓ | $\times$

Finally, in V4, instead of considering the pivot point given by the priorly
computed position from IR tracking, we initialize and iteratively refine the
pose purely based on visual odometry. The result is an essential step for the
goal an eternal inside-out tracking.

###  VI-B Error Analysis and Choice of Robustification Functions

The residuals themselves are composed of errors in sensor readings and the
deviation of the current pose estimate. Aiming to minimize only the part
indicating a pose deviation, we examine errors of sensor readings and their
propagation to the residual term to decide whether to use a Huber loss. The
Huber loss for the reprojection residual (8) is justified, since we identified
an outlier characteristic in the three error sources: A drift error in the
Lucas-Kanade tracked 2D locations
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{u}$
, that can be modeled as a random walk, after $n$ tracking steps, contributes
with a Gaussian probability distribution of
$\boldsymbol{\epsilon}\sim\mathcal{N}(\boldsymbol{0},n\sigma^{2}\mathbf{I})$
linearly to the residual. The second error source is a badly reconstructed
landmark position
$\,{}^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}\boldsymbol{x}_{t^{\prime}}$,
which occurs although the epipolar error threshold reduces the probability of
outliers. The third error arises from a failing filtering of dynamic feature
points. For residuals that penalize solutions not going through the trocar
entry point (3) the use of the Huber loss is not justified, since there is no
measurement error, and we aim to minimize large residuals. The error of IMU
measurements (4, 5, 9) is approximately Gaussian, which leads us to not
choosing the Huber loss.

###  VI-C Weighting Residuals

In this subsection we discuss the factor $\alpha_{k}$ that strongly influences
the optimization (14). Firstly, a task of the parameter is to adjust the
residuals, making them independent on their measurement unit. We ran Algorithm
1 and built a statistic for every residual type, as shown in Tab II.

Table II: Statistic of Residuals | Eq. | Expectation $\operatorname{E}[\boldsymbol{r}]$ | Variance $\operatorname{Var}[\boldsymbol{r}]$ | Unit
---|---|---|---|---
$\boldsymbol{r}_{2D}$ | (6) | 1.6263 | 9.3767 | $\mathrm{Pixel}$
$\boldsymbol{r}_{g}$ | (9) | 0.0003 | 0.0051 | $\mathrm{rad}$
$\boldsymbol{r}_{pivot}$ | (3) | 0.0511 | 0.4374 | $\mathrm{cm}$
$\boldsymbol{r}_{a}$ | (4) | 0.0235 | 0.0405 | $9.81\frac{\mathrm{m}}{\mathrm{s}^{2}}$
$\boldsymbol{r}_{m}$ | (5) | -0.2730 | 4.2801 | $\mu\mathrm{T}$

The expectation of the residuals is close to zero, which simplifies the
normalization
$\frac{\boldsymbol{r}-\operatorname{E}[\boldsymbol{r}_{k}]}{\sqrt{\operatorname{Var}[\boldsymbol{r}_{k}]}}$.
The scaling factor splits into
$\alpha_{k}=\frac{1}{N_{k}}\frac{1}{\sqrt{\operatorname{Var}[\boldsymbol{r}_{k}]}}\beta_{k}$,
with $\beta_{k}$ being an additional term to adjust the importance between
residuals, as discussed in the next Subsection.

(a)

(b)

(c)

Figure 6: Errors and cost reduction of V3 with optimization for different
values of tradeoff value $\gamma$

Hereby, we additionally normalize by the number of residuals of the same type
$N_{k}$ e.g., there exist $N_{k}=1$ gyroscope readings per keyframe, but in
general $N_{k}>1$ reprojection residuals per keyframe.

###  VI-D Evaluation of Trade-off and Optimization

A trade-off between gyroscope and camera-based residuals is adjusted by
$\beta_{k}=\gamma$ for residual $k$ belonging to gyroscope measurements and
$\beta_{k}=(1-\gamma)$ for reprojection residuals, with $\gamma\in[0,1]$. We
conduct experimental runs with optimization being triggered every 10 keyframes
and a window size of 10. Fig. 6(a) and 6(b) show the deviation of the SLAM
algorithm from the ground-truth IR tracked poses for different values of
$\gamma$. Since we want to jointly optimize over all residuals, the cost
emerging from different sensors should simultaneously decrease. The cost
function difference $\Delta f$ is

$\Delta
f=f(\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{n:m}^{(\text{after})})-f(\,^{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{C}}}_{\text{\tiny\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}}\mathbf{T}_{n:m}^{(\text{before})}),$
(20)

with evaluations before and after optimization. Negative values indicate that
this $k$-th part of the cost function was minimized. In Fig. 6(c) we show the
average cost difference for different residual types in one optimization step
for different values of $\gamma$. It exhibits a simultaneous cost decrease for
values close to $\gamma=0,4$.

On the other hand, we observe by inspecting the error trajectories in Fig.
6(a) and 6(b) that for values close to $0.5$ the optimization successfully
incorporates the advantages of both sensors: The rotational error is low at
the beginning, which is typical for camera-based tracking, but also stays
down, which is possible due to utilizing the information of IMU sensors. The
distribution of weights shows great potential in influencing the overall
outcome of the algorithm. An optimal adjustment however, is a dynamic weight
distribution, depending on the current conditions that influence the
reliability of the sensor data.

##  VII DISCUSSION

Movement constraints and additional sensor information are necessary to
provide stable pose estimates in dynamic environments. Absolute references are
provided by the earth’s gravity and magnetic field, whereas gyroscope and
Lucas-Kanade tracks provide the odometry. For a more robust setup, however,
descriptor matching of rigid features can improve the results of long-term
SLAM.

##  VIII CONCLUSIONS

This work shows that visual-inertial algorithms can be applied to dynamically
challenging fields such as minimally invasive surgery. This sets a new
paradigm in the applicability of SLAM. The gain of reliability is due to the
fact that we use inherent motion constraints. Our results show for the use-
case of minimally invasive surgery, that self localization is indeed feasible
and hope that researchers can profit from our publicly available dataset.

## References

  * [1] R. Hartwig, D. Ostler, J.-C. Rosenthal, H. Feussner, D. Wilhelm, and D. Wollherr, “Miti: Slam benchmark for laparoscopic surgery,” TUM, 2021. [Online]. Available: https://mediatum.ub.tum.de/1621941
  * [2] C. Cadena, L. Carlone, H. Carrillo, Y. Latif, D. Scaramuzza, J. Neira, I. Reid, and J. J. Leonard, “Past, present, and future of simultaneous localization and mapping: Toward the robust-perception age,” _Trans. on robotics_ , vol. 32, no. 6, pp. 1309–1332, 2016.
  * [3] G. Klein and D. Murray, “Parallel tracking and mapping for small ar workspaces,” in _Int. symposium on mixed and augmented reality_ , 2007, pp. 225–234.
  * [4] J. Janai, F. Güney, A. Behl, A. Geiger _et al._ , “Computer vision for autonomous vehicles: Problems, datasets and state of the art,” _Foundations and Trends in Computer Graphics and Vision_ , vol. 12, no. 1–3, pp. 1–308, 2020.
  * [5] C. Forster, Z. Zhang, M. Gassner, M. Werlberger, and D. Scaramuzza, “Svo: Semidirect visual odometry for monocular and multicamera systems,” _Trans. on Robotics_ , vol. 33, no. 2, pp. 249–265, 2016.
  * [6] C. Kerl, J. Sturm, and D. Cremers, “Robust odometry estimation for rgb-d cameras,” in _Int. Conf. on Robotics and Automation (ICRA)_ , 2013.
  * [7] S. Leutenegger, S. Lynen, M. Bosse, R. Siegwart, and P. Furgale, “Keyframe-based visual–inertial odometry using nonlinear optimization,” _The International Journal of Robotics Research_ , vol. 34, no. 3, pp. 314–334, 2015.
  * [8] R. Mur-Artal and J. D. Tardós, “Orb-slam2: An open-source slam system for monocular, stereo, and rgb-d cameras,” _Trans. on Robotics_ , vol. 33, no. 5, pp. 1255–1262, 2017.
  * [9] S. Khan and D. Wollherr, “Ibuild: Incremental bag of binary words for appearance based loop closure detection,” in _Int. Conf. on Robotics and Automation (ICRA)_ , 2015, pp. 5441–5447.
  * [10] D. Gálvez-López and J. D. Tardos, “Bags of binary words for fast place recognition in image sequences,” _Trans. on Robotics_ , vol. 28, no. 5, pp. 1188–1197, 2012.
  * [11] D. Nister and H. Stewenius, “Scalable recognition with a vocabulary tree,” in _Int. Conf. on Computer Vision and Pattern Recognition_ , vol. 2. IEEE, 2006, pp. 2161–2168.
  * [12] L. Maier-Hein, A. Groch, A. Bartoli, S. Bodenstedt, G. Boissonnat, P.-L. Chang, N. T. Clancy, D. S. Elson, S. Haase, E. Heim, J. Hornegger, P. Jannin, H. Kenngott, T. Kilgus, B. Müller-Stich, D. Oladokun, S. Röhl, T. R. dos Santos, H.-P. Schlemmer, A. Seitel, S. Speidel, M. Wagner, and D. Stoyanov, “Comparative validation of single-shot optical techniques for laparoscopic 3-d surface reconstruction,” _Trans. on Medical Imaging_ , vol. 33, no. 10, pp. 1913–1930, 2014.
  * [13] Q. Fu, H. Yu, X. Wang, Z. Yang, H. Zhang, and A. Mian, “Fastorb-slam: A fast orb-slam method with coarse-to-fine descriptor independent keypoint matching,” _arXiv preprint arXiv:2008.09870_ , 2020.
  * [14] J. Engel, V. Koltun, and D. Cremers, “Direct sparse odometry,” _Trans. on pattern analysis and machine intelligence_ , vol. 40, no. 3, pp. 611–625, 2017\.
  * [15] L. Von Stumberg, V. Usenko, and D. Cremers, “Direct sparse visual-inertial odometry using dynamic marginalization,” in _Int. Conf. on Robotics and Automation (ICRA)_ , 2018, pp. 2510–2517.
  * [16] T. Qin, P. Li, and S. Shen, “Vins-mono: A robust and versatile monocular visual-inertial state estimator,” _Trans. on Robotics_ , vol. 34, no. 4, pp. 1004–1020, 2018.
  * [17] T. Qin, J. Pan, S. Cao, and S. Shen, “A general optimization-based framework for local odometry estimation with multiple sensors,” _arXiv preprint arXiv:1901.03638_ , 2019.
  * [18] C. Mei, G. Sibley, M. Cummins, P. Newman, and I. Reid, “Rslam: A system for large-scale mapping in constant-time using stereo,” _Int. journal of computer vision_ , vol. 94, no. 2, pp. 198–214, 2011.
  * [19] M. Salerno, K. Zhang, A. Menciassi, and J. S. Dai, “A novel 4-dofs origami enabled, sma actuated, robotic end-effector for minimally invasive surgery,” in _Int. Conf. on Robotics and Automation (ICRA)_ , 2014, pp. 2844–2849.
  * [20] F. Vasconcelos, E. Mazomenos, J. Kelly, and D. Stoyanov, “Rcm-slam: Visual localisation and mapping under remote centre of motion constraints,” in _2019 International conference on robotics and automation (ICRA)_. IEEE, 2019, pp. 9278–9284.
  * [21] B. Lin, Y. Sun, X. Qian, D. Goldgof, R. Gitlin, and Y. You, “Video-based 3d reconstruction, laparoscope localization and deformation recovery for abdominal minimally invasive surgery: a survey,” _The International Journal of Medical Robotics and Computer Assisted Surgery_ , vol. 12, no. 2, pp. 158–178, 2016.
  * [22] S. Baker, D. Scharstein, J. Lewis, S. Roth, M. J. Black, and R. Szeliski, “A database and evaluation methodology for optical flow,” _Int. journal of computer vision_ , vol. 92, no. 1, pp. 1–31, 2011, available at https://vision.middlebury.edu/flow/. [Online]. Available: https://vision.middlebury.edu/flow/
  * [23] A. Geiger, P. Lenz, and R. Urtasun, “Are we ready for autonomous driving? the kitti vision benchmark suite,” in _Int. Conf. on computer vision and pattern recognition_ , 2012, pp. 3354–3361.
  * [24] J. Sturm, N. Engelhard, F. Endres, W. Burgard, and D. Cremers, “A benchmark for the evaluation of rgb-d slam systems,” in _Int. Conf. on Intelligent Robots and Systems_ , 2012, pp. 573–580.
  * [25] M. Allan, J. Mcleod, C. Wang, J. C. Rosenthal, Z. Hu, N. Gard, P. Eisert, K. X. Fu, T. Zeffiro, W. Xia _et al._ , “Stereo correspondence and reconstruction of endoscopic data challenge,” _arXiv preprint arXiv:2101.01133_ , 2021.
  * [26] M. D. Blikkendaal, S. R. Driessen, S. P. Rodrigues, J. P. Rhemrev, M. J. Smeets, J. Dankelman, J. J. van Den Dobbelsteen, and F. W. Jansen, “Surgical flow disturbances in dedicated minimally invasive surgery suites: an observational study to assess its supposed superiority over conventional suites,” _Surgical endoscopy_ , vol. 31, no. 1, pp. 288–298, 2017.
  * [27] D. Dey, D. Gobbi, P. Slomka, K. Surry, and T. Peters, “Automatic fusion of freehand endoscopic brain images to three-dimensional surfaces: creating stereoscopic panoramas,” _Trans. on Medical Imaging_ , vol. 21, no. 1, pp. 23–30, 2002.
  * [28] N. Mahmoud, T. Collins, A. Hostettler, L. Soler, C. Doignon, and J. M. M. Montiel, “Live tracking and dense reconstruction for handheld monocular endoscopy,” _Trans. on Medical Imaging_ , vol. 38, no. 1, pp. 79–89, 2019\.
  * [29] P. Mountney, D. Stoyanov, and G.-Z. Yang, “Three-dimensional tissue deformation recovery and tracking,” _Signal Processing Magazine_ , vol. 27, no. 4, pp. 14–24, 2010.
  * [30] L. Maier-Hein, P. Mountney, A. Bartoli, H. Elhawary, D. Elson, A. Groch, A. Kolb, M. Rodrigues, J. Sorger, S. Speidel _et al._ , “Optical techniques for 3d surface reconstruction in computer-assisted laparoscopic surgery,” _Medical image analysis_ , vol. 17, no. 8, pp. 974–996, 2013.
  * [31] T. D. Barfoot, _State estimation for robotics_. Cambridge University Press, 2017.
  * [32] J. Sola, J. Deray, and D. Atchuthan, “A micro lie theory for state estimation in robotics,” _arXiv preprint arXiv:1812.01537_ , vol. 5, p. 1371, 2018.
  * [33] C. Sommer, V. Usenko, D. Schubert, N. Demmel, and D. Cremers, “Efficient derivative computation for cumulative b-splines on lie groups,” in _Int. Conf. on Computer Vision and Pattern Recognition_ , 2020, pp. 11 148–11 156.
  * [34] S. Madgwick, “An efficient orientation filter for inertial and inertial/magnetic sensor arrays,” _Report x-io and University of Bristol (UK)_ , vol. 25, pp. 113–118, 2010.
  * [35] E. Rublee, V. Rabaud, K. Konolige, and G. Bradski, “Orb: An efficient alternative to sift or surf,” in _Int. Conf. on computer vision_. IEEE, 2011, pp. 2564–2571.
  * [36] D. Forsyth and J. Ponce, _Computer vision: A modern approach._ Prentice hall, 2011.
  * [37] A. Bruhn, J. Weickert, and C. Schnörr, “Lucas/kanade meets horn/schunck: Combining local and global optic flow methods,” _Int. journal of computer vision_ , vol. 61, no. 3, pp. 211–231, 2005.
  * [38] V. Lepetit, F. Moreno-Noguer, and P. Fua, “Epnp: An accurate o (n) solution to the pnp problem,” _International journal of computer vision_ , vol. 81, no. 2, pp. 155–166, 2009.
  * [39] M. A. Fischler and R. C. Bolles, “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography,” _Communications of the ACM_ , vol. 24, no. 6, pp. 381–395, 1981.
  * [40] L. Kneip and P. Furgale, “Opengv: A unified and generalized approach to real-time calibrated geometric vision,” in _Int. Conf. on Robotics and Automation (ICRA)_. IEEE, 2014, pp. 1–8.
  * [41] G. Grisetti, R. Kümmerle, H. Strasdat, and K. Konolige, “g2o: A general framework for (hyper) graph optimization,” in _Int. Conf. on Robotics and Automation (ICRA)_ , 2011, pp. 9–13. [Online]. Available: http://ais.informatik.uni-freiburg.de/publications/papers/kuemmerle11icra.pdf
  * [42] G. Fattori, A. J. Lomax, D. C. Weber, and S. Safai, “Technical assessment of the ndi polaris vega optical tracking system,” _Radiation oncology_ , vol. 16, no. 1, pp. 1–4, 2021.
  * [43] Z. Zhang, “A flexible new technique for camera calibration,” _Trans. on pattern analysis and machine intelligence_ , vol. 22, no. 11, pp. 1330–1334, 2000\.
  * [44] V. Renaudin, M. H. Afzal, and G. Lachapelle, “Complete triaxis magnetometer calibration in the magnetic domain,” _Journal of sensors_ , vol. 2010, pp. 1–10, 2010.
  * [45] S. Agarwal and K. Mierle, “Ceres solver: Tutorial & reference,” _Google Inc_ , vol. 2, no. 72, p. 8, 2012.
  * [46] H. Strasdat, S. Lovegrove, and Others, “Sophus library,” https://github.com/strasdat/Sophus.


#  Balmer Decrement Anomalies in Galaxies at $z\sim 6$ Found by JWST
Observations:
Density-Bounded Nebulae or Excited Hi Clouds?

Hiroto Yanagisawa  Institute for Cosmic Ray Research, The University of Tokyo,
5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Department of Physics,
Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo,
Tokyo 113-0033, Japan  Masami Ouchi National Astronomical Observatory of
Japan, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  Institute for Cosmic Ray Research, The University of Tokyo,
5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Department of Astronomical
Science, SOKENDAI (The Graduate University for Advanced Studies), Osawa
2-21-1, Mitaka, Tokyo, 181-8588, Japan  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba
277-8583, Japan  Kimihiko Nakajima National Astronomical Observatory of Japan,
National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588,
Japan  Hidenobu Yajima Center for Computational Sciences, University of
Tsukuba, Ten-nodai, 1-1-1 Tsukuba, Ibaraki 305-8577, Japan  Hiroya Umeda
Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha,
Kashiwa, Chiba 277-8582, Japan  Department of Physics, Graduate School of
Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan
Shunsuke Baba Institute of Space and Astronautical Science (ISAS), Japan
Aerospace Exploration Agency (JAXA), 3-1-1 Yoshinodai, Chuo-ku, Sagamihara,
Kanagawa 252-5210, Japan  Takao Nakagawa Institute of Space and Astronautical
Science (ISAS), Japan Aerospace Exploration Agency (JAXA), 3-1-1 Yoshinodai,
Chuo-ku, Sagamihara, Kanagawa 252-5210, Japan  Minami Nakane Institute for
Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa,
Chiba 277-8582, Japan  Department of Physics, Graduate School of Science, The
University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan  Akinori
Matsumoto  Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Department of Physics, Graduate
School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo
113-0033, Japan  Yoshiaki Ono Institute for Cosmic Ray Research, The
University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Yuichi
Harikane Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Yuki Isobe Institute for Cosmic
Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba
277-8582, Japan  Department of Physics, Graduate School of Science, The
University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan  Yi Xu
Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha,
Kashiwa, Chiba 277-8582, Japan  Department of Astronomy, Graduate School of
Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan
Yechi Zhang National Astronomical Observatory of Japan, National Institutes of
Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

###### Abstract

We investigate the physical origins of the Balmer decrement anomalies in GS-
NDG-9422 (Cameron et al. 2023) and RXCJ2248-ID (Topping et al. 2024) galaxies
at $z\sim 6$ whose $\mathrm{H}\alpha/\mathrm{H}\beta$ values are significantly
smaller than $2.7$, the latter of which also shows anomalous
$\mathrm{H}\gamma/\mathrm{H}\beta$ and $\mathrm{H}\delta/\mathrm{H}\beta$
values beyond the errors. Because the anomalous Balmer decrements are not
reproduced under the Case B recombination, we explore the nebulae with the
optical depths smaller and larger than the Case B recombination by physical
modeling. We find two cases quantitatively explaining the anomalies; 1)
density-bounded nebulae that are opaque only up to around Ly$\gamma$-Ly8
transitions and 2) ionization-bounded nebulae partly/fully surrounded by
optically-thick excited Hi clouds. The case of 1) produces more H$\beta$
photons via Ly$\gamma$ absorption in the nebulae, requiring fine tuning in
optical depth values, while this case helps ionizing photon escape for cosmic
reionization. The case of 2) needs the optically-thick excited Hi clouds with
$N_{2}\simeq 10^{12}-10^{13}$ $\mathrm{cm^{-2}}$, where $N_{2}$ is the column
density of the hydrogen atom with the principal quantum number of $n=2$.
Interestingly, the high $N_{2}$ values qualitatively agree with the recent
claims for GS-NDG-9422 with the strong nebular continuum requiring a number of
$2s$-state electrons and for RXCJ2248-ID with the dense ionized regions likely
coexisting with the optically-thick clouds. While the physical origin of the
optically-thick excited Hi clouds is unclear, these results may suggest gas
clouds with excessive collisional excitation caused by an amount of accretion
and supernovae in the high-$z$ galaxies.

High-redshift galaxies (734), Nebulae (1095), Interstellar medium (847),
Photoionization (2060)

††software: Astropy (Astropy Collaboration et al., 2022), Cloudy (Chatzikos et
al., 2023), matplotlib (Hunter, 2007), NumPy (Harris et al., 2020), PyNeb
(Luridiana et al., 2015), SciPy (Virtanen et al., 2020), and YMCMC (Hsyu et
al., 2020).

##  1 Introduction

Balmer decrements (i.e., deviations of hydrogen Balmer line ratios from
theoretical values) are known as tracers of dust extinction in galaxies.
However, recent observations with James Webb Space Telescope (JWST) have
identified anomalous Balmer decrements in high-$z$ galaxies. Cameron et al.
(2023) have reported that GS-NDG-9422 at $z=5.943$ shows
$\mathrm{H}\alpha/\mathrm{H}\beta$ significantly smaller than the theoretical
limit with no dust extinction at the $>1\sigma$ level. Similarly, Topping et
al. (2024) have found that RXCJ2248-ID at $z=6.11$ has the anomalously small
$\mathrm{H}\alpha/\mathrm{H}\beta$ value together with the
$\mathrm{H}\gamma/\mathrm{H}\beta$ and $\mathrm{H}\delta/\mathrm{H}\beta$
values significantly larger than the theoretical limits with no dust
extinction. These Balmer decrements cannot be explained by the dust
extinction, since the dust extinction raises
$\mathrm{H}\alpha/\mathrm{H}\beta$ and reduces
$\mathrm{H}\gamma/\mathrm{H}\beta$ and $\mathrm{H}\delta/\mathrm{H}\beta$.
Although these studies do not discuss the physical origins of these anomalies
quantitatively, there may be an important physical mechanism that changes
Balmer line ratios beyond the well-known theoretical limits that are
calculated under the Case B recombination corresponding to the optically-thick
Lyman series. To understand the Balmer decrement anomalies, our study aims to
quantitatively explore the physical origins of the Balmer decrement anomalies,
allowing the optical depths smaller and larger than the Case B recombination.

The structure of this paper is as follows. In Section 2, we present data used
in this paper. In Section 3, we explore the possible physical origins of the
Balmer decrement anomalies. In Section 4, we present some implications of our
results. Section 5 summarizes our results.

##  2 Data

In this paper, we investigate two galaxies, GS-NDG-9422 and RXCJ2248-ID, at
$z\sim 6$ with the Balmer decrement anomalies reported in the previous
studies. We summarize the Balmer line ratios of the GS-NDG-9422 and
RXCJ2248-ID in Table 1.

Table 1: Balmer line ratios of our sample

| GS-NDG-9422 | RXCJ2248-ID
---|---|---
$\mathrm{H}\alpha/\mathrm{H}\beta$ | $2.65\pm 0.03$ | $2.547\pm 0.022^{\dagger}$
$\mathrm{H}\gamma/\mathrm{H}\beta$ | $0.48\pm 0.02$ | $0.512\pm 0.008$
$\mathrm{H}\delta/\mathrm{H}\beta$ | $0.25\pm 0.03$ | $0.283\pm 0.012$

Note. — ($\dagger$) The reported flux refers to the narrow component (M.
Topping, private communication).

###  2.1 GS-NDG-9422

GS-NDG-9422 was observed with JWST NIRSpec via the JWST Advanced Deep
Extragalactic Survey (JADES; PID: 1210, PI: Luetzgendorf; Bunker et al. 2023;
Eisenstein et al. 2023). Cameron et al. (2023) have analysed the spectra and
measured fluxes of the rest-UV and optical emission lines. Cameron et al.
(2023) have reported the Balmer line ratios of
$\mathrm{H}\alpha/\mathrm{H}\beta=2.65\pm 0.03$,
$\mathrm{H}\gamma/\mathrm{H}\beta=0.48\pm 0.02$, and
$\mathrm{H}\delta/\mathrm{H}\beta=0.25\pm 0.03$. Because all of these flux
ratios are calculated from the Prism data covering all of these Balmer
emission lines through the one slit position, slit-loss effect is not a cause
of the Balmer decrement anomalies. None of these Balmer emission lines show
multiple velocity components including a broad component.

###  2.2 RXCJ2248-ID

RXCJ2248-ID was observed with JWST NIRSpec in General Observers (GO) program
(PID: 2478, PI: Stark). Topping et al. (2024) have investigated the spectra
and presented the Balmer line ratios of
$\mathrm{H}\alpha/\mathrm{H}\beta=2.547\pm 0.022$,
$\mathrm{H}\gamma/\mathrm{H}\beta=0.512\pm 0.008$, and
$\mathrm{H}\delta/\mathrm{H}\beta=0.283\pm 0.012$. Because the
$\mathrm{H}\alpha$ and $\mathrm{H}\beta$ emission lines are covered by the
single G395M spectrum, the difference of the slit-loss effect is not
problematic in the $\mathrm{H}\alpha/\mathrm{H}\beta$ ratio. Although the
$\mathrm{H}\gamma$ and $\mathrm{H}\delta$ are covered by both G235M and G395M
gratings, whether the reported fluxes are taken from G235M or G395M is not
clarified in the literature. The $\mathrm{H}\alpha$ line shows a broad
component with the broad-to-narrow line flux ratio of 0.22, while the other
hydrogen lines do not. The $\mathrm{H}\alpha/\mathrm{H}\beta$ value listed in
Table 1 refers to the narrow-line $\mathrm{H}\alpha$ fluxes (M. Topping,
private communication).

Figure 1: Comparisons of the observed and Case B Balmer line ratios. The
orange and red squares represent the values of GS-NDG-9422 (Cameron et al.,
2023) and RXCJ2248-ID (Topping et al., 2024). The gray shaded regions denote
the Case B values, where we change the electron temperature from
$15,000\,\mathrm{K}$ to $25,000\,\mathrm{K}$ and the electron density from
$10^{1}\,\mathrm{cm^{-3}}$ to $10^{5}\,\mathrm{cm^{-3}}$. In the right panel,
the Case B values of a $3\times 3$ grid of $T_{\mathrm{e}}$ and
$n_{\mathrm{e}}$ are shown. The cyan, blue, and magenta points represent the
values of $T_{\mathrm{e}}=15,000\,\mathrm{K}$, $20,000\,\mathrm{K}$, and
$25,000\,\mathrm{K}$, respectively. The circles, triangles, and inverted
triangles show the values of $n_{\mathrm{e}}=10^{1}\,\mathrm{cm^{-3}}$,
$10^{3}\,\mathrm{cm^{-3}}$, and $10^{5}\,\mathrm{cm^{-3}}$, respectively. The
black arrows represent the directions of the effect of the dust extinction and
the Balmer decrement anomalies. Table 2: Parameters used in Cloudy

Parameter | Value | Description
---|---|---
Blackbody | 20,000, 40,000 | Incident radiation [K]
$\log U$ | $-2,-1$ | Ionization parameter
$\log n$ | 2, 5 | Density of the hydrogen gas [$\mathrm{cm^{-3}}$]
Stop efrac(†) | $-1$ | Stopping criterion
Stop thickness(††) | $10^{10-20}$ | Stopping criterion [$\mathrm{cm}$]
Geometry | Plane-parallel | Geometry of the gas
Abundance | $0.1$ | Chemical composition of the gas [$Z_{\odot}$]

Note. — ($\dagger$) When log of the electron fraction is below this value, the
calculation stops. ($\dagger\dagger$) When the thickness of the nebula reaches
this value, the calculation stops.

##  3 Explanation for the Balmer Decrement Anomalies

Here we explore the origin of the Balmer decrement anomalies presented in
Section 2. In Section 3.1, we compare the anomalies with the values in the
Case B recombination, where we define the Case B as optically-thick in the
Lyman series and optically-thin in the other hydrogen lines. We then
investigate the two scenarios with optical depths smaller and larger than the
one of the Case B recombination that are optically-thin in the Lyman series
(Section 3.2) and optically-thick in the Balmer series (Section 3.3),
respectively.

###  3.1 Comparison with the Case B Recombination

Figure 1 compares the observed and Case B values of Balmer line ratios. The
Case B values are calculated with PyNeb (version 1.1.18, Luridiana et al.
2015) for electron temperatures $T_{\mathrm{e}}=15,000-25,000\,\mathrm{K}$ and
electron densities $n_{\mathrm{e}}=10^{1}-10^{5}\,\mathrm{cm^{-3}}$. GS-
NDG-9422 presents $\mathrm{H}\delta/\mathrm{H}\beta$ and
$\mathrm{H}\gamma/\mathrm{H}\beta$ consistent with the Case B values, while
the $\mathrm{H}\alpha/\mathrm{H}\beta$ is tentatively lower than the Case B
value at the $>1\sigma$ level. RXCJ2248-ID shows
$\mathrm{H}\delta/\mathrm{H}\beta$ and $\mathrm{H}\gamma/\mathrm{H}\beta$
larger than those of Case B at the $\sim 1.5\sigma$ and $\sim 4\sigma$ levels,
respectively, whereas the $\mathrm{H}\alpha/\mathrm{H}\beta$ is smaller than
the Case B value at the $\sim 8\sigma$ level. These deviations are opposite
from the effect of the dust extinction. We therefore cannot explain these
Balmer line ratios by the combination of Case B and dust extinction.

Figure 2: From top to bottom, $\mathrm{H}\alpha/\mathrm{H}\beta$,
$\mathrm{H}\gamma/\mathrm{H}\beta$, $\mathrm{H}\delta/\mathrm{H}\beta$, and
$\tau_{\mathrm{Ly}\beta}$ as a function of thickness of the nebula. The left
and right panels correspond to the blackbody temperatures of
$20,000\,\mathrm{K}$ and $40,000\,\mathrm{K}$, respectively. The dased and
solid lines represent the values calculated at $n_{\mathrm{e}}=10^{2}$ and
$10^{5}\,\mathrm{cm^{-3}}$, respectively. The gray and black lines indicate
the values calculated at $\log U=-2$ and $-1$, respectively.
$\tau_{\mathrm{Ly}\gamma}$, $\tau_{\mathrm{Ly}\delta}$, and
$\tau_{\mathrm{Ly}8}$ are lower than $\tau_{\mathrm{Ly}\beta}$ by amounts
corresponding to the length of the arrows in the bottom-left panel. The orange
(red) line and shaded region denote the observed Balmer line ratios and their
uncertainties of GS-NDG-9422 (RXCJ2248-ID), respectively. The horizontal lines
in the bottom panels denote $\tau_{\lambda}=1$, above which the line becomes
optically thick.

###  3.2  Optical Depth smaller than the Case B Recombination

Here we investigate the scenario with optically-thin Lyman series lines, which
have smaller optical depths than those assumed in the Case B recombination.
This physical condition is realized in the density-bounded nebulae. We
simulate the density-bounded nebulae using version 22.02 of Cloudy (Chatzikos
et al., 2023). We vary the stopping thickness of the calculation from
$10^{10}$ to $10^{20}\,\mathrm{cm}$. We change hydrogen densities from
$10^{2}$ to $10^{5}\,\mathrm{cm^{-3}}$, which cover the electron densities of
$\sim 10^{3}$ and $\sim 10^{5}\,\mathrm{cm^{-3}}$ for GS-NDG-9422 and
RXCJ2248-ID reported by Cameron et al. (2023) and Topping et al. (2024),
respectively. We also vary ionization parameter $U$ from $\log U=-2$ to $-1$,
motivated by the high ionization parameters in GS-NDG-9422 and RXCJ2248-ID
implied by Cameron et al. (2023) and Topping et al. (2024), respectively. The
other parameters are summarized in Table 2. The calculated Balmer line ratios
and optical depths of Lyman series as a function of thickness of the nebula
are shown in Figure 2. Figure 2 indicates that the
$\mathrm{H}\alpha/\mathrm{H}\beta$ ratios follow N-shaped changes as optical
depths of Lyman series increase. This may be explained as follows (For detail,
see Figure 4.3 of Osterbrock & Ferland 2006). First, as the thickness
increases, the optical depth of Ly$\beta$ ($\tau_{\mathrm{Ly}\beta}$;
hereafter $\tau_{\lambda}$ refers to an optical depth of a line $\lambda$)
approaches to $1$. This leads $\mathrm{Ly}\beta$ to be optically thick and a
conversion of $\mathrm{Ly}\beta$ to $\mathrm{H}\alpha$ occurs. When
$\tau_{\mathrm{Ly}\gamma}$ approaches to 1, a conversion of
$\mathrm{Ly}\gamma$ to $\mathrm{H}\beta$ occurs, which lowers the
$\mathrm{H}\alpha/\mathrm{H}\beta$ ratio. Conversions of higher order Lyman
lines take place similarly for large thickness. When the optical depth around
the one of Ly8 (transition from a principal quantum number of 8 to 1) is
larger than 1, $\mathrm{H}\alpha/\mathrm{H}\beta$ gradually approaches to the
Case B value. For $\mathrm{H}\gamma/\mathrm{H}\beta$, the first effect of
increasing the thickness is the conversion of Ly$\gamma$ to $\mathrm{H}\beta$
due to $\tau_{\mathrm{Ly}\gamma}$ approaching to 1, which lowers the
$\mathrm{H}\gamma/\mathrm{H}\beta$ ratio. As $\tau_{\mathrm{Ly}\delta}$
approaches to 1, a conversion of Ly$\delta$ to $\mathrm{H}\gamma$ occurs,
which leads to the increase of $\mathrm{H}\gamma/\mathrm{H}\beta$. Further
increase of the thickness makes the higher order Lyman lines optically thick.
When $\tau_{\mathrm{Ly}8}$ is above 1, $\mathrm{H}\gamma/\mathrm{H}\beta$
approaches to the Case B value. The $\mathrm{H}\delta/\mathrm{H}\beta$ ratios
can also be explained similarly. For example, in the case of the blackbody
temperature of $20,000\,\mathrm{K}$, $\log n=2$, and $\log U=-1$ (the gray
solid lines in the left panels), the Balmer line ratios of GS-NDG-9422 and
RXCJ2248-ID are reproduced when the thickness is $\sim 10^{13.4}$ and $\sim
10^{13.8}\,\mathrm{cm}$, respectively, where $\tau_{\mathrm{Ly}8}$ is $\sim
1$. Fine-tuning of the thickness is thus required to reproduce the Balmer line
ratios of these galaxies. For the other $U$ and $n$, the observed Balmer line
ratios are reproduced at different thicknesses, while the necessity of the
fine-tuning does not change. In the case of the blackbody temperature of
$40,000\,\mathrm{K}$, the values of Balmer line ratios are different from
those of $20,000\,\mathrm{K}$ because of the temperature dependence of line
emissivities. The Balmer line ratios of GS-NDG-9422 are reproduced within the
wide ranges of thickness. However, for RXCJ2248-ID,
$\mathrm{H}\alpha/\mathrm{H}\beta$ and $\mathrm{H}\gamma/\mathrm{H}\beta$ are
not reproduced within the 1$\sigma$ level at any of the thickness, ionization
parameter and density, while $\mathrm{H}\delta/\mathrm{H}\beta$ are reproduced
for $\tau_{\mathrm{Ly}\beta}<-1$ or $\tau_{\mathrm{Ly}8}\sim 0.5$. We further
investigate Balmer line ratios for other ranges of parameters (blackbody
temperatures from $20,000$ to $50,000\,\mathrm{K}$, ionization parameters from
$\log U=-2$ to $0$, and densities from $\log n=2$ to $5$). The calculated
Balmer line ratios show the similar trends to those of the parameter sets
discussed above.

Since the Balmer line ratios observed in GS-NDG-9422 and RXCJ2248-ID are
determined with the small errors, only very narrow ranges of thickness are
allowed to explain the Balmer line ratios. Fine-tuning of the thickness is
thus necessary at given $U$, $n$, and blackbody temperature, and such
situations are realized with a very low probability.

###  3.3  Optical Depth larger than the Case B Recombination

We investigate the scenario with optically-thick Balmer series lines, which
have larger optical depths than those assumed in the Case B recombination. We
first investigate a simple optically-thick ionized nebula similar to those of
Yano et al. (2021), who worked on nearby ultraluminous infrared galaxies for
the Brackett lines, by Cloudy modeling similar to those in Section 3.2.
However, we find no solutions explaining the Balmer decrement anomalies of GS-
NDG-9422 and RXCJ2248-ID with any physical parameters in optically-thick
ionized nebula.

Instead, we explore the possibility that the ionization-bounded nebulae are
observed through the thick excited neutral hydrogen gas. In other words,
ionized nebulae are partly (or fully) surrounded by thick excited H i gas.

####  3.3.1 Analytic Calculations

The optical depth of a transition
$\mathcal{N}^{\prime}\rightarrow\mathcal{N}$, where $\mathcal{N}^{\prime}$ and
$\mathcal{N}$ are principal quantum numbers of a hydrogen atom that satisfy
$\mathcal{N}^{\prime}>\mathcal{N}$, is given by

$\tau_{\mathcal{N}^{\prime},\mathcal{N}}=\alpha_{\mathcal{N},\mathcal{N}^{\prime}}N_{\mathcal{N}},$
(1)

where $\alpha_{\mathcal{N},\mathcal{N}^{\prime}}$ is an absorption cross
section of a transition $\mathcal{N}\rightarrow\mathcal{N}^{\prime}$ and
$N_{\mathcal{N}}$ is a column density of hydrogen atoms at the principal
quantum number $\mathcal{N}$ 111In this section, we follow the mathematical
notations of Yano et al. (2021). Assuming that the line velocity profile is a
Gaussian function, $\alpha_{\mathcal{N},\mathcal{N}^{\prime}}$ is given by

$\alpha_{\mathcal{N},\mathcal{N}^{\prime}}=\frac{hc}{4\pi^{3/2}}\frac{B_{\mathcal{N},\mathcal{N}^{\prime}}}{v_{\mathrm{Dop}}},$
(2)

where $h$ is the Planck constant, $c$ is the speed of light,
$B_{\mathcal{N},\mathcal{N}^{\prime}}$ is a Einstein B coefficient of the
transition $\mathcal{N}\rightarrow\mathcal{N}^{\prime}$, and
$v_{\mathrm{Dop}}$ is a Doppler velocity half width. If the velocity of the
gas is determined by thermal motion, the Doppler velocity half width is
expressed as $v_{\mathrm{Dop}}=\sqrt{2kT/m_{\mathrm{H}}}$, where $k$ is the
Boltzmann constant, $T$ is the temperature of the gas, and $m_{\mathrm{H}}$ is
the mass of the hydrogen atom. For $T=20,000\,\mathrm{K}$, we obtain
$v_{\mathrm{Dop}}\sim 18\,\mathrm{km\,s^{-1}}$. We assume the same
$v_{\mathrm{Dop}}$ also for the neutral gas. Substituting Equation (2) into
Equation (1) and using the Einstein coefficients taken from Wiese & Fuhr
(2009), we obtain

$\begin{split}\tau_{\mathrm{H}\alpha}&\sim 1.0\left(\frac{N_{2}}{2.9\times
10^{13}\,\mathrm{cm^{-3}}}\right)\left(\frac{v_{\mathrm{Dop}}}{18\,\mathrm{km\,s^{-1}}}\right)^{-1}\\\
\tau_{\mathrm{H}\beta}&\sim 0.14\left(\frac{N_{2}}{2.9\times
10^{13}\,\mathrm{cm^{-3}}}\right)\left(\frac{v_{\mathrm{Dop}}}{18\,\mathrm{km\,s^{-1}}}\right)^{-1}\\\
\tau_{\mathrm{H}\gamma}&\sim 0.05\left(\frac{N_{2}}{2.9\times
10^{13}\,\mathrm{cm^{-3}}}\right)\left(\frac{v_{\mathrm{Dop}}}{18\,\mathrm{km\,s^{-1}}}\right)^{-1}\\\
\tau_{\mathrm{H}\delta}&\sim 0.02\left(\frac{N_{2}}{2.9\times
10^{13}\,\mathrm{cm^{-3}}}\right)\left(\frac{v_{\mathrm{Dop}}}{18\,\mathrm{km\,s^{-1}}}\right)^{-1},\end{split}$
(3)

We find that if $v_{\mathrm{Dop}}=18\,\mathrm{km\,s^{-1}}$ is assumed, $N_{2}$
as large as $2.9\times 10^{13}\,\mathrm{cm^{-2}}$ makes $\mathrm{H}\alpha$
optically thick, while the other Balmer lines remain optically thin. In this
case, observed Balmer line ratios are calculated as

$\frac{\mathrm{H}n}{\mathrm{H}\beta}_{\mathrm{obs}}=\frac{\mathrm{H}n}{\mathrm{H}\beta}_{\mathrm{int}}\times\exp\left({-\tau_{\mathrm{H}n}}+\tau_{\mathrm{H}\beta}\right),$
(4)

where $n$ represents $\alpha$, $\gamma$, or $\delta$. We calculate
$(\mathrm{H}n/\mathrm{H}\beta)_{\mathrm{int}}$ using PyNeb and optical depths
are taken from Equation (3), where we assume
$v_{\mathrm{Dop}}=18\,\mathrm{km\,s^{-1}}$. In Figure 3, the Balmer line
ratios attenuated by neutral hydrogen are plotted for several $N_{2}$ and dust
extinction $A_{\mathrm{V}}$ values. Note that we treat $N_{2}$ and
$A_{\mathrm{V}}$ as the independent free parameters, because the gas-to-dust
ratio is unknown for the high-$z$ metal-poor galaxies. Here an SMC extinction
law is assumed (Gordon et al., 2003). RXCJ2248-ID shows broad components in
$\mathrm{H}\alpha$ and [O iii] $\lambda$5007 with FWHMs of 607 and 1258
$\mathrm{km\,s^{-1}}$, respectively, which implies existence of interstellar
shock. The effect of interstellar shock thus needs to be incorporated in the
model. In Figure 3, the cyan, blue, and purple shaded regions present the
effect of interstellar shock on $\mathrm{H}\alpha/\mathrm{H}\beta$ suggested
by Shull & McKee (1979), who present that interstellar shock increase
$\mathrm{H}\alpha/\mathrm{H}\beta$ by up to $\sim 1.6$ times. Figure 3
presents that larger $N_{2}$ leads to smaller
$\mathrm{H}\alpha/\mathrm{H}\beta$ and larger
$\mathrm{H}\gamma/\mathrm{H}\beta$ and $\mathrm{H}\delta/\mathrm{H}\beta$.
This is because the line with the longer wavelength has larger probability of
absorption (Equation 3). Comparing the models with the observational results,
we find that Balmer line ratios of GS-NDG-9422 and RXCJ2248-ID are reproduced
for $\log N_{2}\sim 12-13$ and $A_{\mathrm{V}}\sim 0$ within the $\sim
1\sigma$ levels. Although the absolute values of $N_{2}$ depend on the
assumption of the $v_{\mathrm{Dop}}$ value, this ambiguity does not
qualitatively change the following discussion.

Figure 3: Comparisons between the observed Balmer line ratios and theoretical
Balmer line ratios including the effects of the optical depths of the Balmer
lines. The squares are the same as Figure 1. The black crosses indicate the
Case B values for $T_{\mathrm{e}}=20,000\,\mathrm{K}$ and
$n_{\mathrm{e}}=10^{3}\,\mathrm{cm^{-3}}$. The cyan, blue, and purple lines
and circles denote the Balmer line ratios including the effects of the optical
depths with $A_{\mathrm{V}}=0$, $1$, and $3$, respectively. The numbers shown
beside the circles represent the values of $\log N_{2}$. The cyan, blue, and
purple shaded region in the left panel indicate the regions that can be
reproduced by the interstellar shock (Shull & McKee, 1979) in the case of
$A_{\mathrm{V}}=0,1$, and $3$, respectively.

####  3.3.2 MCMC Analyses

To estimate $N_{2}$ self-consistently with other physical parameters, we use
Markov Chain Monte Carlo (MCMC) algorithm. We utilized ymcmc code developed by
Hsyu et al. (2020). Although ymcmc have originally been developed to derive
helium abundances of galaxies, we rewrite the code to model hydrogen line
fluxes with the optical depth effects. Hereafter this code is referred to as
modified ymcmc. We set 5 free parameters of $T_{\mathrm{e}}$,
$n_{\mathrm{e}}$, a reddening parameter $c(\mathrm{H}\beta)$, a ratio of
neutral to ionized hydrogen density $\xi$, and $N_{2}$. At each step of the
MCMC, model flux ratios are calculated by

$\frac{F(\lambda)}{F(\mathrm{H}\beta)}_{\mathrm{mod}}=\frac{E(\lambda)}{E(\mathrm{H}\beta)}\frac{1+\frac{C}{R}(\lambda)}{1+\frac{C}{R}(\mathrm{H}\beta)}10^{-f(\lambda)\,c(\mathrm{H}\beta)}\\\
\times\exp{\left(-\tau_{\lambda}+\tau_{\mathrm{H}\beta}\right)},$ (5)

where $F(\lambda)$, $E(\lambda)$, $C/R(\lambda)$, and $f(\lambda)$ represent
the flux, emissivity, collisional-to-recombination coefficient, and reddening
law, respectively. See Hsyu et al. (2020) for the details of these parameters.
We minimize a log-likelihood function given by

$\log\left(\mathcal{L}\right)=\sum_{\lambda}\frac{\left(\frac{F(\lambda)}{F\left(\mathrm{H}\beta\right)}_{\mathrm{obs}}-\frac{F(\lambda)}{F\left(\mathrm{H}\beta\right)}_{\mathrm{mod}}\right)^{2}}{\sigma(\lambda)^{2}},$
(6)

where the subscript ‘obs’ denotes the observed value and $\sigma(\lambda)$
represents the $1\sigma$ uncertainty of the observed flux of the line
$\lambda$. We use flat priors of

$\begin{gathered}10,000\leq T_{\mathrm{e}}\leq 30,000\\\
1\leq\log_{10}n_{\mathrm{e}}\leq 5\\\ 0\leq c(\mathrm{H}\beta)\leq 0.5\\\
-6\leq\log_{10}\xi\leq-1\\\ 9\leq\log_{10}N_{2}\leq 15.\end{gathered}$ (7)

We do not use $\mathrm{H}\alpha$ for RXCJ2248-ID to avoid contamination by the
interstellar shock. The probability distributions for GS-NDG-9422 and
RXCJ2248-ID are presented in Figure 4 and 5, respectively. The best-fit
parameters are shown in Table 3. We obtain
$\log\left(N_{2}/\mathrm{cm^{-2}}\right)=12.5^{+0.3}_{-0.3}$ and
$13.7^{+0.2}_{-0.2}$ for GS-NDG-9422 and RXCJ2248-ID, respectively. The
probability distributions for $c(\mathrm{H}\beta)$ prefer values close to 0,
implying dust-poor conditions. These results are consistent with those seen in
Figure 3. The electron temperatures and densities are not constrained because
the Balmer line ratios are only weakly dependent on these parameters. In
Appendix A, we apply Gaussian priors for $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$
that are determined by the line diagnostics of the previous studies, and
obtain constraints on these parameters. Whether the Gaussian priors are
applied or not, our main results do not change.

Table 3: Best-fit parameters recovered with modified ymcmc

ID | $T_{\mathrm{e}}$ | $\log\left(n_{\mathrm{e}}/\mathrm{cm^{-3}}\right)$ | $c\left(\mathrm{H}\beta\right)$ | $\log\xi$ | $\log\left(N_{2}/\mathrm{cm^{-2}}\right)$
---|---|---|---|---|---
| [K] |  |  |  |
GS-NDG-9422 | … | … | $0.06^{+0.09}_{-0.05}$ | $-4.53^{+1.42}_{-0.99}$ | $12.49^{+0.27}_{-0.30}$
RXCJ2248-ID | … | … | $0.18^{+0.20}_{-0.12}$ | $-3.96^{+1.83}_{-1.42}$ | $13.67^{+0.20}_{-0.22}$

Note. — Median and the 16th and 84th percentiles of the probability
distributions are shown. Values for $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$ are
not presented because the obtained distributions are almost constant.

Figure 4: Probability distributions for GS-NDG-9422 recovered with modified
ymcmc. The diagonal panels show the 1D probability distribution for each
parameter. The off-diagonal panels present the 2D contours. The blue solid and
black dashed lines indicate the median and the 16th and 84th percentiles of
the probability distributions, respectively. Figure 5: Same as Figure 4, but
for RXCJ2248-ID.

##  4 Discussion

In the case where the Lyman lines are optically thin discussed in Section 3.2,
the nebula is density-bounded. This indicates the strong ionizing radiation in
the high-$z$ galaxies, which is in line with the large O32 value reported in
RXCJ2248-ID (Topping et al., 2024). Also in this case, the nebula can provide
a large number of $\mathrm{Ly}\alpha$ photons, which is compatible with the
strong $\mathrm{Ly}\alpha$ emission observed in GS-NDG-9422 (Cameron et al.,
2023) and RXCJ2248-ID (Mainali et al., 2017). The density-bounded nebulae can
provide ionizing photon with high escape fraction, which may play a key role
for cosmic reionization (Nakajima & Ouchi, 2014; Naidu et al., 2022). In the
density-bounded situation, however, it is difficult to explain the strong two-
photon emission observed in GS-NDG-9422 (Cameron et al., 2023). The two-photon
emission occurs in the forbidden $2s\rightarrow 1s$ transition. The large
$N_{2}$ value is needed to produce the strong two-photon emission, which is
difficult to achieve in the density-bounded nebula due to the low column
density. An existence of a damped Ly$\alpha$ (DLA) system with 30% leakage is
also suggested as an alternative explanation for the spectral feature of GS-
NDG-9422, while Cameron et al. (2023) regard it as unlikely. Although the
density-bounded nebulae can coexist with the DLA, scattering or absorption of
the Balmer line photons by DLA are likely to occur, which change the observed
Balmer line ratios. The condition with optically-thin Lyman lines does not
contradict the dense situation observed in RXCJ2248-ID (Topping et al., 2024).

In the case where the Balmer lines are optically-thick discussed in Section
3.3, $N_{2}$ needs to be large to make Balmer lines optically thick. To
increase $N_{2}$ in the neutral hydrogen gas, a large number of hydrogen atoms
at the ground state need to be excited to the second energy level. Although a
physical origin of this excessive excitation is not constrained, our result
may suggest that supernovae or accretion from the large scale structure to a
halo kinematically cause collisional excitation of the hydrogen gas in the
high-$z$ galaxies (Dekel et al., 2009; Yajima et al., 2015; Shimoda & Laming,
2019). Our result of large $N_{2}$ qualitatively compatible with the strong
two-photon emission observed in GS-NDG-9422 (Cameron et al., 2023), because a
large number of hydrogen atoms at the $2s$ state produce strong two-photon
emission via the $2s\rightarrow 1s$ transition. The neutral gas with
optically-thick Balmer lines can coexist with the dense ionized regions that
are reported in RXCJ2248-ID (Topping et al., 2024).

##  5 Summary

In this paper, the physical mechanisms that cause the observed Balmer
decrement anomalies are investigated. Our main findings are summarized below:

  1. 1.

GS-NDG-9422 and RXCJ2248-ID found with JWST at $z\sim 6$ show
$\mathrm{H}\alpha/\mathrm{H}\beta$ significantly smaller than the Case B
values as reported by Cameron et al. (2023) and Topping et al. (2024),
respectively. RXCJ2248-ID also presents $\mathrm{H}\gamma/\mathrm{H}\beta$ and
$\mathrm{H}\delta/\mathrm{H}\beta$ significantly larger than the Case B
values. The observed Balmer line ratios are not explained by the combination
of the dust extinction and Case B, since the observed deviation from the Case
B values are opposite to the effect of the dust extinction. Changing the
electron temperature and density also does not reproduce the observed Balmer
line ratios.

  2. 2.

Using Cloudy, we model the density-bounded nebula where the Lyman lines are
optically thin. We calculate the Balmer line ratios and the optical depth of
Ly$\alpha$ with several stopping thicknesses. We find that the optical depths
of Lyman series changes the Balmer line ratios due to the conversions of the
Lyman lines to Balmer lines. We find that when the nebula is optically thick
only up to Ly$\gamma-$Ly8, the observed Balmer line ratios are reproduced.
However, this situation requires fine-tuning of the thickness of the nebula.
The density-bounded situation suggests the strong ionizing radiation, which
qualitatively agree with the large O32 observed in RXCJ2248-ID. The nebula
that is optically-thin to Lyman lines is also compatible with the strong
Ly$\alpha$ emission observed in GS-NDG-9422. The density-bounded nebulae
supply ionizing photon with high escape fraction, which may contribute to the
cosmic reionization. However, the density-bounded situation is difficult to
explain the strong two-photon continuum observed in GS-NDG-9422.

  3. 3.

We also model the neutral hydrogen gas that surround the ionized region and is
optically thick to the Balmer lines. We find that, if $N_{2}$ is larger than
$\sim 10^{13}\,\mathrm{cm^{-2}}$, the Balmer lines become optically thick, and
that the observed Balmer line ratios are reproduced. We develop nebular models
with neutral gas that surround the ionized gas and have large $N_{2}$. Our
nebular models extending the ymcmc models (Hsyu et al., 2020) suggest
$\log\left(N_{2}/\mathrm{cm^{-2}}\right)=12.5^{+0.3}_{-0.3}$ and
$13.7^{+0.2}_{-0.2}$ for GS-NDG-9422 and RXCJ2248-ID, respectively, as well as
the dust-poor conditions. This results suggest that the excited H i clouds
that partly (or fully) surround the ionized region may be the cause of the
Balmer decrement anomalies. Although a physical process that enhances $N_{2}$
in the neutral gas around the ionized gas is unclear, there may be accretion
or supernovae that cause excessive collisional excitation of the second energy
level of hydrogen atoms. Our result of large $N_{2}$ for GS-NDG-9422 in line
with the strong two-photon continuum suggested by Cameron et al. (2023).
Because a large number of hydrogen atoms at the $2s$ state produce the two-
photon emission via the $2s\rightarrow 1s$ transition. The neutral hydrogen
nebula with high $N_{2}$ can coexist the dense ionized nebula reported in
RXCJ2248-ID (Topping et al., 2024).

## acknowledgments

We thank Michael W. Topping for providing the information about the
$\mathrm{H}\alpha$ flux reported by Topping et al. (2024). We are thankful to
Kentaro Nagamine for helpful comments. This work use the results reported by
Alex J. Cameron and Michael W. Topping. We acknowledge Alex J. Cameron and
Michael W. Topping for publicly releasing their analysed data. This work is
based on the observations with the NASA/ESA/CSA James Webb Space Telescope
associated with programs of GTO-1210 (JADES) and GO-2478. We are grateful to
GTO-1210 and GO-2478 teams led by Nora Luetzgendorf and Daniel P. Stark,
respectively. This publication is based upon work supported by the World
Premier International Research Center Initiative (WPI Initiative), MEXT,
Japan, KAKENHI (20H00180, 21H04467, 21H04489, 21H04496, and 23H05441) through
Japan Society for the Promotion of Science, and JST FOREST Program (JP-
MJFR202Z). This work was supported by the joint research program of the
Institute for Cosmic Ray Research (ICRR), University of Tokyo.

##  Appendix A Gaussian Priors for $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$ in
MCMC Analyses

The flat priors (7) give the almost constant probability distributions for
$T_{\mathrm{e}}$ and $n_{\mathrm{e}}$, since the dependence of Balmer
decrements on these parameters are weak. However, providing independent
measurements of $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$ may give physically more
reasonable distributions. Here, using $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$
derived from emission line diagnostics, we put weak Gaussian priors on these
parameters. This procedure is proposed by Aver et al. (2011), who have used
mock data and demonstrated that the input model parameters are recovered more
accurately when the weak Gaussian prior is imposed. Given the $T_{\mathrm{e}}$
and $n_{\mathrm{e}}$ values derived from the emission line diagnostics, log-
likelihood is calculated as

$\log\left(\mathcal{L}\right)=\sum_{\lambda}\frac{\left(\frac{F(\lambda)}{F\left(\mathrm{H}\beta\right)}_{\mathrm{obs}}-\frac{F(\lambda)}{F\left(\mathrm{H}\beta\right)}_{\mathrm{mod}}\right)^{2}}{\sigma(\lambda)^{2}}+\frac{(T_{\mathrm{e,obs}}-T_{\mathrm{e,mod}})^{2}}{2\sigma_{T_{\mathrm{e}}}^{2}}+\frac{(n_{\mathrm{e,obs}}-n_{\mathrm{e,mod}})^{2}}{2\sigma_{n_{\mathrm{e}}}^{2}},$
(A1)

where $T_{\mathrm{e,obs}}$ and $n_{\mathrm{e,obs}}$ are the electron
temperature and density derived from observation, respectively, and
$T_{\mathrm{e,mod}}$ and $n_{\mathrm{e,mod}}$ are those of the model
parameters. $\sigma_{T_{\mathrm{e}}}$ and $\sigma_{n_{\mathrm{e}}}$ are the
standard deviations of the Gaussian distribution. Here we take
$\sigma_{T_{\mathrm{e}}}$ and $\sigma_{n_{\mathrm{e}}}$ to be conservative
values of $0.2\,T_{\mathrm{e,obs}}$ and $0.4\,n_{\mathrm{e,obs}}$,
respectively, the former of which is proposed by Aver et al. (2011). The
obtained probability distributions for GS-NDG-9422 and RXCJ2248-ID are
presented in Figure 6 and 7, respectively. The $T_{\mathrm{e,obs}}$ and
$n_{\mathrm{e,obs}}$ values and best-fit parameters are summarized in Table 4.
Note that $n_{\mathrm{e,obs}}$ is not used for GS-NDG-9422 (the last term in
right hand side of Equation (A1) is omitted) because the only weak upper limit
on $n_{\mathrm{e}}$ is provided (Cameron et al., 2023). The best-fit $N_{2}$
values are $\sim 12$ and $\sim 13$ for GS-NDG-9422 and RXCJ2248-ID,
respectively, indicating that imposing the weak Gaussian priors does not
change our main result.

Table 4: Best-fit parameters recovered with modified ymcmc

ID | $T_{\mathrm{e}}$ | $\log\left(n_{\mathrm{e}}/\mathrm{cm^{-3}}\right)$ | $c\left(\mathrm{H}\beta\right)$ | $\log\xi$ | $\log\left(N_{2}/\mathrm{cm^{-2}}\right)$
---|---|---|---|---|---
| [K] |  |  |  |
GS-NDG-9422 | $18278^{+2454}_{-2514}$ | $\dots$ | $0.06^{+0.10}_{-0.05}$ | $-4.38^{+1.34}_{-1.14}$ | $12.49^{+0.27}_{-0.28}$
RXCJ2248-ID | $24049^{+3132}_{-3284}$ | $4.01^{+0.67}_{-1.03}$ | $0.18^{+0.19}_{-0.14}$ | $-4.10^{+1.69}_{-1.28}$ | $13.45^{+0.19}_{-0.25}$

Note. — Same as Table 3, but with Gaussian priors for $T_{\mathrm{e}}$ and
$n_{\mathrm{e}}$. For GS-NDG-9422, since $n_{\mathrm{e}}$ derived from line
diagnostics does not exist, no constraint on $n_{\mathrm{e}}$ is obtained.

Figure 6:  Same as Figure 4, but for GS-NDG-9422 with Gaussian priors for
$T_{\mathrm{e}}$.  Figure 7: Same as Figure 4, but for RXCJ2248-ID with
Gaussian priors for $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$.

## References

  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Aver et al. (2011) Aver, E., Olive, K. A., & Skillman, E. D. 2011, J. Cosmology Astropart. Phys, 2011, 043, doi: 10.1088/1475-7516/2011/03/043
  * Bunker et al. (2023) Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2306.02467, doi: 10.48550/arXiv.2306.02467
  * Cameron et al. (2023) Cameron, A. J., Katz, H., Witten, C., et al. 2023, arXiv e-prints, arXiv:2311.02051, doi: 10.48550/arXiv.2311.02051
  * Chatzikos et al. (2023) Chatzikos, M., Bianchi, S., Camilloni, F., et al. 2023, Rev. Mexicana Astron. Astrofis., 59, 327, doi: 10.22201/ia.01851101p.2023.59.02.12
  * Dekel et al. (2009) Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451, doi: 10.1038/nature07648
  * Eisenstein et al. (2023) Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv e-prints, arXiv:2306.02465, doi: 10.48550/arXiv.2306.02465
  * Gordon et al. (2003) Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279, doi: 10.1086/376774
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hsyu et al. (2020) Hsyu, T., Cooke, R. J., Prochaska, J. X., & Bolte, M. 2020, ApJ, 896, 77, doi: 10.3847/1538-4357/ab91af
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Luridiana et al. (2015) Luridiana, V., Morisset, C., & Shaw, R. A. 2015, A&A, 573, A42, doi: 10.1051/0004-6361/201323152
  * Mainali et al. (2017) Mainali, R., Kollmeier, J. A., Stark, D. P., et al. 2017, ApJ, 836, L14, doi: 10.3847/2041-8213/836/1/L14
  * Naidu et al. (2022) Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022, ApJ, 940, L14, doi: 10.3847/2041-8213/ac9b22
  * Nakajima & Ouchi (2014) Nakajima, K., & Ouchi, M. 2014, MNRAS, 442, 900, doi: 10.1093/mnras/stu902
  * Osterbrock & Ferland (2006) Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Shimoda & Laming (2019) Shimoda, J., & Laming, J. M. 2019, MNRAS, 485, 5453, doi: 10.1093/mnras/stz758
  * Shull & McKee (1979) Shull, J. M., & McKee, C. F. 1979, ApJ, 227, 131, doi: 10.1086/156712
  * Topping et al. (2024) Topping, M. W., Stark, D. P., Senchyna, P., et al. 2024, arXiv e-prints, arXiv:2401.08764, doi: 10.48550/arXiv.2401.08764
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wiese & Fuhr (2009) Wiese, W. L., & Fuhr, J. R. 2009, Journal of Physical and Chemical Reference Data, 38, 565, doi: 10.1063/1.3077727
  * Yajima et al. (2015) Yajima, H., Li, Y., Zhu, Q., & Abel, T. 2015, ApJ, 801, 52, doi: 10.1088/0004-637X/801/1/52
  * Yano et al. (2021) Yano, K., Baba, S., Nakagawa, T., et al. 2021, ApJ, 922, 272, doi: 10.3847/1538-4357/ac26be