Alternative Minimal U(1)B-L

Kento Asai(a), Kazunori Nakayama(a,b), Shih-Yen Tseng(a)

(a)Department of Physics, Faculty of Science,
The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan
(b)Kavli IPMU (WPI), The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

###### Abstract

We provide a minimal alternative gauged U(1)B-L model in which three right-
handed neutrinos with charges of $(5,-4,-4)$ and only one ${\rm B-L}$ Higgs
field with charge $1$ are introduced. Consistent active neutrino masses and
mixings can be obtained if the $Z_{2}$ symmetry on two of the right-handed
neutrinos is introduced. It predicts two heavy degenerate right-handed
neutrinos, which may realize resonant leptogenesis scenario, and one
relatively light sterile neutrino, which is a good dark matter candidate.

##  1 Introduction

Observed active neutrino masses and mixings are considered to be naturally
explained by the seesaw mechanism [1, 2, 3]. In the seesaw mechanism right-
handed neutrinos (RHNs) are introduced which are heavy and have yukawa
couplings to the Standard Model (SM) leptons. After integrating out RHNs tiny
active neutrino masses are obtained. RHNs not only explains the neutrino
masses but also can create the cosmological baryon asymmetry through
leptogenesis [4].

On the other hand, the existence of RHNs is required if the U(1)B-L symmetry
is gauged, in order to cancel the gauge anomaly. Conventionally three RHNs
$N_{i}$ $(i=1,2,3)$ are introduced with the ${\rm B-L}$ charge $(-1,-1,-1)$.
We call this case as “conventional U(1)B-L”. Phenomenology of such a scenario
is well studied. However, there is another possibility for the charge
assignments on three RHNs consistent with the anomaly cancellation:
$(5,-4,-4)$ [5]. We call this case as the “alternative U(1)B-L”. In this paper
we focus on the alternative U(1)B-L.

In the alternative U(1)B-L model, it is nontrivial to have appropriate RHN
masses and yukawa couplings for the seesaw mechanism because of the unusual
${\rm B-L}$ charge assignments. In order to overcome this difficulty, several
${\rm B-L}$ Higgs fields with nontrivial charges are often introduced.
Phenomenology of such alternative U(1)B-L has been studied in, e.g. Refs. [6,
7, 8, 9, 10, 11, 12, 13, 14, 15].

In this paper, we propose a alternative minimal U(1)B-L model in which there
is only one ${\rm B-L}$ Higgs field with the ${\rm B-L}$ charge $+1$. It is
minimal in the sense that we do not introduce any additional field except for
three RHNs and one ${\rm B-L}$ Higgs field. The only extra assumption is that
we will impose a $Z_{2}$ symmetry under which two RHNs are exchanged:
$N_{2}\leftrightarrow N_{3}$. It drastically affects the yukawa structure and
it can lead to correct active neutrino masses and mixings through the seesaw
mechanism. An interesting prediction of this model is the characteristic
structure of RHN masses: two RHNs are heavy and nearly degenerate in mass and
the other RHN is hierarchically light. The light RHN is stable and a candidate
of dark matter (DM) of the universe.

In Sec. 2 we introduce the minimal alternative U(1)B-L model without $Z_{2}$
symmetry and explain its difficulty to generate active neutrino masses. In
Sec. 3 we introduce a $Z_{2}$ symmetry and see that it solves the problem of
the model without $Z_{2}$. In Sec. 4 we discuss several phenomenological
implications.

##  2 Alternative minimal U(1)B-L without $Z_{2}$

First we introduce a minimal alternative U(1)B-L without $Z_{2}$ and explain
its basic structure. Eventually it will fail to explain the active neutrino
masses, but it is instructive for later discussion.

Minimal alternative U(1)B-L charge assignments are given in Table. 1, where
$Q_{Li},u_{Ri},d_{Ri}$ with $i=1$–$3$ are the SM quarks, $L_{\alpha}$,
$e_{R\alpha}$ $(\alpha=e,\mu,\tau)$ are the SM leptons, $N_{i}$ $(i=1,2,3)$
are RHNs, $H$ is the SM Higgs doublet and $\phi$ is the ${\rm B-L}$ Higgs
field. Since we introduce only one ${\rm B-L}$ Higgs field, we need higher
dimensional operators to obtain yukawa couplings and some elements of RHN mass
matrix. For example, $L_{\alpha}$ and $N_{1}$ can have yukawa couplings of the
form
$\mathcal{L}\sim(\phi^{*}/\Lambda)^{6}\overline{L}_{\alpha}N_{1}\widetilde{H}$
with $\Lambda$ being the cutoff scale. After taking the vacuum expectation
value (VEV) of $\phi$, the relevant Lagrangian is given by

$\displaystyle\mathcal{L}_{Y}=\sum_{\alpha=e,\mu,\tau}\sum_{i=1,2,3}y_{\alpha
i}\overline{L}_{\alpha}N_{i}\widetilde{H}+{\rm h.c.},$ (1)
$\displaystyle\mathcal{L}_{N}=-\sum_{i,j=1,2,3}\frac{1}{2}M_{ij}N_{i}N_{j}+{\rm
h.c.},$ (2)

where

$\displaystyle y_{\alpha
i}\sim\begin{pmatrix}\epsilon^{6}&\epsilon^{3}&\epsilon^{3}\\\
\epsilon^{6}&\epsilon^{3}&\epsilon^{3}\\\
\epsilon^{6}&\epsilon^{3}&\epsilon^{3}\end{pmatrix},~{}~{}~{}~{}~{}~{}~{}~{}~{}M_{ij}\sim
v_{\phi}\begin{pmatrix}\epsilon^{9}&1&1\\\ 1&\epsilon^{7}&\epsilon^{7}\\\
1&\epsilon^{7}&\epsilon^{7}\end{pmatrix},$ (3)

omitting all $\mathcal{O}(1)$ coefficients. Here $\epsilon\equiv
v_{\phi}/\Lambda$ with $v_{\phi}\equiv\left<|\phi|\right>$. The detailed
structure of the scalar potential $V(\phi,H)$ is not relevant for our study,
but for completeness we give it here:

$\displaystyle
V(\phi,H)=\lambda_{H}\left(|H|^{2}-v_{H}^{2}\right)^{2}+\lambda_{\phi}\left(|\phi|^{2}-v_{\phi}^{2}\right)^{2}+\lambda_{H\phi}\left(|H|^{2}-v_{H}^{2}\right)\left(|\phi|^{2}-v_{\phi}^{2}\right).$
(4)

If $\lambda_{H}>0$ and $\lambda_{\phi}>0$, the potential minimum is
$\left<|H|^{2}\right>=v_{H}^{2}$ and $\left<|\phi|\right>=v_{\phi}$ for
$\lambda_{H\phi}^{2}<4\lambda_{H}\lambda_{\phi}$. It is possible to cure the
SM Higgs vacuum instability if
$\lambda_{H\phi}^{2}/(4\lambda_{\phi})\gtrsim\lambda_{H}$ and
$\sqrt{\lambda_{\phi}}v_{\phi}\lesssim 10^{10}\,{\rm GeV}$ [17, 18], although
we do not go into details of this subject.

$Q_{Li}$ | $u_{Ri}$ | $d_{Ri}$ | $L_{\alpha}$ | $e_{R\alpha}$ | $N_{1}$ | $N_{2}$ | $N_{3}$ | $H$ | $\phi$
---|---|---|---|---|---|---|---|---|---
$1/3$ | $1/3$ | $1/3$ | $-1$ | $-1$ | $5$ | $-4$ | $-4$ | $0$ | $1$
Table 1: Alternative U(1)B-L charge assignments on the SM quarks
($Q_{Li},u_{Ri},d_{Ri}$ with $i=1$–$3$), SM leptons ($L_{\alpha},e_{R\alpha}$
with $\alpha=e,\mu,\tau$), RHNs ($N_{1},N_{2},N_{3}$), SM Higgs ($H$) and
${\rm B-L}$ Higgs ($\phi$).

From the RHN mass matrix shown in (3), one finds that the mass eigenvalues
$\widetilde{M}_{i}$ $(i=1,2,3)$ and eigenstates $\widetilde{N}_{i}$ are given
by

$\displaystyle\widetilde{M}_{1}\simeq\widetilde{M}_{2}\simeq v_{\phi}$
$\displaystyle~{}~{}~{}{\rm
for}~{}~{}~{}\widetilde{N}_{1,2}\simeq\frac{1}{\sqrt{2}}(N_{1}\pm N_{\rm H}),$
(5) $\displaystyle\widetilde{M}_{3}\simeq\epsilon^{7}v_{\phi}$
$\displaystyle~{}~{}~{}{\rm for}~{}~{}~{}\widetilde{N}_{3}\simeq N_{\rm L},$
(6)

to the leading order in $\epsilon$, where

$\displaystyle N_{\rm
H}\equiv\frac{1}{\sqrt{2}}(N_{2}+N_{3}),~{}~{}~{}~{}~{}~{}N_{\rm
L}\equiv\frac{1}{\sqrt{2}}(-N_{2}+N_{3}).$ (7)

Therefore, for $\epsilon\lesssim 0.1$, one of the RHNs becomes very light and
the split mass spectrum is naturally realized. Note that $N_{1}$ and $N_{\rm
H}$ forms a Dirac mass term and hence $\widetilde{M}_{1}$ and
$\widetilde{M}_{2}$ are almost degenerate.

However, this model fails to explain the observed neutrino masses. In (3), we
see that $N_{2}$ and $N_{3}$ have yukawa couplings of similar order and they
are much larger than the ones of $N_{1}$, only $N_{\rm L}$ dominantly
contributes to the active neutrino masses after the seesaw mechanism. Thus two
active neutrinos become almost massless while only one active neutrino can
obtain an appropriate mass. Actually one can show that the active neutrino
masses are given by

$\displaystyle m_{\nu}\sim\frac{v_{H}^{2}}{\epsilon
v_{\phi}}(\epsilon^{9},\epsilon^{9},1).$ (8)

It cannot explain the observed neutrino mass differences unless $\epsilon$ is
very close to $1$, which, however, is just a limit that approaches the
conventional U(1)B-L model and not of our interest here.

##  3 Alternative minimal U(1)B-L with $Z_{2}$

In the previous section we have seen that the minimal alternative U(1)B-L is
not consistent with the observed neutrino masses. The problem is that one
light RHN $(N_{\rm L})$ dominantly contributes to the seesaw mechanism and
hence there appear two (almost) massless active neutrino species. As we shall
see below, one can decouple $N_{\rm L}$ from the SM yukawa sector by
introducing an appropriate $Z_{2}$ parity so that the remaining two heavy RHNs
contribute to the seesaw mechanism.

Let us introduce $Z_{2}$ parity under which $N_{2}$ and $N_{3}$ are exchanged:

$\displaystyle Z_{2}~{}:~{}N_{2}\leftrightarrow N_{3}.$ (9)

In terms of the irreducible representation of $Z_{2}$, it is equivalent to

$\displaystyle Z_{2}~{}:~{}N_{\rm H}\rightarrow N_{\rm
H},~{}~{}~{}~{}~{}~{}N_{\rm L}\rightarrow-N_{\rm L},$ (10)

where

$\displaystyle N_{\rm
H}\equiv\frac{1}{\sqrt{2}}(N_{2}+N_{3}),~{}~{}~{}~{}~{}~{}N_{\rm
L}\equiv\frac{1}{\sqrt{2}}(-N_{2}+N_{3}).$ (11)

The Lagrangian is of the form (1)–(2), but the structure of the yukawa and RHN
mass matrix are constrained due to the $Z_{2}$ symmetry. First let us see the
RHN mass matrix. Due to the $Z_{2}$ symmetry, the mass matrix of RHNs have the
structure of

$\displaystyle M_{ij}=v_{\phi}\begin{pmatrix}c_{1}\epsilon^{9}&c_{2}&c_{2}\\\
c_{2}&c_{3}\epsilon^{7}&c_{4}\epsilon^{7}\\\
c_{2}&c_{4}\epsilon^{7}&c_{3}\epsilon^{7}\end{pmatrix},$ (12)

where $c_{1},c_{2},c_{3},c_{4}$ are $\mathcal{O}(1)$ coefficients. The mass
eigenvalues of the RHNs are similar to the case without $Z_{2}$:

$\displaystyle\widetilde{M}_{1}\simeq\left(\sqrt{2c_{2}^{2}}+\frac{c_{3}+c_{4}}{2}\epsilon^{7}\right)v_{\phi}$
$\displaystyle~{}~{}~{}{\rm
for}~{}~{}~{}\widetilde{N}_{1}\simeq\frac{1}{\sqrt{2}}(N_{1}+N_{\rm H}),$ (13)
$\displaystyle\widetilde{M}_{2}\simeq\left(-\sqrt{2c_{2}^{2}}+\frac{c_{3}+c_{4}}{2}\epsilon^{7}\right)v_{\phi}$
$\displaystyle~{}~{}~{}{\rm
for}~{}~{}~{}\widetilde{N}_{2}\simeq\frac{1}{\sqrt{2}}(N_{1}-N_{\rm H}),$ (14)
$\displaystyle\widetilde{M}_{3}=(c_{3}-c_{4})\epsilon^{7}v_{\phi}$
$\displaystyle~{}~{}~{}{\rm for}~{}~{}~{}\widetilde{N}_{3}=N_{\rm L}.$ (15)

For later convenience, we define $M_{\rm H1}\equiv|\widetilde{M}_{1}|$,
$M_{\rm H2}\equiv|\widetilde{M}_{2}|$ and $M_{\rm
L}\equiv|\widetilde{M}_{3}|$. Note again that $M_{\rm H1}$ and $M_{\rm H2}$
are almost degenerate: the mass difference is $\mathcal{O}(\epsilon^{7})$. In
this case, $N_{\rm L}$ is exactly the mass eigenstate while the mixing of
$N_{\rm H}$ and $N_{1}$ are slightly modified by $\mathcal{O}(\epsilon^{7})$
terms. Similar to the previous case, two heavy RHNs are almost degenerate and
the other one $(N_{\rm L})$ is hierarchically light.

Next let us see the yukawa sector. Due to the $Z_{2}$ symmetry, yukawa
couplings must satisfy $y_{\alpha 2}=y_{\alpha 3}(\equiv y_{\alpha\rm
H}/\sqrt{2})$. Thus the yukawa terms are written as

$\displaystyle\mathcal{L}_{Y}=\sum_{\alpha=e,\mu,\tau}\left[y_{\alpha
1}\overline{L}_{\alpha}N_{1}\widetilde{H}+y_{\alpha\rm
H}\overline{L}_{\alpha}N_{\rm H}\widetilde{H}\right]+{\rm h.c.}.$ (16)

Here is the crucial observation: it is only $N_{1}$ and $N_{\rm H}$ that have
yukawa couplings to the active neutrinos. Since both $N_{1}$ and $N_{\rm H}$
are heavy and have almost the same mass, it effectively reduces to the seesaw
with two RHNs [19, 20]. Remarkably, the light RHN, $N_{\rm L}$, decouples from
the SM yukawa sector. This is due to the $Z_{2}$ symmetry under which only
$N_{\rm L}$ flips the sign and all the remaining fields do not. As far as the
$Z_{2}$ symmetry is exact, $N_{\rm L}$ is stable and a good DM candidate.

Now the active neutrino mass matrix is obtained through the seesaw formula:

$\displaystyle m_{\nu}=m_{\rm D}(M_{N})^{-1}m_{\rm D}^{\rm
T}\simeq\frac{\epsilon^{9}v_{H}^{2}}{\sqrt{2}c_{2}v_{\phi}}\begin{pmatrix}2d_{1}d_{4}&d_{2}d_{4}+d_{1}d_{5}&d_{3}d_{4}+d_{1}d_{6}\\\
d_{2}d_{4}+d_{1}d_{5}&2d_{2}d_{5}&d_{3}d_{5}+d_{2}d_{6}\\\
d_{3}d_{4}+d_{1}d_{6}&d_{3}d_{5}+d_{2}d_{6}&2d_{3}d_{6}\end{pmatrix},$ (17)

where

$\displaystyle m_{\rm
D}=v_{H}\begin{pmatrix}d_{1}\epsilon^{6}&d_{4}\epsilon^{3}\\\
d_{2}\epsilon^{6}&d_{5}\epsilon^{3}\\\
d_{3}\epsilon^{6}&d_{6}\epsilon^{3}\end{pmatrix},~{}~{}~{}~{}~{}~{}M_{N}=v_{\phi}\begin{pmatrix}c_{1}\epsilon^{9}&\sqrt{2}c_{2}\\\
\sqrt{2}c_{2}&(c_{3}+c_{4})\epsilon^{7}\end{pmatrix},$ (18)

with $d_{1}$–$d_{6}$ being $\mathcal{O}(1)$ coefficients.#1#1#1 The seesaw
structure is similar to those studied in Refs. [21, 22] in the context of
sneutrino chaotic inflation.  Note that the neutrino mass matrix is rank 2,
since the matrix $M_{N}$ is rank 2 and there is one massless neutrino. The
neutrino mass eigenvalues are given by

$\displaystyle
m_{\nu}\simeq\frac{\epsilon^{9}v_{H}^{2}}{\sqrt{2}c_{2}v_{\phi}}\begin{pmatrix}0,&d_{A}+d_{B},&d_{A}-d_{B}\end{pmatrix},$
(19)

where $d_{A}$ and $d_{B}$ are $\mathcal{O}(1)$ coefficients defined by

$\displaystyle
d_{A}\equiv\sum_{i=1,2,3}d_{i}d_{i+3},~{}~{}~{}~{}~{}~{}d_{B}\equiv\sqrt{\left(\sum_{i=1,2,3}d_{i}^{2}\right)\left(\sum_{j=4,5,6}d_{j}^{2}\right)}.$
(20)

Using (15) to eliminate $\epsilon$, the typical active neutrino mass is
estimated as

$\displaystyle m_{\nu}$ $\displaystyle=\left|\frac{d_{A}\pm
d_{B}}{\sqrt{2}c_{2}(c_{3}-c_{4})^{9/7}}\right|\frac{v_{H}^{2}M_{\rm
L}^{9/7}}{v_{\phi}^{16/7}}\equiv C_{\pm}\frac{v_{H}^{2}M_{\rm
L}^{9/7}}{v_{\phi}^{16/7}}$ (21) $\displaystyle\simeq 0.01\,{\rm eV}\times
C_{\pm}\left(\frac{M_{\rm L}}{10^{4}\,{\rm
GeV}}\right)^{9/7}\left(\frac{10^{9}\,{\rm GeV}}{v_{\phi}}\right)^{16/7}.$
(22)

Fig. 1 shows the contours of $\mathcal{O}(1)$ coefficient $C$ (one of the
$C_{\pm}$) on the $(v_{\phi},M_{\rm L})$ plane, which reproduces the observed
neutrino mass difference $|\Delta m_{23}|\sim 0.05\,{\rm eV}$ [23]. The other
observed neutrino mass difference $|\Delta m_{21}|$ is also reproduced if the
other coefficient of $C_{\pm}$ is slightly smaller.

In calculating the active neutrino mass matrix, all the $3\times 2$ elements
of the yukawa matrix equally contribute to the final expression of the
neutrino mass matrix and hence there are enough degrees of freedom to
reproduce the observed neutrino masses, mixings and CP phase. Yukawa couplings
(or $d_{1}$–$d_{6}$) are in general complex and hence contain $12$ real
parameters. Among them, $3$ phases can be absorbed by the phase rotation of
$L_{\alpha}$ and hence there are $9$ parameters remained. On the other hand,
observable quantities in the neutrino sector are $2$ neutrino masses (since
the lightest neutrino is massless), $3$ mixing angles, $1$ Dirac CP phase and
$1$ Majorana phase (one more Majorana phase is not physical when one of the
neutrinos is massless). Thus $9$ parameters in the model are enough to explain
the observed data.

---
Figure 1: The predictions on the $v_{\phi}$-$M_{\rm L}$ plane are shown in the
plot. A band of the predicted coefficient $C$ is given in the range ($0.1,10$)
corresponding to the observed mass difference $\left|\Delta
m_{23}\right|\simeq 0.05$ eV. We also present some $\epsilon$ contours and the
excluded region with $\epsilon>1$ is shaded.

##  4 Discussion

### Dark matter

In our $Z_{2}$ symmetric model, the light RHN ($N_{\rm L}$) is stable and a
candidate of DM. Production mechanisms $N_{\rm L}$ particles are highly
dependent on properties of DM and inflation models. There are various types of
production scenarios, depending on the reheating temperature $T_{\rm R}$ and
whether DM is in thermal bath in the early universe or not. In the thermal
freeze-in scenario [24, 25, 26, 27], $N_{\rm L}$ particles are produced
through the scattering of SM fermions in thermal bath with the ${\rm B-L}$
gauge boson exchange. The abundance is dominated by those created at the
temperature $T\sim T_{\rm R}$ ($T_{\rm R}$ denotes the reheating temperature),
assuming $T_{\rm R}\gtrsim M_{\rm L}$. The abundance in terms of the density
parameter is given by [24]#2#2#2 Since the ${\rm B-L}$ charge of light RHN is
$-4$, the abundance is $16$ times larger than the estimate in Ref. [24].

$\displaystyle\Omega_{N_{\rm L}}h^{2}\sim 0.5\left(\frac{M_{\rm
L}}{10^{3}\,{\rm GeV}}\right)\left(\frac{T_{\rm R}}{10^{3}\,{\rm
GeV}}\right)^{3}\left(\frac{10^{9}\,{\rm GeV}}{v_{\phi}}\right)^{4}.$ (23)

The ${\rm B-L}$ gauge boson mass is assumed to be much larger than the
temperature. If $T_{\rm R}\lesssim M_{\rm L}$, on the other hand, thermal
production is dominated by $T\sim M_{\rm L}$ before the completion of the
reheating, since there is high temperature dilute plasma and the maximum
temperature is much higher than $T_{\rm R}$ and possibly than $M_{\rm L}$. The
abundance in this case is given by

$\displaystyle\Omega_{N_{\rm L}}h^{2}\sim 0.5\left(\frac{10^{3}\,{\rm
GeV}}{M_{\rm L}}\right)^{3}\left(\frac{T_{\rm R}}{10^{3}\,{\rm
GeV}}\right)^{7}\left(\frac{10^{9}\,{\rm GeV}}{v_{\phi}}\right)^{4}.$ (24)

It can reproduce the observed DM abundance for reasonable choice of the
reheating temperature. Here it is implicitly assumed that $N_{\rm L}$ is not
thermalized, although for higher $T_{\rm R}$ and/or smaller $v_{\phi}$ it is
possible that $N_{\rm L}$ takes part in thermal bath. In such a case the final
relic abundance is determined by the $N_{\rm L}$ self-annihilation cross
section through the freeze-out mechanism [28]. Another possible source is the
decay of a scalar [29, 30, 31, 32]. In our model it may be possible to obtain
a correct relic abundance through the $\phi$ decay, but it strongly depends on
the dynamics of $\phi$ in the early universe. The detailed discussion will be
done in a future work.

### Leptogenesis

In our model, two heavy RHNs are almost degenerate in mass. The mass
difference is of the order of $\left(M_{\rm H1}-M_{\rm H2}\right)/M_{\rm
H1}\sim\mathcal{O}(\epsilon^{7})\sim M_{\rm L}/M_{\rm H1}$. If they are
produced thermally or non-thermally, their CP-violating decay can generate a
lepton asymmetry. Due to the degeneracy of two RHNs, the so-called resonant
leptogenesis may occur [33, 34, 35, 36, 37, 38, 39]. The CP asymmetry
parameter per RHN decay is given by

$\displaystyle\eta_{i}=\frac{{\rm
Im}\left[(\widetilde{y}_{i\alpha}\widetilde{y}^{\dagger}_{\alpha
j})^{2}\right]}{8\pi(\widetilde{y}_{i\alpha}\widetilde{y}^{\dagger}_{\alpha
i})}\frac{M_{{\rm H}i}M_{{\rm H}j}(M_{{\rm H}i}^{2}-M_{{\rm
H}j}^{2})}{(M_{{\rm H}i}^{2}-M_{{\rm H}j}^{2})^{2}+R^{2}},$ (25)

for $i=1,2$, where $\widetilde{y}_{i\alpha}$ denotes the yukawa matrix in the
RHN mass eigenstate basis and $R$ denotes the regulator, which is of the order
of $M_{\rm H}\Gamma$ with $\Gamma$ being the RHN decay width [40, 41, 42, 43,
44]. As a very rough estimate, $\Gamma\sim\epsilon^{6}M_{\rm H}/(8\pi)$, which
is more or less the same order as $M_{\rm L}$ and hence $R\sim M_{{\rm
H}i}^{2}-M_{{\rm H}j}^{2}$. Remarkably, in this model the CP asymmetry is
automatically maximized. Note that typically
$\widetilde{y}_{i\alpha}\sim\epsilon^{3}$, we have the CP parameter of
$\eta_{i}\sim\mathcal{O}(0.1)$.

The sphaleron effect converts the lepton asymmetry into the baryon asymmetry
[45]. The final baryon asymmetry is $Y_{\rm
B}\simeq-(28/79)\kappa\sum_{i}\eta_{i}Y_{\widetilde{N}_{i}}$, where
$Y_{\widetilde{N}_{i}}$ denotes the RHN number density divided by the cosmic
entropy density and $\kappa$ represents the washout factor [46]. In our case,
for relatively low reheating temperature $T_{\rm R}$ as seen from (23), the
initial RHN abundance can be highly suppressed compared with thermal
abundance, which can lead to small $Y_{\widetilde{N}_{i}}$ for reproducing the
observed baryon asymmetry $Y_{\rm B}\sim 10^{-10}$. One should also be careful
that the cutoff scale of the theory is $\Lambda\sim M_{\rm H}/\epsilon$, and
the temperature should be lower than it for reliable calculation. We leave
detailed analyses of this issue for a future work.

### Comparison with conventional U(1)B-L

One may introduce a $Z_{2}$ parity in the conventional U(1)B-L model in which
all three RHNs have ${\rm B-L}$ charge $-1$ [28] (see e.g. Refs. [47, 48, 26,
49, 50, 51] for phenomenology of such a model). For example, let us assume
only $N_{1}$ is $Z_{2}$ odd. Then $N_{1}$ is decoupled from the SM yukawa
sector and becomes stable, and hence it is a DM candidate. The seesaw
mechanism works with two RHNs: $N_{2}$ and $N_{3}$. In this case, there does
not appear a special mass structure for three RHNs. Naturally all the three
RHNs should have masses of the same order unless some additional flavor
symmetry is introduced. On the other hand, in an alternative U(1)B-L with
$Z_{2}$, due to the special charge assignments of ${\rm B-L}$, there appears a
special RHN mass structure: two heavy degenerate RHNs and one light RHN. It is
a characteristic feature of minimal alternative U(1)B-L that is distinguished
from conventional models.

### Constraints

So far we have not discussed the size of U(1)B-L gauge coupling. If it is not
very small, the ${\rm B-L}$ gauge boson mass is not far from $v_{\phi}$ and it
is heavier than $\sim{\rm TeV}$ in the range of Fig. 1. The LHC constraint as
well as other experimental constrains are safely avoided for such a heavy
gauge boson [47, 48, 49]. If the gauge coupling is very small, on the other
hand, the ${\rm B-L}$ gauge boson becomes light and there are several
constraints from collider experiments and stellar physics [26, 50], which
would open up a possibility that the signal of ${\rm B-L}$ gauge boson will be
detected in future.

### Additional operator

We started with the Lagrangian (1) and (2) and integrated out RHNs to obtain
active neutrino masses. However, one can introduce the following operator by
hand before the integration of RHNs,

$\displaystyle\mathcal{L}_{\nu}=\sum_{\alpha,\beta=e,\mu,\tau}c_{\alpha\beta}\frac{\phi^{*2}}{\Lambda^{3}}\left(\overline{L}_{\alpha}\widetilde{H}\right)\left(\overline{L}_{\beta}\widetilde{H}\right)+{\rm
h.c.},$ (26)

where $c_{\alpha\beta}$ are $\mathcal{O}(1)$ coefficients. It leads to
neutrino masses of $m_{\nu}\sim\epsilon^{3}v_{H}^{2}/v_{\phi}$. On the other
hand, the neutrino mass that arises after the integration of RHNs is
$m_{\nu}\sim\epsilon^{-1}v_{H}^{2}/v_{\phi}$ for the heaviest one, if there is
no $Z_{2}$ symmetry ($N_{2}\leftrightarrow N_{3}$) as shown in Sec. 2. Thus
the latter one is dominant and the discussion that we need to introduce the
$Z_{2}$ symmetry to suppress the latter contribution remains intact. Once the
$Z_{2}$ symmetry is introduced, the neutrino mass that arises after
integrating out RHNs is $m_{\nu}\sim\epsilon^{9}v_{H}^{2}/v_{\phi}$ as shown
in Sec. 3, which is much suppressed compared with the contribution (26). It is
phenomenologically viable, although the relation between RHN masses and active
neutrino masses is modified from the one that we discussed in the main text.
In this case we have

$\displaystyle m_{\nu}\sim\epsilon^{3}\frac{v_{H}^{2}}{v_{\phi}}\sim 0.2\,{\rm
eV}\left(\frac{M_{\rm L}}{1\,{\rm GeV}}\right)^{3/7}\left(\frac{10^{10}\,{\rm
GeV}}{v_{\phi}}\right)^{10/7},$ (27)

instead of (22). Although the prediction that the lightest active neutrino is
massless is lost, the unique structure of the theory, i.e. the two heavy
degenerate RHNs and one light stable RHN, is unchanged.

On the other hand, it is also possible that the contribution from the operator
like (26) is subdominant. For example, let us introduce an approximate global
$Z_{3}$ symmetry under which only the ${\rm B-L}$ Higgs field transforms like
$\phi\to e^{2\pi i/3}\phi$ and all other fields are neutral. The yukawa terms
(1) is allowed by this $Z_{3}$ symmetry, but the RHN mass term (2) as well as
the operator (26) violate it. Thus they should be suppressed by some small
parameter $\eta$ that represents the violation of $Z_{3}$ symmetry:

$\displaystyle\mathcal{L}_{N}=-\sum_{i,j=1,2,3}\frac{\eta}{2}M_{ij}N_{i}N_{j}+{\rm
h.c.},$ (28)
$\displaystyle\mathcal{L}_{\nu}=\sum_{\alpha,\beta=e,\mu,\tau}\eta
c_{\alpha\beta}\frac{\phi^{*2}}{\Lambda^{3}}\left(\overline{L}_{\alpha}\widetilde{H}\right)\left(\overline{L}_{\beta}\widetilde{H}\right)+{\rm
h.c.},$ (29)

where $M_{ij}$ is given by (12). The contribution to the neutrino masses from
the operator (29) is $m_{\nu}\sim\eta\epsilon^{3}v_{H}^{2}/v_{\phi}$. The
neutrino masses that arise after integrating out RHNs are calculated in the
same way as Sec. 3 just by multiplying RHN masses by a factor $\eta$:
$m_{\nu}\sim\epsilon^{9}v_{H}^{2}/(\eta v_{\phi})$. Therefore, the latter
contribution is dominant if $\eta\lesssim\epsilon^{3}$. In such a case, the
analysis of Sec. 3 remains almost the same: the only change is to reinterpret
$v_{\phi}$ appearing in Sec. 3 as $\widetilde{v}_{\phi}\equiv\eta v_{\phi}$.
For example, Eq. (22) is unchanged once $v_{\phi}$ is interpreted as
$\widetilde{v}_{\phi}$.#3#3#3 Note that $v_{\phi}$ in Eqs. (23) and (24)
should not be reinterpreted as $\widetilde{v}_{\phi}$.

## Acknowledgments

This work was supported by JSPS KAKENHI Grant (Nos. JP19J13812 [KA], 18K03609
[KN], 17H06359 [KN] and 20J22214 [ST]).

## References

  * [1] T. Yanagida, Conf. Proc. C 7902131 (1979), 95-99
  * [2] M. Gell-Mann, P. Ramond and R. Slansky, Conf. Proc. C 790927 (1979), 315-321 [arXiv:1306.4669 [hep-th]].
  * [3] P. Minkowski, Phys. Lett. B 67 (1977), 421-428
  * [4] M. Fukugita and T. Yanagida, Phys. Lett. B 174 (1986), 45-47
  * [5] J. C. Montero and V. Pleitez, Phys. Lett. B 675 (2009), 64-68 [arXiv:0706.0473 [hep-ph]].
  * [6] B. L. Sánchez-Vega and E. R. Schmitz, Phys. Rev. D 92 (2015), 053007 [arXiv:1505.03595 [hep-ph]].
  * [7] E. Ma, N. Pollard, R. Srivastava and M. Zakeri, Phys. Lett. B 750 (2015), 135-138 [arXiv:1507.03943 [hep-ph]].
  * [8] T. Nomura and H. Okada, Eur. Phys. J. C 78 (2018) no.3, 189 [arXiv:1708.08737 [hep-ph]].
  * [9] C. Q. Geng and H. Okada, Phys. Dark Univ. 20 (2018), 13-19 [arXiv:1710.09536 [hep-ph]].
  * [10] S. Singirala, R. Mohanta, S. Patra and S. Rao, JCAP 11 (2018), 026 [arXiv:1710.05775 [hep-ph]].
  * [11] A. Das, N. Okada and D. Raut, Eur. Phys. J. C 78 (2018) no.9, 696 [arXiv:1711.09896 [hep-ph]].
  * [12] N. Okada, S. Okada and D. Raut, Phys. Rev. D 100 (2019) no.3, 035022 [arXiv:1811.11927 [hep-ph]].
  * [13] A. Das, N. Okada, S. Okada and D. Raut, Phys. Lett. B 797 (2019), 134849 [arXiv:1812.11931 [hep-ph]].
  * [14] A. Das, P. S. B. Dev and N. Okada, Phys. Lett. B 799 (2019), 135052 [arXiv:1906.04132 [hep-ph]].
  * [15] S. Mahapatra, N. Narendra and N. Sahu, [arXiv:2002.07000 [hep-ph]].
  * [16] D. Choudhury, K. Deka, T. Mandal and S. Sadhukhan, JHEP 06 (2020), 111 [arXiv:2002.02349 [hep-ph]].
  * [17] O. Lebedev, Eur. Phys. J. C 72 (2012), 2058 [arXiv:1203.0156 [hep-ph]].
  * [18] J. Elias-Miro, J. R. Espinosa, G. F. Giudice, H. M. Lee and A. Strumia, JHEP 06 (2012), 031 [arXiv:1203.0237 [hep-ph]].
  * [19] P. H. Frampton, S. L. Glashow and T. Yanagida, Phys. Lett. B 548 (2002), 119-121 [arXiv:hep-ph/0208157 [hep-ph]].
  * [20] A. Ibarra and G. G. Ross, Phys. Lett. B 591 (2004), 285-296 [arXiv:hep-ph/0312138 [hep-ph]].
  * [21] K. Nakayama, F. Takahashi and T. T. Yanagida, Phys. Lett. B 757 (2016), 32-38 [arXiv:1601.00192 [hep-ph]].
  * [22] K. Nakayama, F. Takahashi and T. T. Yanagida, Phys. Lett. B 773 (2017), 179-185 [arXiv:1705.04796 [hep-ph]].
  * [23] P. A. Zyla et al. [Particle Data Group], PTEP 2020 (2020) no.8, 083C01
  * [24] S. Khalil and O. Seto, JCAP 10 (2008), 024 [arXiv:0804.0336 [hep-ph]].
  * [25] A. Kusenko, F. Takahashi and T. T. Yanagida, Phys. Lett. B 693 (2010), 144-148 [arXiv:1006.1731 [hep-ph]].
  * [26] K. Kaneta, Z. Kang and H. S. Lee, JHEP 02 (2017), 031 [arXiv:1606.09317 [hep-ph]].
  * [27] A. Biswas and A. Gupta, JCAP 09 (2016), 044 [arXiv:1607.01469 [hep-ph]].
  * [28] N. Okada and O. Seto, Phys. Rev. D 82 (2010), 023507 [arXiv:1002.2525 [hep-ph]].
  * [29] M. Shaposhnikov and I. Tkachev, Phys. Lett. B 639 (2006), 414-417 [arXiv:hep-ph/0604236 [hep-ph]].
  * [30] A. Kusenko, Phys. Rev. Lett. 97 (2006), 241301 [arXiv:hep-ph/0609081 [hep-ph]].
  * [31] K. Petraki and A. Kusenko, Phys. Rev. D 77 (2008), 065014 [arXiv:0711.4646 [hep-ph]].
  * [32] A. Merle, V. Niro and D. Schmidt, JCAP 03 (2014), 028 [arXiv:1306.3996 [hep-ph]].
  * [33] L. Covi and E. Roulet, Phys. Lett. B 399, 113 (1997) [hep-ph/9611425].
  * [34] A. Pilaftsis, Nucl. Phys. B 504, 61 (1997) [hep-ph/9702393].
  * [35] W. Buchmuller and M. Plumacher, Phys. Lett. B 431, 354 (1998) [hep-ph/9710460].
  * [36] A. Pilaftsis and T. E. J. Underwood, Nucl. Phys. B 692, 303 (2004) [hep-ph/0309342].
  * [37] A. Pilaftsis and T. E. J. Underwood, Phys. Rev. D 72, 113001 (2005) [hep-ph/0506107].
  * [38] A. Anisimov, A. Broncano and M. Plumacher, Nucl. Phys. B 737, 176 (2006) [hep-ph/0511248].
  * [39] B. Dev, M. Garny, J. Klaric, P. Millington and D. Teresi, Int. J. Mod. Phys. A 33 (2018), 1842003 [arXiv:1711.02863 [hep-ph]].
  * [40] M. Garny, A. Kartavtsev and A. Hohenegger, Annals Phys. 328 (2013), 26-63 [arXiv:1112.6428 [hep-ph]].
  * [41] B. Garbrecht and M. Herranen, Nucl. Phys. B 861, 17 (2012) [arXiv:1112.5954 [hep-ph]].
  * [42] S. Iso, K. Shimada and M. Yamanaka, JHEP 1404, 062 (2014) [arXiv:1312.7680 [hep-ph]].
  * [43] S. Iso and K. Shimada, JHEP 1408 (2014) 043 [arXiv:1404.4816 [hep-ph]].
  * [44] B. Garbrecht, F. Gautier and J. Klaric, JCAP 1409, no. 09, 033 (2014) [arXiv:1406.4190 [hep-ph]].
  * [45] V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, Phys. Lett. B 155 (1985), 36
  * [46] G. F. Giudice, A. Notari, M. Raidal, A. Riotto and A. Strumia, Nucl. Phys. B 685 (2004), 89-149 [arXiv:hep-ph/0310123 [hep-ph]].
  * [47] N. Okada and S. Okada, Phys. Rev. D 93 (2016) no.7, 075003 [arXiv:1601.07526 [hep-ph]].
  * [48] N. Okada and S. Okada, Phys. Rev. D 95 (2017) no.3, 035025 [arXiv:1611.02672 [hep-ph]].
  * [49] S. Okada, Adv. High Energy Phys. 2018 (2018), 5340935 [arXiv:1803.06793 [hep-ph]].
  * [50] N. Okada, S. Okada and Q. Shafi, Phys. Lett. B 810 (2020), 135845 [arXiv:2003.02667 [hep-ph]].
  * [51] D. Borah, S. Jyoti Das and A. K. Saha, [arXiv:2005.11328 [hep-ph]].


# Ada-Detector: Adaptive Frontier Detector for Rapid Exploration

Zezhou Sun†, Banghe Wu†, Chengzhong Xu‡, Hui Kong∗ ${\dagger}$ School of
Computer Science and Engineering, Nanjing University of Science and
Technology, Nanjing, Jiangsu, China.${\ddagger}$ The State Key Laboratory of
Internet of Things for Smart City (SKL-IOTSC), Department of Computer Science,
University of Macau, Macau, China.$*$ The State Key Laboratory of Internet of
Things for Smart City (SKL-IOTSC), Department of Electromechanical Engineering
(EME), University of Macau, Macau, China. This work is supported by National
Key Research and Development Program of China (No. 2019YFB2102100), the
Science and Technology Development Fund of Macau SAR (File no. 0015/2019/AKP),
Guangdong-Hong Kong-Macao Joint Laboratory of Human-Machine Intelligence-
Synergy Systems (No. 2019B121205007), and the National Natural Science
Foundation of China (No.61803083).

###### Abstract

In this paper, we propose an efficient frontier detector method based on
adaptive Rapidly-exploring Random Tree (RRT) for autonomous robot exploration.
Robots can achieve real-time incremental frontier detection when they are
exploring unknown environments. First, our detector adaptively adjusts the
sampling space of RRT by sensing the surrounding environment structure. The
adaptive sampling space can greatly improve the successful sampling rate of
RRT (the ratio of the number of samples successfully added to the RRT tree to
the number of sampling attempts) according to the environment structure and
control the expansion bias of the RRT. Second, by generating non-uniform
distributed samples, our method also solves the over-sampling problem of RRT
in the sliding windows, where uniform random sampling causes over-sampling in
the overlap area between two adjacent sliding windows. In this way, our
detector is more inclined to sample in the latest explored area, which
improves the efficiency of frontier detection and achieves incremental
detection. We validated our method in three simulated benchmark scenarios. The
experimental comparison shows that we reduce the frontier detection runtime by
about $40\%$ compared with the SOTA method, DSV Planner.

##  I INTRODUCTION

The frontier detection module is a key module for mobile robots to
autonomously explore unknown environments. Frontiers are the boundaries
separating known space from unknown area. The frontier detection algorithms
detect the frontiers according to the map established in real-time during the
movement of the robot and send the frontier with the highest priority to the
path planning module. The ideal frontier detection module deployed to an
autonomous exploratory robot is expected to detect new frontiers within the
interval between two LiDAR scans and empowers robots to replace humans to
achieve efficient and robust autonomous exploration in various environments.

Since the RRT algorithm [1] is heavily biased to grow towards unknown regions
of the map [2] and can easily be extended to high-dimensional spaces, it is
widely used for frontier detection tasks in the field of autonomous
exploration. The original RRT-based frontier detector [3] extends an RRT in
the free space (the space that has been observed by sensors and is not
occupied by obstacles) and treats the nodes in the RRT as viewpoints. The
exploration gain of each viewpoint is calculated according to the sensor
coverage. The viewpoints whose exploration gain is greater than the preset
threshold are called frontiers. Thereafter, the robot sorts the frontiers
according to the exploration gain and drives to the frontier with the highest
priority.

Compared with path planning, since frontier detection does not have a clear
expansion target, it can only use the RRT method to expand to the surrounding
area uniformly instead of using RRT improvement methods such as the Informed
RRT* [4], Batch Informed trees (BIT*) [5], and Regionally accelerated batch
informed trees (RABIT*) [6] etc. Therefore, the problem of how to improve the
efficiency of frontier detection has attracted the attention of researchers.
The Graph-Based subterranean exploration path Planner (GBPlanner) [7] uses a
periodic sliding window to restrict the sampling space of RRT to the
surrounding region of the robot. Compared to covering the entire exploration
area, RRT can fill the sliding window more quickly. Thus, this method
significantly improves the efficiency of frontier detection.

Figure 1:  An illustration of the dynamically expanded frontier detector in
DSV Planner [8]. (a) shows the RRT-based frontier detection performed in a
sliding window. (b) After the sliding window is updated with the movement of
the robot, DSV Planner [8] discards the node that is beyond the updated
sliding window and continues expanding the RRT tree. (c) In two adjacent
sliding windows, the overlap area is repeatedly sampled twice and the non-
overlap area is the new observation area.

Usually, the size of the sliding window is set as the sensing range limit of
the sensor on the robot. However, with the improvement of sensor manufacturing
technology, the sensing range limit of a modern multi-channel LiDAR has
exceeded 100m. Thus, the sliding window bounded by the sensing range limit is
too large for real-time frontier detection. In addition, the use of a fixed-
sized sliding window cannot fully reflect the robot’s perception of the
surrounding structure of the environment. For example, when passing through a
narrow corridor, it is preferred that the RRT should expand in the direction
of the corridor. When passing an intersection, it should pay more attention to
the left and right sides. When entering a room, the sampling space of RRT
should be limited by walls. In contrast, the fixed-sized sliding window
enables the RRT to expand uniformly in all directions. Through the sensor’s
perception of the surrounding environment structure, we can adaptively reduce
the size of the sliding window to the smallest circumscribed rectangle of
LiDAR scan.

Furthermore, there is inevitably an overlap area between two adjacent sliding
windows. The non-overlap area of the current sliding window actually
corresponds to the new observed (explored) area of the LiDAR sensor. The
GBPlanner [7] uniformly generates sampling points in each sliding window.
Thus, the overlap areas are twice as dense as the non-overlap area due to
repeated sampling. We call this phenomenon an over-sampling problem, as shown
in Fig. 1(c). However, it is preferred that the detector should focus more on
detecting frontiers in the new observed area. Obviously, the over-sampling
problem can reduce the tendency of extending the RRT tree to a new observed
area.

Very recently, the Dual-Stage Viewpoint Planner (DSV Planner) [8] proposed a
dynamically extended frontier detector. Figure 1 shows the dynamic expansion
process. In Fig. 1(a), the blue box is the sliding window, representing the
sampling space of RRT. It expands an RRT tree to detect the frontiers in the
sliding window. When the robot moves to the new position in Fig. 1(b), the
sliding window is moved accordingly to the new position of the robot (the blue
box in (b)). Nodes in the RRT tree that is beyond the sliding window are
discarded. Iteratively, the robot continues expanding the RRT tree in the
updated sliding window. Once detecting the frontiers, DSV Planner [8] selects
the frontiers with the highest priority and additionally increases the
sampling probability within the neighborhood of these frontiers. When the
number of nodes in the RRT tree exceeds a certain threshold or the number of
sampling attempts reaches the upper limit, the robot considers that the
frontier detection in this sliding window is completed. Then the detector
sends the frontier with the highest priority to the path planning module.

However, the DSV Planner [8] still makes lots of sampling attempts in the
overlap area, which is hardly helpful for the robot to travel to the newly
explored frontiers. In addition, this method only increases the bias of the
RRT expansion toward some specific frontiers, causing the detector to pay
insufficient attention to the newly explored area that is beyond the
neighborhood of selected frontiers. It can also make RRT trapped into obstacle
regions in a complex environment. To deal with these issues, we counteract
over-sampling and achieve incremental detection by generating non-uniform
samples.

Figure 2:  An illustration of the improvement of our method from the previous
method. Compared with the previous methods that detect frontiers in fixed-
sized sliding windows, our method adaptively adjusts the sliding windows
through the perception of the surrounding structural environment. Our detector
also generates non-uniform distributed samples to solve the over-sampling
problem in the overlap area between two adjacent sliding windows.

The main contributions of this paper are as follows. (1) We propose an
adaptive frontier detector with a varying-sized sliding window, which can
improve the efficiency of frontier detection by limiting the sampling space of
RRT. (2) We generate non-uniform samples to avoid over-sampling between two
adjacent sliding windows, making RRT biased to detect frontiers in newly
explored areas and achieve incremental detection.

The remainder of this paper is organized as follows. Section II outlines the
related works. The details of our method are given in Section III.
Experimental results are shown in Section IV, followed by conclusions in
Section V.

##  II Related works

Frontier detection divides the created occupancy map into free area, unknown
area, and occupied area according to the confidence and searches for the
boundaries between free area and unknown area in the map [9]. According to two
different path planning algorithms, Dijkstra [10] and RRT, frontier detection
is mainly divided into two categories: search-based and sampling-based.

####  II-1 The search-based frontier detectors

Yamauchi et al. [9] proposed their seminal search-based exploration work based
on frontier detection and tracking. In order to detect frontiers more
efficiently, the following-up works are devoted to reducing the search space
of frontiers. The Wavefront Frontier Detection (WFD) [11] reduces the search
space from the entire map to the known area. The Expanding-Wavefront Frontier
Detection (EWFD) [12] only incrementally searches for the latest explored
areas. The Dense Frontier Detection method (DFD) [13] performs frontier
detection in the submaps of Cartographer [14], reducing the search space to
the total frontier length of all previous submaps, so that the detected
frontiers are not affected by the graph optimization to the map. The improved
DFD method based on the sliding window [15] reduce the search space to the
frontier pool of the latest modified submaps. These methods are designed to
detect all frontiers in the exploration area, however, they are mostly used
for 2D exploration and are difficult to extend to high-dimensional spaces.

####  II-2 The sampling-based frontier detectors

The sampling-based frontier detection methods are more suitable for robots to
explore in the 3D environment. Since the “next best view” method was proposed
[16], Oriolo [17] and Freda [18] introduced a probabilistic planning method
called Sensor-based Random Tree (SRT), a variant of RRT, to effectively sample
the map. The Multiple Rapidly-exploring Randomized Trees (multi-layer RRT)[2]
improves the efficiency of frontier detection by using global and local RRT
trees. With more applications of unmanned aerial vehicles (UAV), RRT-based
frontier detection methods have been widely used in UAV exploration, such as
the Receding Horizon Next-Best-View method (NBVP) [3] and its variants [19]
[20]. However, the sampling space of these methods is the explored map that
expands with the progress of robot exploration, which means that the
efficiency of frontier detection continues to decrease.

Recently, the Graph-Based exploration path Planner (GBPlanner) [7] and Motion-
primitive Based exploration path Planner (MBPlanner) [21] implement a fast
frontier detector based on periodic sliding windows. It limits the sampling
space from the free space of all explored areas to that of the local map
around the robot within the LiDAR sensing range. However, GBPlanner‘s detector
[7] cannot instantly detect all frontiers within the LiDAR sensing range
during the interval of the periodic sliding window. Therefore, robots can miss
partial frontiers and leave some unexplored areas to be ignored.

On this basis, the Dual-Stage Viewpoint Planner (DSV Planner) [8] proposed an
efficient exploration algorithm based on a dynamically extended frontier
detector. DSV Planner[8] achieves SOTA frontier detection efficiency by
reusing part of the nodes in the previous sliding window and additionally
expanding toward specific frontiers.

For a fair comparison of different algorithms, Zhang et al. [22] established a
public simulation environment to evaluate exploration methods and provided the
best human practice results [22]. The simulation environment provides 5
different scenarios including multi-layer garage, indoor, forest, tunnel, and
campus environments. Zhang et al. tested their DSV Planner [8] and TARE [23]
with NBVP [3], GBPlanner [7], MBPlanner [21], and the best human practice
results [22] in these simulation environments.

##  III Methodology

Figure 3: The system diagram of our Ada-detector.

###  III-A Adaptive Sliding Window

In the simulation environment, since it is necessary to detect slopes so that
the robot can explore in 3D space, the maximum range of LiDAR sensing that can
be used for mapping is set to $15$m. Therefore, as did in DSV Planner [8], we
do not consider the LiDAR points beyond this limit. Specifically, we first
align the filtered LiDAR scan with the map coordinates. Then, by traversing
the LiDAR points, we calculate the four corners of the smallest circumscribed
rectangle of the LiDAR scan. The rectangle formed by the corners is the
current sliding window and the RRT regards it as the sampling space. We update
the sliding window by updating the corners.

Generally, the successful sampling rate (the ratio of the sampled points that
can be successfully added to the RRT tree to the number of sampling attempts)
is directly related to the proportion of free space to the total sampling
space. It can be seen from Fig. 4 that the successful sampling rate is
generally much higher in an adaptive sliding window than in a fixed-sized one.

Figure 4:  (a) shows the sampling space of a fixed-sized sliding window and
(b) shows the sampling space of an adaptive sliding window of our method. With
the same number of sampling attempts (green points), our method has a higher
successful sampling rate.

###  III-B Non-uniform Sampling within an Adaptive Sliding Window

Figure 5:  An illustration of non-uniform sampling process. (a) Just after the
sliding window is updated, the nodes that exceed the sliding window in the RRT
tree are discarded and the remaining $N_{o}$ nodes are all located in the
overlap area. (b) According to Alg. 1, we calculate how many points are needed
to be sampled in the overlap area and the non-overlap area, respectively, to
make the density of RRT nodes in the two areas equal and both larger than
$\theta$. We then perform a non-uniform sampling based on this. (c) We expand
the RRT tree based on the non-uniform sampling.

To deal with the over-sampling problem of DSV Planner [8] as mentioned before,
we give the non-uniform sampling strategy. In each new adaptive sliding
window, let $\theta$ be the ratio of the number of successful samples (samples
successfully added to the RRT tree) to the size of the sampling space, called
the successful sampling density. Let $N_{total}$ be the number of RRT tree
nodes, $S_{o}$ and $S_{n}$ be the area of the overlap and non-overlap region,
respectively, $N_{o}$ and $N_{n}$ be the number of RRT nodes in the overlap
and non-overlap region, respectively. The motivation of the proposed non-
uniform sampling is that, for an efficient detector, the value of $\theta$
should be the same no matter in the non-overlap or overlap area in each
adaptive sliding window, i.e.,
$\frac{N_{o}}{S_{o}}=\frac{N_{n}}{S_{n}}=\theta$. Let $\tau$ be the maximal
number of successful samples in the whole adaptive sliding window, and thus we
have $\theta=\frac{\tau}{S_{total}}$, where $S_{total}$ is the area of the
sliding window. Therefore, once the condition
$\theta=\frac{\tau}{S_{total}}=\frac{N_{o}}{S_{o}}=\frac{N_{n}}{S_{n}}$ is
satisfied, it is guaranteed to solve the over-sampling problem in the DSV
Planner [8].

Note that $\tau$ is a constant for each adaptive sliding window, but
$S_{total}$, $S_{o}$ and $S_{n}$ are not, meaning that $N_{o}$ and $N_{n}$ are
both variables. Generally, $S_{o}$ and $S_{n}$ are different. Therefore,
$N_{o}$ and $N_{n}$ are not equal as well. This explains why we need a non-
uniform sampling within each new adaptive sliding window, where we treat the
sampling in the overlap and non-overlap region differently.

Specifically, given the area of a new updated sliding window $S_{total}$, we
can obtain $\theta$. We can get the overlap area $S_{o}$ and non-overlap area
$S_{n}$ based on the four corners of the previous and the current sliding
window, respectively, (lines 2-3 of Alg. 1). In turn, we can get the number of
RRT nodes that should be successfully sampled from the overlap area, $N_{eo}$,
and from the non-overlap area, $N_{en}$, respectively, (lines 4 of Alg.1), as
shown in Fig. 5. Therefore, we can get the sampling probability $\eta_{n}$ of
non-overlap area and $\eta_{o}$ of overlap area in the current sliding window
to be $\frac{N_{en}}{N_{eo}+N_{en}}$ and $\frac{N_{eo}}{N_{eo}+N_{en}}$ (line
4-5 of Alg.1), respectively . Once again, when the number of RRT nodes in the
current sliding window equals $\tau$ or the successful sampling density is
larger than $\theta$, we update the sliding window. Meanwhile, we discard the
nodes in the RRT tree that are beyond the new sliding window and reserve the
ones in the overlap area (line 1 of Alg.1).

Input:

corners of the previous sliding window $Pre_{corners}$

corners of the new sliding window $P_{corners}$

list of RRT tree nodes $RRT_{nodes}$

expected sampling density $\theta$

Output:

sampling probability in the non-overlap area $\eta_{n}$

sampling probability in the overlap area $\eta_{o}$

1

2($RRT_{nodes}$, $N_{o}$) $\leftarrow$ Discard($RRT_{nodes}$,$P_{corners}$)

3 ($S_{o}$, $S_{n}$) $\leftarrow$ overlapArea($Pres_{corners}$, $P_{corners}$)

4 $Pre_{corners}$ = $P_{corners}$

5 $N_{eo}$ = $\theta*S_{o}-N_{o}$; $N_{en}$ = $\theta*S_{n}$

6 $\eta_{n}=\frac{N_{en}}{N_{eo}+N_{en}}$;
$\eta_{o}=\frac{N_{eo}}{N_{eo}+N_{en}}$

Algorithm 1 Non-uniform sampling within an adaptive sliding window

###  III-C Frontier Detection

Given $\eta_{n}$ and $\eta_{o}$, we expand RRT in the current adaptive sliding
window in a non-uniform sampling mode (Alg.2). We randomly generate a number
$r$, with $r\in[0,1]$. If $r>\eta_{n}$, we sample in the overlap area.
Instead, we sample in the non-overlap area (lines 3-6). After obtaining the
sample $P_{rand}$, we expand the RRT tree. Lines 7-11 of Alg. 2 are the
standard RRT process [1].

If the sample can be successfully added to the RRT tree, we calculate the
exploration gain of this node and increase the number of successful samples by
one (line 12-14). The method of calculating exploration gain is the same as
the one used in DSV Planner [8]. Since we want to keep the successful sampling
density greater than or equal to $\theta$ in each adaptive sliding window.
Once an adaptive sliding window is updated, the successful sampling density of
the overlap area of the new adaptive sliding window is approximately equal to
$\theta$. From this point of view, our method is about equivalent to an
incremental frontier detection.

Input:

probability of sampling in non-overlap area $\eta_{n}$

list of RRT tree nodes after discarding $RRT_{nodes}$

the adaptive sliding window area $S_{total}$

number of RRT tree nodes after discarding $N_{o}$

occupancy map $M$

min exploration gain $MinExplorationGain$

expected sampling density $\theta$

the upper limit of the RRT nodes in each sliding window $\tau$

Output:

detected frontiers $P_{frontier}$

1 $N_{total}$ = 0

2 while  _$N_{total} <\tau$ and $N_{total}/S_{total}<\theta$ _ do

3 if  _rand() $>\eta_{n}$_ then

4 $P_{rand}$ $\leftarrow$ Sample in overlap area

5 else

6 $P_{rand}$ $\leftarrow$ Sample in non-overlap area

7 $P_{nearest}$ = nearest($RRT_{nodes}$, $P_{rand}$)

8 $P_{steer}$ = steer($P_{nearest}$, $P_{rand}$)

9 if  _checkingObstacle( $P_{nearest}$, $P_{steer}$, $M$)_ then

10 $RRT_{nodes}$.push_back($P_{steer}$)

11 $N_{total}$ = $N_{total}$ \+ 1

12 if  _explorationGain( $P_{steer}$, $M$) $>MinExplorationGain$_ then

13 $P_{frontier}\leftarrow P_{steer}$

Algorithm 2 Frontier Detection

###  III-D Time Complexity Comparison with DSV Planner [8]

If $N$ points are successfully sampled in the non-overlap area, we suppose
$k*N$ successful samples are needed in the whole sliding window. We suppose a
successful sampling in a sliding window requires $p$ sampling attempts.
Therefore, the total number of required sampling attempts is $pk*N$. The time
cost of expanding $N$ nodes in the non-overlap area (regardless of the
initialization step) can be calculated as $pk*N$ times the cost of a single
while-loop in Alg. 2.

$\displaystyle T(pk*N)=T_{sample}(pk*N)+T_{nearest}(pk*N)+$ (1) $\displaystyle
T_{steer}(pk*N)+T_{checkingObstacle}(pk*N)+T_{add}(k*N)$

where $T_{sample}$, $T_{steer}$ are simple operations in linear time with a
complexity of $O(pk*N)$. According to [24], the complexity of finding the
nearest vertex in a kd-tree and adding a vertex to the tree is
$O(n^{1-\frac{1}{d}})$ and $O(\log n)$, respectively, where $d$ represents the
search dimension. Therefore, the complexity of $T_{nearest}$ and $T_{add}$ is
about $O((pk*N)^{\frac{3}{2}})$ and $O((k*N)\log(k*N))$, respectively.
$T_{add}$ only needs to be executed $k*N$ times because the remaining
unsuccessful sampling attempts are discarded in the “checkingObstacle” step.
$T_{checkingObstacle}$ is related to the “steer” step, which takes constant
time. Thus, the complexity of $T_{checkingObstacle}$ is $O(pk*N)$. In
addition, the preprocessing of the LiDAR point cloud takes constant time and
can be easily integrated into the point cloud processing module, and is
ignored. Therefore, the overall complexity of our approach is

$\displaystyle O(pk*N)+O((pk*N)^{\frac{3}{2}})+O(pk*N)+O(pk*N)+$ (2)
$\displaystyle
O((k*N)\log(k*N))\approx(pk)^{\frac{3}{2}}*O(N^{\frac{3}{2}})+3pk*O(N)+$
$\displaystyle(k*N)*O(\log(k*N)).$

Since $N$ is not extremely large, the first-order items of $N$ are not
neglectable, we thus keep them. Note that the complexity of the DSV Planner
[8] can also represented as above. But $p$ and $k$ are different in our method
and the DSV Planner. Generally, the successful sampling rate is proportional
to the ratio of the free space to the total space in the sliding window. From
Section III.A, it is known that the sliding windows of both our method and the
DSV Planner cover the same free space (the space that has been observed by
sensors and is not occupied by obstacles). Let the sliding window area of our
method and the DSV Planner is $S_{total}$ and $S_{dsv}$, respectively, then
the ratio of the successful sampling rate of our method to that of the DSV
Planner [8] is equal to $S_{dsv}$/$S_{total}$. When $N$ points are
successfully sampled in the non-overlapping area, our method needs to sample
$\frac{N}{\eta_{n}}$ points successfully in the whole sliding window, while
DSV Planner [8] needs to sample $0.2*N+0.8*N*\frac{S_{total}}{S_{n}}$ points,
where $0.2$ and $0.8$ is the bias sampling parameter of DSV Planner [8].

Based on the analysis above, the $p$ value is significantly larger in the DSV
Planner [8] than the one in our method, especially in narrow passageways. In
wide open regions with obstacles, our method is still more favorable with
quite much improved efficiency. In the worst case where there is no obstacle
around, the adaptive sliding window method is equivalent to the fixed-sized
sliding window scheme. For $k$, we mainly save the calculation of sampling in
overlap and redundant regions. Our method almost never performs sampling in
overlap or redundant area that beyond the current LiDAR scan range, as shown
in Fig. 2(a).

##  IV Experiments

We obtain our evaluation results using the default settings of the Autonomous
Exploration Development Environment [22], which are the same as in the DSV
Planner [8]. We compare our method with DSV Planner [8] in three scenarios:
the indoor corridors, forest, and multi-storage garage scenes 111A
representative exploration result of the indoor corridors environment can be
found at: https://youtu.be/jEY3mSrxRoU. Our algorithm’s run-time is evaluated
based on a 2.6GHz i7-9750H CPU. We evaluate the average exploration efficiency
$\varepsilon$, the average exploration volume $V$, the average exploration
distance $L$, the average sampling attempts in each sliding window, the
average successful samples (the samples successfully added to the RRT tree) in
each sliding window, and the average sliding window duration of two methods in
Table. I.

() Indoor Corridors Environment

() Forest Environment

() Multi-storage Garage Environment

Figure 6: In (a), upper left: the view of the indoor corridor environment.
upper right: the chart of instant exploration efficiency of two methods, where
the vertical axis is the exploration space volume and the horizontal axis is
the exploration time. lower left: the chart of instant exploration distance of
two methods, where the vertical axis is the traveling distance and the
horizontal axis is the exploration time. lower right: the chart of instant
algorithm runtime (including the frontier detection module and the path
planning module) of two methods, where the vertical axis is total running time
and the horizontal axis is the exploration time. (b) and (c) show the results
in the forest and multi-storage garage scenarios, respectively.

() Indoor Corridors Environment

() Forest Environment

() Multi-storage Garage Environment

Figure 7: In (a), upper left: the best exploration trajectory of our method in the indoor corridors simulated environments. upper right: the chart of the average sampling attempts of two methods, where the vertical axis is the number of sampling attempts in each sliding window and the horizontal axis is the sliding window update times. lower left: the chart of the average successful samples of two methods, where the vertical axis is the number of successful samples added into the RRT tree in each sliding window and the horizontal axis is the sliding window update times. lower right: the chart of the average sliding window duration of two methods, where the vertical axis is the time of sliding window duration and the horizontal axis is the sliding window update times. (b) and (c) show the results in the forest and multi-storage garage scenarios, respectively.  TABLE I: The efficiency of different methods in the simulated benchmark environments.  | DSV Planner [8] | Ada-detector
---|---|---
| $\varepsilon$ | $V$ | $L$ |  | Sampling
---
Attemps
| Successful
---
Samples
| Sliding window
---
duration
$\varepsilon$ | $V$ | $L$ |  | Sampling
---
Attemps
| Successful
---
Samples
| Sliding window
---
duration
Indoor | 6.8 | 5711 | 1384 | 89062 | 83 | 0.64 | 7.0 | 5703 | 1380 | 61943 | 41 | 0.29
Forest | 37.7 | 43647 | 2159 | 56798 | 98 | 0.47 | 37.7 | 42571 | 2124 | 41123 | 86 | 0.37
Garage | 16.4 | 43311 | 4952 | 51558 | 63 | 0.31 | 16.3 | 43216 | 4946 | 25212 | 43 | 0.18

The vehicle navigates at $2m/s$. Like the DSV Planner [8], we run our methods
10 times. The average exploration efficiency $\varepsilon(m^{3}/s)$ = mean of
(total explored volume of a run)/(total time of a run) of 10 runs. The average
exploration volume $V$ and distance $L$ are the mean of the final exploration
volume and distance of 10 runs, respectively. The exploration data
$\varepsilon,V$ and $L$ of DSV Planner [8] are provided by the author [22].
The average sampling attempts, the average successful samples, and the average
sliding window duration in each sliding window of the two methods are the mean
values of 10 runs in each sliding window, respectively. We draw them and the
best exploration trajectories of three scenes in Fig. 7.

Table. I shows that the successful sampling rate of our method in each sliding
window is $73\%$ (indoor), $116\%$ (forest), and $139\%$ (garage) of that of
the DSV Planner [8], respectively. In the indoor environment, due to dense
obstacles, the fixed-sized sliding windows of DSV Planner [8] often cover the
free space of other rooms (Fig. 7(a)), i.e., the free space covered by the
fixed-sized sliding window of DSV Planner [8] is larger than that covered by
the adaptive sliding window. Although higher in this case, these samples are
not helpful for frontier detection because they are in other rooms, still
reducing the frontier detection efficiency. In the forest and garage
environments, our method adaptively adjusts the size of the sliding window
according to the surrounding environment structure to obtain a higher
successful sampling rate.

The average number of RRT nodes added in each adaptive sliding window of three
environments is $49.4\%$ (indoor), $87.8\%$ (forest), and $68.3$ (garage) of
that of the DSV Planner [8], respectively. In three different types of
scenarios, our exploration efficiency is very close to that of DSV Planner
[8], as shown in Fig. 6. It indicates that both our method and the DSVP method
can adequately fill the respective sliding windows. We believe that for a
fully filled sliding window, both methods using the same exploration gain
reach the upper limit of exploration efficiency.

Higher successful sampling rate and fewer added RRT nodes lead to higher
frontier detection efficiency. The average duration of each sliding window is
$0.29$, $0.37$ and $0.18$ seconds, which is $45.3\%$, $78.7\%$ and $58.1\%$ of
that of DSV Planner [8] in the indoor, forest, and garage environments,
respectively. We reduce the frontier detection runtime by about $40\%$
compared with the DSV Planner [8].

##  V Conclusion

In this paper, we propose an adaptive RRT-based frontier detector for
autonomous exploration. Our main contributions include an adaptive sliding
window of RRT to improve successful sampling rate and a non-uniform sampling
strategy to solve the over-sampling problem of RRT between two adjacent
sliding windows. We validated our method in three simulated benchmark
scenarios. The experimental comparison shows that we reduce the frontier
detection runtime by about $40\%$ compared with the SOTA method.

## References

  * [1] Steven M LaValle et al.  Rapidly-exploring random trees: A new tool for path planning.  1998\.
  * [2] Hassan Umari and Shayok Mukhopadhyay.  Autonomous robotic exploration based on multiple rapidly-exploring randomized trees.  In 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 1396–1402. IEEE, 2017.
  * [3] Andreas Bircher, Mina Kamel, Kostas Alexis, Helen Oleynikova, and Roland Siegwart.  Receding horizon” next-best-view” planner for 3d exploration.  In 2016 IEEE international conference on robotics and automation (ICRA), pages 1462–1468. IEEE, 2016.
  * [4] Jonathan D Gammell, Siddhartha S Srinivasa, and Timothy D Barfoot.  Informed rrt*: Optimal sampling-based path planning focused via direct sampling of an admissible ellipsoidal heuristic.  In 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 2997–3004. IEEE, 2014.
  * [5] Jonathan D Gammell, Siddhartha S Srinivasa, and Timothy D Barfoot.  Batch informed trees (bit*): Sampling-based optimal planning via the heuristically guided search of implicit random geometric graphs.  In 2015 IEEE international conference on robotics and automation (ICRA), pages 3067–3074. IEEE, 2015.
  * [6] Sanjiban Choudhury, Jonathan D Gammell, Timothy D Barfoot, Siddhartha S Srinivasa, and Sebastian Scherer.  Regionally accelerated batch informed trees (rabit*): A framework to integrate local information into optimal path planning.  In 2016 IEEE International Conference on Robotics and Automation (ICRA), pages 4207–4214. IEEE, 2016.
  * [7] Mihir Dharmadhikari, Tung Dang, Lukas Solanka, Johannes Loje, Huan Nguyen, Nikhil Khedekar, and Kostas Alexis.  Motion primitives-based path planning for fast and agile exploration using aerial robots.  In 2020 IEEE International Conference on Robotics and Automation (ICRA), pages 179–185. IEEE, 2020.
  * [8] Hongbiao Zhu, Chao Cao, Yukun Xia, Sebastian Scherer, Ji Zhang, and Weidong Wang.  Dsvp: Dual-stage viewpoint planner for rapid exploration by dynamic expansion.  In 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). IEEE.
  * [9] Brian Yamauchi.  A frontier-based approach for autonomous exploration.  In Proceedings 1997 IEEE International Symposium on Computational Intelligence in Robotics and Automation CIRA’97.’Towards New Computational Principles for Robotics and Automation’, pages 146–151. IEEE, 1997\.
  * [10] Philip L Frana and Thomas J Misa.  An interview with edsger w. dijkstra.  Communications of the ACM, 53(8):41–47, 2010.
  * [11] Matan Keidar and Gal A Kaminka.  Robot exploration with fast frontier detection: theory and experiments.  In Proceedings of the 11th International Conference on Autonomous Agents and Multiagent Systems-Volume 1, pages 113–120. International Foundation for Autonomous Agents and Multiagent Systems, 2012.
  * [12] Phillip Quin, Alen Alempijevic, Gavin Paul, and Dikai Liu.  Expanding wavefront frontier detection: An approach for efficiently detecting frontier cells.  In Australasian Conference on Robotics and Automation, ACRA, 2014\.
  * [13] Juraj Oršulić, Damjan Miklić, and Zdenko Kovačić.  Efficient dense frontier detection for 2-d graph slam based on occupancy grid submaps.  IEEE Robotics and Automation Letters, 4(4):3569–3576, 2019.
  * [14] Wolfgang Hess, Damon Kohler, Holger Rapp, and Daniel Andor.  Real-time loop closure in 2d lidar slam.  In 2016 IEEE International Conference on Robotics and Automation (ICRA), pages 1271–1278. IEEE, 2016.
  * [15] Zezhou Sun, Banghe Wu, Cheng-Zhong Xu, Sanjay E. Sarma, Jian Yang, and Hui Kong.  Frontier detection and reachability analysis for efficient 2d graph-slam based active exploration.  In IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS 2020, Las Vegas, NV, USA, October 24, 2020 - January 24, 2021, pages 2051–2058. IEEE, 2020.
  * [16] Cl Connolly.  The determination of next best views.  In Proceedings. 1985 IEEE international conference on robotics and automation, volume 2, pages 432–435. IEEE, 1985.
  * [17] Giuseppe Oriolo, Marilena Vendittelli, Luigi Freda, and Giulio Troso.  The srt method: Randomized strategies for exploration.  In IEEE International Conference on Robotics and Automation, 2004\. Proceedings. ICRA’04. 2004, volume 5, pages 4688–4694. IEEE, 2004.
  * [18] Luigi Freda and Giuseppe Oriolo.  Frontier-based probabilistic strategies for sensor-based exploration.  In Proceedings of the 2005 IEEE International Conference on Robotics and Automation, pages 3881–3887. IEEE, 2005.
  * [19] Tung Dang, Christos Papachristos, and Kostas Alexis.  Visual saliency-aware receding horizon autonomous exploration with application to aerial robotics.  In 2018 IEEE International Conference on Robotics and Automation (ICRA), pages 2526–2533. IEEE, 2018.
  * [20] Christos Papachristos, Shehryar Khattak, and Kostas Alexis.  Uncertainty-aware receding horizon exploration and mapping using aerial robots.  In 2017 IEEE international conference on robotics and automation (ICRA), pages 4568–4575. IEEE, 2017.
  * [21] Tung Dang, Marco Tranzatto, Shehryar Khattak, Frank Mascarich, Kostas Alexis, and Marco Hutter.  Graph-based subterranean exploration path planning using aerial and legged robots.  Journal of Field Robotics, 37(8):1363–1388, 2020.
  * [22] Chao Cao, Hongbiao Zhu, Fan Yang, Yukun Xia, Howie Choset, Jean Oh, and Ji Zhang.  Autonomous exploration development environment and the planning algorithms.  https://www.cmu-exploration.com/development-environment.
  * [23] Chao Cao, Hongbiao Zhu, Howie Choset, and Ji Zhang.  Exploring large and complex environments fast and efficiently.  In 2021 IEEE International Conference on Robotics and Automation (ICRA). IEEE.
  * [24] Howie M Choset, Kevin M Lynch, Seth Hutchinson, George Kantor, Wolfram Burgard, Lydia Kavraki, Sebastian Thrun, and Ronald C Arkin.  Principles of robot motion: theory, algorithms, and implementation.  MIT press, 2005.