text:
\section{Introduction\qquad \qquad}

Finite temperature effects are important in extremely high temperature
environments, such as in the early Universe a few seconds after the Big
Bang, astrophysical environments etc., where they are significant enough and
can not be ignored in comparison with the vacuum contribution. The high
temperature and density effects in ultra-relativistic plasma need to be
incorporated, for example, in Quantum Electrodynamics (QED) plasma, quark
gluon plasma and the core of dense stars. More recently renewed interest in
hot and dense QED plasmas has been generated due to the possibility of
creating ultra-relativistic electron positron plasmas with high-intensity
lasers ($\approx 10^{18}W/cm^{2}$) [1-3]. Two opposite laser pulses hitting
a thin gold foil can heat up the electrons in the foil up to several MeV
(\sim 10^{10}K$ $)$.

The particles propagating in vacuum can be assumed to be the ones with
interactions switched off. When these particles propagate through a medium,
several kinds of interaction processes take place. This makes the properties
of the system different from that in which all\ the particles are assumed to
be completely independent of each other, behaving as freely propagating bare
particles. When dealing with finite temperature environments in QED, where
the particles propagate in statistical background at energies around the
thresholds for particle antiparticle pair production, the temperature
effects need to be appropriately taken into account. These effects arise due
to continuous electron and photon exchanges between particles during the
physical interactions that take place in a heat bath containing hot
particles and antiparticles. The net statistical effects of the background
electrons and photons enter the theory through the fermion and boson
distributions respectively. Finite temperature calculations also provide a
guideline to estimate the density corrections, through chemical potential
effects of the background plasma, for higher order loop corrections.

The thermal background effects are included through the radiative
corrections [4,5]. Self energies and the wavefunctions of the propagating
particles acquire temperature corrections in this environment due to
exchanges of energy and momentum with real particles. The exact state of all
these background particles is unknown since they continually fluctuate
between different configurations. Temperatures of interest in such a
situation are in the range of a few MeV. Thermal propagators in real time
formalism comprise of temperature dependent terms added to the particle
propagators in vacuum theory [6]. In finite temperature electrodynamics,
electric fields are screened due to such interactions.

We use the real time formulation [7] for calculations of wavefunction
renormalization as a second order perturbative correction in $\alpha $ due
to the ease of obtaining the temperature corrections as additive terms to
the usual contribution in vacuum. Here we prefer that the loops with
temperature dependent momenta are integrated before temperature independent
variables in the relevant order $\alpha ^{2}$\ loops in QED, and therefore,
review the electron self energy in Ref. [8,9] earlier. This makes the
calculations of integrations over loop momenta much more simpler and easier
to handle. The results are obtained in a generalized form such that
intermediate temperatures $T\sim m$ are also included while the ranges of
high temperature $T>>m$ and low temperature $T<<m$, are retrieved\ from them
as the limiting cases.

In literature, the ways to compute finite temperature effects on
phase-space, vertex, mass corrections and photon emission or absorption
[10-19 ] are extensively discussed. The finite temperature wave function
renormalization has been dealt with several approaches [12--21],
specifically in the context of weak decay rates during primordial
nucleosynthesis. They agree on using finite temperature Dirac spinors to
obtain the corresponding effective projection operator. Differences in the
spinors presented in Refs. [20] and [21] were also pointed out [22].
However, their results for the case of $\beta $-decay and related processes
agreed with the ones which can be obtained using the approaches that had
already existed for wave function renormalization, except for the case of a
scalar boson decay in fermion-antifermion pair. We work out the two loop
corrections to the finite temperature wave function renormalization in the
generalized temperature framework, for the first time. Section 2 is based on
the re-examined and simplified calculations of loop correction upto two
orders in $\alpha $ that contribute to electron self energy in this
background. The expression for the relative correction in electron mass is
redone and the wave function renormalization constant is calculated in
section 3. Section 4 gives discussion of the results.

\section{Loop corrections to electron self energy\qquad \qquad \qquad \qquad}

At the one loop level, Feynman diagrams are calculated in the usual way by
substituting the finite temperature electron, positron and photon
propagators in place of those in vacuum. In real time formalism, the finite
temperature terms remain separate at order $\alpha $ since the terms
depending on temperature (hot) are additive to temperature independent
(cold) terms in the propagators. Therefore, at the one loop level the hot
and cold loop momenta are integrated separately.

Due to interactions with the background, the electron and positron masses
are known to get enhanced at one-loop and higher loop levels [4-9]. The
photons also acquire dynamically generated mass due to plasma screening
effect [23, 24, 27]. The presence of effective mass implies the fact that
the propagating particles constantly interact with the background.The
radiatively generated thermal mass creates a mass shift and in physical
quantities this acts as a kinematical cut-off, e.g., while determining the
production rate of particles in the heat bath.

Higher order loop corrections are required to study the perturbative
behavior at finite temperature. The two loop integrals comprise a
combination of cold and hot momenta which appear due to an overlap of
temperature dependent and temperature independent terms from the particle
propagators. The loop integrations involve an overlap of finite and
divergent terms due to which these become analytically more complicated. In
such situations, a prefered approach needs to be adopted for integrating
overlapping hot and cold loop momenta, i.e., to integrate over hot loop
momenta before integrating over cold ones, even at the two loop level [24].
This not only helps to simplify the loop integrations but allows one to
handle the statistical effects more appropriately.

The problem of renormalization in finite temperature field theories is
somewhat different from that at zero temperature due to the additional hot
infrared divergences at finite temperatures. The temperature, however, acts
as a regularization parameter for the hot ultraviolet divergences. The
infrared divergences introduced in this framework are also appropriately
removable in particle decay processes via bremstrahlung emission and
absorption effects [5, 25]. This was studied in detail at order $\alpha $
for all the possible ranges of temperature valid in QED including $T\sim m$
[25, 27]. The renormalization of QED\ was also established at the one loop
level for all the relevant ranges in temperatures and chemical potential
[23, 25--27]. \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad

\[
\FRAME{itbpF}{1.9692in}{0.9556in}{0in}{}{}{Figure}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "T";width 1.9692in;height 0.9556in;depth
0in;original-width 2.5832in;original-height 1.2393in;cropleft "0";croptop
"1";cropright "1";cropbottom "0";tempfilename
'LJSD6I00.wmf';tempfile-properties "XPR";}}
\
\[
\text{{\small Fig. 1. Two loop electron self energy diagrams and counter
terms }}
\
The electron self energy is once again calculated here at the two loop
level, from the point of view of renormalization, by integrating over loop
momenta in the electron self energy diagrams in Fig. 1(a) and 1(b). The
relevant counter terms required to cancel the divergences in these two loop
self energies are included in Fig. 1(c) and 1(d). Removal of overlapping
divergences has been even checked in QED up to the two loop level for $T\sim
m$, $T<<m$ and $T>>m\ $[9]. The integrations over the temperature dependent
momenta are re-examined here and wherever needed are re-done before
temperature independent momentum variables in loops for all the ranges of
temperature that are relevant in QED. The overlapping loops in Fig.1(a)
gives nonzero real contribution to b
\begin{eqnarray}
\Sigma ^{a}(p) &=&\frac{\alpha ^{2}}{2}{\LARGE [}\frac{1}{2\pi ^{3}}{\Large
\{}\frac{1}{\varepsilon }[2\NEG{I}-(\NEG{p}+6m)I_{A}+3(2m-\NEG{p})J_{A}+2\NE
{J}_{B}]-3(\NEG{p}+4m)I_{A}  \nonumber \\
&&+4\NEG{I}+(12m-5\NEG{p})J_{A}+2\NEG{J}_{B}{\Large \}}+(-1)^{r+1}e^{-r\beta
E}\dsum\limits_{n,r,s=1}^{\infty }{\Large \{[}\frac{3T^{2}}{4}{\LARGE \{
f_{+}(n,r)  \nonumber \\
&&\times {\Large [}mf_{+}(s,r){\Large (}\frac{\mathbf{\gamma }.\mathbf{p}}{
\mathbf{p}|^{2}}{\Large )}^{2}-f_{-}(s,r)h_{-}(p,\gamma )\frac{\mathbf
\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}}{\Large ]}-f_{-}(n,r){\Large [
f_{+}(s,r)h_{-}(p,\gamma )\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2
}  \nonumber \\
&&-\frac{1}{m}{\Large \{}1-\frac{2}{3}\gamma ^{0}\frac{\mathbf{\gamma }
\mathbf{p}}{v|\mathbf{p}|}+{\Large (}\frac{\mathbf{\gamma }.\mathbf{p}}{v
\mathbf{p}|}{\Large )}^{2}{\Large \}}f_{-}(s,r){\Large ]}{\LARGE \}}+2
{\Large [}(4+\NEG{p}\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}
)f_{+}(s,r)  \nonumber \\
&&+\NEG{p}h_{-}(p,\gamma ){\Large \{}f_{-}(s,r)-f_{-}(n,r)\frac{1}{m}\frac
I_{C}}{8\pi }{\Large \}}+\frac{I_{B}}{8\pi }{\Large \{(}4-\NEG{p}\frac
\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}}{\Large )}f_{+}(n,r)  \nonumber
\\
&&+\frac{{\Large (}\NEG{p}+m{\Large )}I_{C}}{8\pi }{\Large \}]}{\LARGE \}
+(-1)^{s}{\LARGE [}T^{2}{\Large \{}f_{+}(n,r){\Large [}\frac{\mathbf{\gamma
.\mathbf{p}}{|\mathbf{p}|^{2}}{\Large (}1-3m\frac{\mathbf{\gamma }.\mathbf{p
}{|\mathbf{p}|^{2}}{\Large )}f_{+}(s,r)  \nonumber \\
&&+3h_{-}(p,\gamma )\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}
f_{-}(s,r){\Large ]}-{\Large [\{}\frac{1}{m}h_{-}(p,\gamma )-3m{\Large (
\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}}{\Large )}^{2}{\Large \}
f_{+}(s,r)  \nonumber \\
&&-\frac{3}{m}h_{+}(p,\gamma )f_{-}(s,r){\Large ]}f_{-}(n,r){\Large \}}-
{\Large \{}{\LARGE [}{\Large \{(}5\NEG{p}+3m^{2}\frac{\mathbf{\gamma }
\mathbf{p}}{|\mathbf{p}|^{2}}{\Large )}\frac{\mathbf{\gamma }.\mathbf{p}}{
\mathbf{p}|^{2}}-12{\Large \}}f_{+}(n,r)  \nonumber \\
&&+\frac{5\NEG{p}}{m}h_{-}(p,\gamma )f_{-}(n,r){\LARGE ]}\func{Ei}_{-}
{\Large [}\frac{1}{m}f_{-}(n,r)h_{+}(p,\gamma )-\frac{\mathbf{\gamma }
\mathbf{p}}{|\mathbf{p}|^{2}}f_{+}(n,r)h_{-}(p,\gamma ){\Large ]}  \nonumber
\\
&&\times 3E\func{Ei}_{+}+{\Large \{[}3h_{-}(p,\gamma )\frac{\mathbf{\gamma }
\mathbf{p}}{v^{2}}f_{+}(n,r)-\frac{1}{m}f_{-}(n,r)\{3|\mathbf{p
|^{2}h_{+}(p,\gamma )  \nonumber \\
&&-5E\NEG{p}h_{-}(p,\gamma )\}{\Large ][}\frac{2e^{-rm\beta }}{m}\sinh
sm\beta +\frac{(r\func{Ei}_{+}+s\func{Ei}_{-})}{T}{\Large ]}-{\Large [}\frac
(r\func{Ei}_{+}-s\func{Ei}_{-})}{T}  \nonumber \\
&&+\frac{2e^{-rm\beta }}{m}\cosh sm\beta {\Large ]}m{\Large [\{
h_{-}(p,\gamma )-3m{\Large (}\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p
|^{2}}{\Large )}^{2}{\Large \}}f_{-}(n,r)-\frac{m}{2}\frac{\mathbf{\gamma }
\mathbf{p}}{|\mathbf{p}|^{2}}  \nonumber \\
&&\times {\Large (}1-3m\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}
{\Large )}f_{+}(n,r){\Large ]}-{\Large [}2\gamma ^{0}T{\LARGE \{}1+\frac{T}{
}f_{+}(s,r){\LARGE \}}+\NEG{p}T{\LARGE \{}\frac{2}{m}h_{-}(p,\gamma )
\nonumber \\
&&-\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}}{\Large \}]
f_{-}(s,r)+\{E\func{Ei}_{+}-m\gamma ^{0}\func{Ei}_{-}\}+\NEG{p}{\Large \{
\frac{m^{2}}{2}\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}}+\frac
2E^{2}}{m}h_{-}(p,\gamma ){\Large \}}  \nonumber \\
&&\times {\Large \{}\frac{e^{-m\beta (s+r)}-e^{-m\beta (r-s)}}{m}+\beta (
\func{Ei}_{+}-s\func{Ei}_{-}){\Large \}]}\frac{I_{C}}{8\pi }{\Large ]\}].}
\end{eqnarray
wher
\[
f_{\pm }(n,r)={\Large \{}\frac{1}{(n+r)}\pm \frac{1}{(n-r)}{\Large \}};\
f_{\pm }(s,r)={\Large \{}\frac{1}{(s+r)}\pm \frac{1}{(s-r)}{\Large \}},
\
\[
h_{\pm }(p,\gamma )=(\gamma ^{0}\pm \frac{\mathbf{\gamma }.\mathbf{p}}{v
\mathbf{p}|}),
\
\[
\func{Ei}_{\pm }=\func{Ei}[-m\beta (r+s)]\pm \func{Ei}[-m\beta (r-s),
\
\[
I_{A}=8\pi \int_{0}^{\infty }\frac{dk}{k}n_{B}(k),
\
\[
I_{B}=8\pi \dsum\limits_{r=1}^{\infty }(-1)^{r}\int_{0}^{\infty }\frac{dk}{k
e^{-r\beta (p-k)}n_{B}(k),
\
\[
I_{C}=8\pi \dsum\limits_{r=1}^{\infty }(-1)^{r}\int_{0}^{\infty }\frac{dk}{k
e^{-r\beta k}n_{B}(k),
\
\[
J^{A}\simeq -8\pi b(m\beta ),
\
\[
\frac{I^{0}}{E}=-\frac{2\pi ^{3}T^{2}}{3E^{2}v}\ln \frac{1-v}{1+v},
\
\[
\frac{\mathbf{I.p}}{\mathbf{p}^{2}}=-\frac{2\pi ^{3}T^{2}}{3E^{2}v^{3}}\{\ln
\frac{1-v}{1+v}+2v\},
\
\[
\frac{J_{B}^{0}}{E}\simeq 4\pi \lbrack \frac{T}{\mathbf{p}E}\ln \frac{1+v}
1-v}\{ma(m\beta )-Tc(m\beta )\}-3b(m\beta )],
\
\[
\frac{\mathbf{J_{B}.p}}{\mathbf{p}^{2}}\simeq \frac{\pi }{v^{2}E^{2}
[\{E^{2}-\frac{2}{3}m^{2}\}b(m\beta )+4T\{\frac{1}{v}\ln \frac{1+v}{1-v
+2\}\{ma(m\beta )-Tc(m\beta )\}],
\
with $v=\frac{|\mathbf{p}|}{p_{_{0}}},$ $(p_{_{0}}=E)$
\[
a(m\beta )=\ln (1+e^{-m\beta }),
\
\[
b(m\beta )=\dsum\limits_{_{n=1}}^{\infty }(-1)^{n}\func{Ei}(-nm\beta ),
\
\[
c(m\beta )=\dsum\limits_{_{n=1}}^{\infty }(-1)^{n}\frac{e^{-nm\beta }}{n^{2}
,
\
\ and $\func{Ei}(-x)$ is the error integral given b
\[
\func{Ei}(-x)=-\int_{x}^{\infty }\frac{dt}{t}e^{-t}.
\
The non vanishing real terms from loop within loop correction in Fig.1(b)
are also re-examined and recalculated, wherever needed, by retaining the\
specific order of integrating loop momenta, i.e., on the integration\ over
the variables in hot momenta before the cold momenta, we ge
\begin{eqnarray}
\Sigma _{\beta }^{b}(p) &=&\alpha ^{2}{\LARGE \{}\frac{2T^{2}}{3m^{2}}(\NE
{p}-m)+mT\dsum\limits_{n,r,s=1}^{\infty }(-1)^{s+r}{\LARGE [}{\Large \{
\beta (r+s)\func{Ei}[-m\beta (r+s)]  \nonumber \\
&&\text{ \ }+\frac{e^{-m\beta (s+r)}}{m}{\Large \}\{}\frac{1}{(n-r)}[\frac{2
}{m}+\gamma ^{0}(\frac{1}{2}-\frac{E^{2}}{m^{2}})+r\beta (m-\frac{\NEG{p}}{2
)]\text{\ \ \ \ \ \ \ \ \ \ \ \ \ }\   \nonumber \\
&&+\frac{1}{(n-s)}{\Large [}h(p,\gamma )+m\frac{\mathbf{\gamma }.\mathbf{p}}
|\mathbf{p}|^{2}}{\Large ]\}}+\frac{1}{2(n-r)}{\Large \{}\frac{2}{m}-\frac
E\gamma ^{0}}{m^{2}}  \nonumber \\
&&+r\beta \lbrack h(p,\gamma )+m\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{
}|^{2}}]\func{Ei}[-m\beta (r+s)]{\Large \}}{\LARGE ]}  \nonumber \\
&&-\frac{\pi T}{6|\mathbf{p}|}\dsum\limits_{n,r,s=1}^{\infty
}(-1)^{r+1}e^{-n\beta E}{\LARGE [}h(p,\gamma ){\Large \{[}1+(-1)^{s}]\frac
e^{-m\beta (r-n-s)}}{r-n-s}{\Large \}}  \nonumber \\
&&-T\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}}\frac{e^{-m\beta
(r-n)}}{r-n}+{\Large (}2-\frac{m\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2
}{\Large )\{}m\beta (r-n)\func{Ei}[-m\beta (r-n)]  \nonumber \\
&&-e^{-m\beta (r-n)}+\beta \lbrack 1+(-1)^{1+s}]\func{Ei}[-m\beta (r-n-s)
{\Large \}}{\LARGE ]\}.}
\end{eqnarray
In Eqs. (1) and (2), the preferred order of integration not only
sufficiently eases the calculations but the results are also simpler as
compared to those in Ref. [9] where this preference was not realized. The
electron self energies at the two loop level in Fig 1(a) and 1(b) are then
combined and rearranged to obtain the temperature corrections to the
electron mass and wavefunction.

\section{The Wavefunction Renormalization}

To incorporate finite temperature effecst on physical processes beyond the
tree level, one needs to have a consistent method of temperature dependent
renormalization. As already mentioned, renormalizability of the electron
mass was done through the order by order cancellation of singularities up to
two loop level. It can be easily checked that the second order\ in $\alpha $
correction is much smaller than the first order contribution so that the
perturbative behavior is valid. In a background, with $T\neq 0$, the Lorentz
invariance is broken and momentum independent renormalization constant is no
longer sufficient. Donoghue and Holstein used the temperature dependent
propagator to modify electron mass as well as the spinors accordingly [5].

The shift in the electron mass due to finite temperature effects is
calculated here from Eqs. (1) and (2). For this all the finite terms in
electron self energy upto second order in $\alpha $ are put together.
Following Ref. [5] the physical mass of the electron at one loop was
obtained in Ref. [25]\ in generalized form, by writin
\[
\Sigma (p)=A(p)E\gamma _{_{0}}-B(p)\vec{p}.\vec{\gamma}-C(p),
\
where $A(p)$, $B(p)$, and $C(p)$ are the relevant coefficients. Taking the
inverse of the propagator with momentum and mass term separated a
\[
S^{-1}(p)=(1-A)E\gamma ^{o}-(1-B)p.\gamma -(m-C),
\
and the physical mass $m_{phy}=m+\delta m^{(1)}+\delta m^{(2)},$ was deduced
by locating the pole of the propagator $\frac{i(\NEG{p}+m)}
p^{2}-m^{2}+i\varepsilon }$ . $\delta m^{(1)}$ and $\delta m^{(2)}$\ is the
shift in electron mass due to temperature effects at one and two loop level
respectively.$~$Using the same procedure, the relative shift in electron
mass at the two loop level was obtained [9]. This is recalculated here,
wherever required, and after recombining the similar summations it becomes
\begin{eqnarray}
\frac{\delta m^{(2)}}{m} &=&2\alpha ^{2}\dsum\limits_{r=1}^{\infty }{\LARGE
}\frac{T^{2}}{m^{2}}{\LARGE \{}\dsum\limits_{n=3}^{r+1}\text{ }(-1)^{n+r+1
\frac{\pi m\text{ }}{6|\mathbf{p}|}\frac{e^{-\beta (rE+mn)}}{n}\   \nonumber
\\
&&-\frac{3}{8}(-1)^{r}\frac{e^{-r\beta E}}{|\mathbf{p}|^{2}}{\Large [}\frac
9E^{2}}{2m^{2}}+6\dsum\limits_{s=3}^{r+1}\frac{1}{s}+4\dsu
\limits_{n,s=3}^{r+1}\frac{1}{ns}+(-1)^{s-r}{\Large \{}\frac{9E}{m
(3+4\dsum\limits_{s=3}^{r+1}\frac{1}{s}{\Large )}  \nonumber \\
&&+2{\Large (}\frac{|\mathbf{p}|^{2}}{m^{2}}-3{\Large )(}9+18\dsu
\limits_{s=3}^{r+1}\frac{1}{s}+8\dsum\limits_{n,s=3}^{r+1}\frac{1}{ns
{\Large )\}]}+\frac{4}{E^{2}v^{2}}{\LARGE \}}-\frac{m^{2}}{\pi ^{2}}c(m\beta
)  \nonumber \\
&&-\frac{T}{m}{\LARGE \{}\frac{\pi \text{ }}{6|\mathbf{p}|
\dsum\limits_{s=2}^{r+1}\dsum\limits_{n=1}^{s+1}\frac{e^{-\beta (rE+mn)}}{n
{\LARGE [}1{\LARGE \ }-{\Large \{}(-1)^{r+n}-(-1)^{s+n}{\Large \}}{\LARGE ]}
\nonumber \\
&&+{\Large [\{}\func{Ei}(-m\beta )-\func{Ei}(-2m\beta ){\Large \}\{}\frac{9
}{4}{\Large (}\frac{E}{|\mathbf{p}|^{2}}-\frac{1}{m}{\Large )}+{\Large (
\frac{5E}{m}-21+\frac{E^{2}}{2m^{2}}{\Large )}\dsum\limits_{n=3}^{r+1}\frac{
}{n}{\Large \}}  \nonumber \\
&&+{\Large \{}\frac{9}{4v^{2}}-\dsum\limits_{n=1}^{s+1}\dsu
\limits_{s=3}^{r+1}{\Large [}1-E^{2}{\Large (}\frac{1}{2m^{2}}+\frac{3}{
\mathbf{p}|^{2}}{\Large )}+\frac{3E}{m}{\Large ]\}}(-1)^{s}\func{Ei
(-sm\beta )]  \nonumber \\
&&+e^{-rm\beta }{\LARGE \{[}\frac{9E}{2v^{2}}+2{\LARGE (}\frac{3E}{v^{2}}
\frac{3|\mathbf{p}|^{2}}{m}-5E{\LARGE )}\dsum\limits_{n=3}^{r+1}\frac{1}{n
{\LARGE ]}\dsum\limits_{s=1}^{\infty }\sinh sm\beta  \nonumber \\
&&-\frac{3m^{3}}{|\mathbf{p}|^{2}}{\Large (}\frac{3}{4}-\dsu
\limits_{n=3}^{r+1}\frac{1}{n}{\LARGE )}\dsum\limits_{s=1}^{\infty }\cosh
sm\beta {\LARGE \}]\}}+{\Large {\LARGE \{}}\frac{9m}{4|\mathbf{p}|^{2}
{\Large (}E^{3}+\frac{m^{3}}{2}{\Large )}  \nonumber \\
&&+{\Large [}\frac{3m}{|\mathbf{p}|^{2}}(E^{3}+m^{3})+5mE-3|\mathbf{p}|^{2
{\Large ]}\dsum\limits_{n=3}^{r+1}\frac{1}{n}{\LARGE \}}{\Large \{}\func{Ei
(-m\beta )-2\func{Ei}(-2m\beta ){\Large \}}  \nonumber \\
&&-\dsum\limits_{n=3}^{r+1}{\LARGE \{}\dsum\limits_{s=1}^{r+1}\frac{(-1)^{s
}{n}{\Large [}\frac{m^{2}r}{2}e^{-sm\beta }+{\LARGE \{}s{\LARGE (}2mE-\frac
E^{3}}{m}{\LARGE )}+\frac{m^{2}(s-r)}{2}{\LARGE \}}\func{Ei}(-sm\beta
{\Large ]}  \nonumber \\
&&-\frac{\pi m^{2}\text{ }}{3|\mathbf{p}|}{\LARGE [}e^{-\beta rE}\text{
(-1)^{n+r}(n+1){\LARGE \ }-\dsum\limits_{s=2}^{r+1}(-1)^{n+s}{\LARGE ]}\func
Ei}(-nm\beta ){\LARGE \}].}
\end{eqnarray}

From the reviewed expression for the electron self energy obtained in Eqs.
(1) and (2) the relation for the wave function renormalization constant is
derived. This comes out to b
\begin{eqnarray}
Z_{2}^{-1} &=&\frac{\partial \Sigma }{\partial \NEG{p}}  \nonumber \\
&=&1-\alpha {\LARGE [}\frac{1}{4\pi }\left( \frac{3}{\varepsilon }-4\right)
\frac{5}{\pi }b(m\beta )+\frac{I_{A}}{4\pi ^{2}}  \nonumber \\
&&-\frac{T^{2}}{\pi vE^{2}}\ln \frac{1+v}{1-v}{\LARGE \{}\frac{\pi ^{2}}{6
-c(m\beta )+m\beta a(m\beta ){\LARGE \}]}  \nonumber \\
&&-\alpha ^{2}{\LARGE [}\frac{1}{4\pi ^{3}}\{\frac{1}{\varepsilon
(I_{A}-3J_{A})+(3I_{A}+5J_{A})\}-\frac{2T^{2}}{3\pi ^{3}m^{2}}  \nonumber \\
&&+\frac{m}{8\pi }\dsum\limits_{n,r,s=1}^{\infty }(-1)^{s+r}{\Large \{}r
\frac{e^{-m\beta (s+r)}}{m}-\beta (r+s)\func{Ei}\{-m\beta (r+s)\}]{\Large \}}
\nonumber \\
&&+\frac{1}{8}\dsum\limits_{n,r,s=1}^{\infty }(-1)^{r}T{\LARGE \{}e^{-r\beta
E}{\LARGE [}f_{+}(s,r)\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}}
\frac{I_{B}I_{C}}{64\pi ^{2}}  \nonumber \\
&&+h(p,\gamma )\{f_{-}(n,r)\frac{I_{C}}{8\pi }-f_{-}(s,r)\}+f_{+}(n,r)\frac
\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}}\frac{I_{B}}{8\pi }{\LARGE ]}
\nonumber \\
&&+{\LARGE [\{}5\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}
f_{+}(n,r)-\frac{5}{m}h(p,\gamma )f_{-}(n,r){\LARGE \}}\func{Ei}_{-}
\nonumber \\
&&+\frac{5E}{m^{2}}h(p,\gamma )f_{-}(n,r){\LARGE \{}\frac{2e^{-rm\beta }}{m
\sinh sm\beta +\beta (r\func{Ei}_{+}+s\func{Ei}_{-}){\LARGE \}}  \nonumber \\
&&+{\LARGE \{}\frac{2}{m}h(p,\gamma )-\frac{\mathbf{\gamma }.\mathbf{p}}{
\mathbf{p}|^{2}}{\LARGE \}}f_{-}(s,r){\LARGE ]\}}+{\LARGE \{}\frac{m^{2}}{2
\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}}+\frac{2E^{2}}{m
h(p,\gamma ){\LARGE \}}  \nonumber \\
&&\times {\LARGE \{}\frac{2e^{-rm\beta }}{m}\sinh sm\beta +\beta (r\func{Ei
_{+}-s\func{Ei}_{-}){\LARGE \}}\frac{I_{C}}{8\pi }{\LARGE ].}
\end{eqnarray
From this expression for $Z_{2}^{-1},$ not only the behavior at intermediate
temperatures $T\sim m$\ can be extracted but the ranges of high temperature
T>>m$, low temperature $T<<m$, can be also retrieved\ from it as limiting
cases.

\section{Results and Discussion}

With the preferred order of integration of hot loops before the cold ones,
the previously calculated self-mass correction terms are redone, wherever
required, for all the possible ranges of temperature. This preference
simplifies the calculations since the statistical effects are taken care of
through hot loop momenta integrations here, before the zero temperature
integration variables are dealt with. Therefore, we have re-written the
electron self energy expressions in QED at the two loop level that were
presented in Ref. [9]. \ This led to the modfied expression for the relative
change in the electron mass at the two loop level in Eq. (3). From these
corrections one can then retrieve the results for all temperature ranges of
interest here, classified as, the high temperature $T>>m$ $($having $m\beta
\longrightarrow 0$ with $e^{-m\beta }$ falling off exponentially as compared
to $\frac{T^{2}}{m^{2}}$), the low temperature $T<<m$ (with fermions
contribution negligible) and the intermediate temperatures $T\sim m$ (by
taking $m\beta \rightarrow 1$).

We calculated here, for the first time in this generalized form, the
wavefunction renormalization constant using the thermal contributions to the
second order self energy diagrams. The divergences get cancelled as usual by
including the counter terms in Fig. 1(c) and 1(d), as already checked [9].
The fermions do not pick any contribution from the heat bath at low
temperature. Therefore, the second order in $\alpha $ corrections to the
electron self energy at low temperature can be retrieved as a limiting case
that contains contribution from hot photons only giving
\[
\Sigma _{\beta }(p)\text{ }_{\longrightarrow }^{^{T<<m}}\frac{\alpha ^{2}}
4\pi ^{3}}\left[ 4\NEG{I}+\frac{8\pi T^{2}}{3m^{2}}(\NEG{p}-m)\right] .
\
The wave function renormalization constant upto two loops at $T<<m$\ from
Eq. (4) is therefor
\begin{eqnarray}
&&Z_{2}^{-1}\bigskip \text{ $_{\longrightarrow }^{^{T<<m}}$ }1+\frac{\alpha
}{4\pi }\left( 4-\frac{3}{\varepsilon }\right) -\frac{\alpha }{4\pi ^{2}
\left( I_{A}-\frac{I^{0}}{E}\right)  \nonumber \\
&&\text{ \ \ \ \ \ \ \ \ \ \ \ }-\frac{\alpha ^{2}}{4\pi ^{2}}\left( 3+\frac
1}{\varepsilon }\right) I_{A}+\frac{2\alpha ^{2}T^{2}}{3\pi ^{2}m^{2}},
\end{eqnarray
which is the same as\ that in Ref. [28]. The high temperature limit for this
constant gives
\begin{eqnarray}
&&Z_{2}^{-1}\text{ }_{\longrightarrow }^{^{T>>m}}\text{ }1-\alpha {\Large [
\frac{2I_{A}}{\pi }+\frac{1}{4\pi }\left( \frac{3}{\varepsilon }-4\right)
\frac{4\pi T^{2}}{3}{\Large ]}  \nonumber \\
&&\text{ \ \ \ \ \ \ \ \ \ \ \ }-\alpha ^{2}{\LARGE [}\frac{1}{4\pi ^{3}
{\LARGE \{}\frac{1}{\varepsilon }(I_{A}-3J_{A})+(3I_{A}+5J_{A})-\frac{8T^{2
}{3m^{2}}{\LARGE \}}  \nonumber \\
&&\text{ \ \ \ \ \ \ \ \ \ \ \ }+\frac{1}{8}\dsum\limits_{n,r,s=1}^{\infty
}(-1)^{r}T{\LARGE \{}e^{-r\beta E}{\LARGE [}f_{+}(s,r)\frac{\mathbf{\gamma }
\mathbf{p}}{|\mathbf{p}|^{2}}-\frac{I_{B}I_{C}}{64\pi ^{2}}  \nonumber \\
&&\text{ \ \ \ \ \ \ \ \ \ \ \ }+h(p,\gamma )\{f_{-}(n,r)\frac{I_{C}}{8\pi
-f_{-}(s,r)\}+f_{+}(n,r)\frac{\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}
\frac{I_{B}}{8\pi }  \nonumber \\
&&\text{ \ \ \ \ \ \ \ \ \ \ \ }+{\LARGE \{}\frac{2}{m}h(p,\gamma )-\frac
\mathbf{\gamma }.\mathbf{p}}{|\mathbf{p}|^{2}}{\LARGE \}}f_{-}(s,r){\LARGE
]\}].}
\end{eqnarray}

It can be seen from Eq.(7) that the leading contribution in this range of
temperature, $T>>m$, at the two loop level is $\frac{2T^{2}}{3m^{2}}$. The
two loop fermion self energy in QED has been calculated in detail recently
[29] using the hard thermal loop resummation introduced by Braaten and
Pisarski [30]. As far as the renormalization is concerned, the resummed hard
thermal loops (HTL) do not affect it [29]. Hence $Z_{2}$ does not get any
contribution from HTL here. Moreover, it is worth noticing that the thermal
corrections, at second order in $\alpha ,$ to the wavefunction
renormalization constant at extreme temperatures ( $T<<m$ and $T>>m$ ) are
still proportional to $\frac{T^{2}}{m^{2}}$ as in case of the selfmass of
electron. However, the expression for the intermediate temperatures is
significantly different from the self-energy expression obtained earlier
[9]. The calculations around $T\sim m$ for the mass and wavefunction
renormalizations are still cumbersome at the two-loop level, but much less
difficult than those in Ref. [9]. The renormalizability of the theory at
finite temperature can be explicitly checked and holds through the order by
order cancellation of singularities. This provides a platform to include the
general effects due to chemical potential in the hot and dense background,
later on. With the experience of including chemical potential at one loop
level, in real time formulation [23,26,27], it is foreseen that two loop
self energies will be much more complicated but is still worth-doing to
develop a calculational technique for high density hot plasmas or even
superfluids inside the cores of neutron stars[31]. The modified wavefunction
is expected to modify the finite temperature contributions to electroweak
processes [32] as well as the neutrino magnetic moments [33] up to the two
loop level.

\begin{quotation}
\textbf{Acknowledgement}
\end{quotation}

{\small One of the authors (MQH) thanks Higher Education Commission,
Pakistan for providing partial funding under a research grant during this
work.}

\begin{quotation}
\textbf{Figure caption}
\end{quotation}

{\small Fig. 1 \ Two loop electron self energy diagrams and counter terms}

\begin{quotation}
\textbf{References}
\end{quotation}

\begin{enumerate}
\item {\small E.P. Liang, S.C. Wilks, and M. Tabak, Pair production by
ultraintense lasers, Phys. Rev. Lett. 81 (1998) 4887--4890.}

\item {\small B. Shen, J. Meyer-ter-Vehn, Pair and }$\gamma ${\small -photon
production from a thin foil confined by two laser pulses, Phys. Rev. E 65
(2001) 016405--016410.}

\item {\small M.H. Thoma, Ultrarelativistic electron-positron plasma, Eur.
Phys. J. D 55 (2009) 271--278; M.H. Thoma, Colloquium: Field theoretic
description of ultrarelativistic electron positron plasmas, Rev. Mod. Phys.
81 (2009) 959--968.}

\item {\small H.A. Weldon, Covariant calculations at finite temperature: The
relativistic plasma, Phys. Rev. D 26 (1982) 1394--1407.}

\item {\small J.F. Donoghue and B.R. Holstein, Renormalization and radiative
corrections at finite-temperature, Phys. Rev. D 28 (1983) 340--348; ibid 29
(1983) 3004(E).}

\item {\small P. Landsman and Ch G. Weert, Real and imaginary time field
theory at finite temperature and density, Phys. Rep. 145 (1987) 141--249.}

\item {\small V.V. Klimov, Spectrum of elementary Fermi excitations in quark
gluon plasma, Sov. J. Nucl. Phys. 33 (1981) 934--935; L. Dolan and R.
Jackiw, Symmetry behavior at finite temperature, Phys. Rev. D 9 (1974)
3320--3341; S. Weinberg, Gauge and global symmetries at high temperature,
Phys. Rev. D 9 (1974) 3357--3378; C. Bernard, Feynman rules for gauge
theories at finite temperature, Phys. Rev. D 9 (1974) 3312-- 3320. }$\ \ \ $

\item {\small Mahnaz Qader (Haseeb), Samina S. Masood, and K. Ahmed,
Second-order electron mass dispersion relation at finite temperature, Phys.
Rev. D 44 (1991) 3322--3327.}

\item {\small Mahnaz Qader (Haseeb), Samina S. Masood, and K. Ahmed,
Second-order electron mass dispersion relation at finite temperature. II.,
Phys. Rev. D 46 (1992) 5633--5647.}

\item {\small D.A. Dicus, E.W. Kolb, A.M. Gleeson, E.C.G. Sudarshan,
Primordial nucleosynthesis including radiative, Coulomb, and
finite-temperature corrections to weak rates, V.L. Teplitz, and M.S. Turner,
Phys. Rev. D 26 (1982) 2694--2709.}

\item {\small J.-L. Cambier, J.R. Primack, and M. Sher, Finite temperature
radiative corrections to neutron decay and related processes, Nucl. Phys.
B209 (1982) 372--388.}

\item {\small J.F. Donoghue, B.R. Holstein, and R.W. Robinett, Quantum
electrodynamics at finite temperature, Ann. Phys. (N.Y.) 164 (1985) 233--276
}

\item {\small A.E. Johansson, G. Peresutti, and B.S. Skagerstam, Quantum
field theory at finite temperature: Renormalization and radiative
corrections, Nucl. Phys. B278 (1986) 324--342.}

\item {\small W. Keil, Radiative corrections renormalization and at
finite-temperature: A real time approach, Phys. Rev. D 40 (1989) 1176--1198.}

\item {\small R. Baier, E. Pilon, B. Pire, and D. Schiff, Finite-temperature
radiative corrections to early universe neutron-proton ratio: Cancellation
of infrared and mass singularities, Nucl. Phys. B336 (1990) 157--183.}

\item {\small W. Keil and R. Kobes, Mass and wave function renormalization
at finite temperature, Physica A 158 (1989) 47--57.}

\item {\small M. LeBellac and D. Poizat, Renormalization of external lines
in relativistic field theories at finite temperature, Z. Phys. C 47 (1990)
125--132.}

\item {\small T. Altherr and P. Aurenche, Fermion self-energy corrections in
perturbative theory at finite temperature, Phys. Rev. D 40 (1989) 4171--4177
}

\item {\small R.L. Kobes and G.W. Semeneff, Discontinuity of Green functions
in field theory at finite temperature and density (I), Nucl. Phys. B260
(1985) 714--746; Discontinuity of Green functions in field theory at finite
temperature and density (II), B272 (1986) 329--364.}

\item {\small R.F. Sawyer, Temperature-dependent wave function
renormalization and weak interaction rates prior to nucleosynthesis, Phys.
Rev. D 53 (1996) 4232--4236.}

\item {\small I.A. Chapman, Finite temperature wave-function
renormalization, a comparative analysis, Phys. Rev. D 55 (1997) 6287--6291.}

\item {\small S. Esposito, G. Mangano, G. Miele, and O. Pisanti, Wave
function renormalization at finite temperature, Phys. Rev. D 58 (1998)
105023-1--5.}

\item {\small K. Ahmed and Samina S. Masood, Vacuum polarization at finite
temperature and density in QED, Ann. Phys. (N.Y.) 207 (1991) 460--473.}

\item {\small Samina S. Masood and Mahnaz Q. Haseeb, Second-order
corrections to QED coupling at low temperature, Int. J. of Mod. Phys. A 23
(2008) 4709--4719.}

\item {\small K. Ahmed and Samina Saleem (Masood), Renormalization and
radiative corrections at finite-temperature reexamined, Phys. Rev. D 35
(1987) 1861--1871.}

\item {\small K. Ahmed and Samina Saleem (Masood), Finite-temperature and
-density renormalization effects in QED, Phys. Rev. D 35 (1987) 4020--4023.}

\item {\small Samina S. Masood, Photon mass in the classical limit of
finite-temperature and -density QED, Phys. Rev. D 44 (1991) 3943--3948;
Samina S. Masood, Renormalization of QED in superdense media, Phys. Rev. D\
47 (1993) 648--652.}

\item {\small Mahnaz Q. Haseeb and Samina S. Masood, Two loop low
temperature corrections to electron self energy, Chin. Phys. C, published
online Nov. 29, 2010 (in press).}

\item {\small E. Mottola and Z. Sz\'{e}p, Systematics of high temperature
perturbation theory: The two-loop electron self energy in QED, Phys. Rev. D
81 (2010) 025014--025052.}

\item {\small E. Braaten and R.\thinspace D. Pisarski, Soft amplitudes in
hot gauge theories: A general analysis, Nucl. Phys. B 337 (1990) 569--634;
E. Braaten and R.\thinspace D. Pisarski, Calculation of the gluon damping
rate in hot QCD, Phys. Rev. D 42 (1990) 2156--2160; R.\thinspace D.
Pisarski, Scattering amplitudes in hot gauge theories, Phys. Rev. Lett. 63
(1989) 1129--1132.}

\item {\small D. Page, et. al., Rapid cooling of the neutron star in
cassiopeia a triggered by neutron superfluidity in dense matter, Phys. Rev.
Lett. 106 (2011) 081101-081108. }

\item {\small Mahnaz Qader (Haseeb) and Samina S. Masood, Finite-temperature
and -density corrections to electroweak processes, Phys. Rev. D 46 (1992)
5110--5116. }

\item {\small Samina S. Masood, Neutrino physics in hot and dense media,
Phys. Rev. D 48 (1993) 3250--3258.}
\end{enumerate}

\end{document}


text:
\section{System Overview}
\label{sec:overview}

The proposed framework operates upon a voxel grid map.
As illustrated in Fig.\ref{fig:overview}, it is composed of an incremental update of the FIS (Sect.\ref{sec:frontier}) and a hierarchical exploration planning approach (Sect.\ref{sec:planning}).
Whenever the map is updated using sensor measurements, it is examined whether any frontier clusters are influenced.
If so, FISs of influenced clusters are removed while new frontiers along with their FISs are extracted (Sect.\ref{sec:frontier}).
After that, the exploration planning is triggered, which finds global exploration tour, refines local viewpoints, and generates trajectory to a selected viewpoint successively (Sect.\ref{sec:planning}).
The exploration is considered finished if no frontier exists.

\end{document}
\section{Introduction}
\label{sec:intro}

\IEEEPARstart{U}{nmanned} aerial vehicles, especially quadrotors have gained widespread popularity in various applications, such as inspection, precision agriculture, and search and rescue.
Among the tasks, autonomous exploration, in which the vehicle explores and maps the unknown environment to gather information, is a fundamental component.

Various exploration planning methods have been proposed in recent years, with some real-world experiments presented \cite{cieslewski2017rapid,schmid2020efficient,meng2017two,selin2019efficient}.
However, most of them demonstrate a low/medium exploration rate, which is unsatisfactory for many large-scale real-world applications.
First of all, many existing planners plan exploring motions in greedy manners, such as maximizing the immediate information gain, or navigating to the closest unknown region.
The greedy strategies ignore global optimality and therefore result in low overall efficiency.
Besides, most methods generate rather conservative motions in order to guarantee information gain and safety simultaneously in previously unknown environments.
Low-speed exploration, however, disallows quadrotors to fully exploit their dynamic capability to fulfill the missions.
Last but not least, many methods suffer from high computational overhead and can not respond quickly and frequently to environmental changes.
However, to enable faster exploration, it is desirable to replan new motions immediately whenever new information of the environment is available.

\begin{figure}[t]
   \begin{center}
      \subfigure[A cluttered environment for the exploration tests.]
      {\includegraphics[width=0.75\columnwidth]{large1} \label{fig:intro1}}
      \subfigure[The online built map and executed trajectory.]
      {\includegraphics[width=0.75\columnwidth]{large1_rviz} \label{fig:intro2}}
      \vspace{-0.4cm}
	\end{center}
   \caption{\label{fig:intro} A quadrotor autonomous exploration test conducted in a complex indoor scene. Video of the experiments is available at: \url{https://www.youtube.com/watch?v=_dGgZUrWk-8}.
   }
   \vspace{-0.7cm}
\end{figure}

Motivated by the above facts, this paper proposes \textbf{FUEL}, a hierarchical framework that can support \textbf{F}ast \textbf{U}AV \textbf{E}xp\textbf{L}oration in complex environments.
We introduce a frontier information structure (FIS), which contains essential information in the entire space required by exploration planning.
The structure can be updated efficiently and incrementally when new information is collected, so it is capable of supporting high-frequency planning.
Based on the FIS, a hierarchical planner generates exploring motion in three coarse-to-fine steps.
It starts by finding a global exploration tour that is optimal in the context of accumulated environment information.
Then, viewpoints on a local segment of the tour are refined, further improving the exploration rate.
Finally, a safe, dynamically feasible and minimum-time trajectory is generated.
The planner produces not only efficient global coverage paths, but also safe and agile local maneuvers.
Moreover, the planner is triggered whenever unvisited regions are explored, so that the quadrotor always responds promptly to any environmental changes, leading to consistently fast exploration.

We compare our method with classic and state-of-the-art methods in simulation.
Results show that in all cases our method achieves complete exploration in much shorter time (3-8 times faster).
What's more, we conduct fully onboard real-world exploration in various challenging environments.
Both the simulation and real-world tests demonstrate an unprecedented performance of our method compared to state-of-the-art ones.
To benefit the community, we will make the source code public.
The contributions of this paper are summarized as follows:

1) An incrementally updated FIS, which captures essential information of the entire explored space and facilitates exploration planning in high frequency.

2) A hierarchical planning method, which generates efficient global coverage paths, and safe and agile local maneuvers for high-speed exploration.

3) Extensive simulation and real-word tests that validate the proposed method. The source code of our system will be made public.

\end{document}
\section{Related Work}
\label{sec:related}

\subsection{Exploration Path Planning}
\label{subs:related_exploration}

Robotic exploration, which uses mobile robots to map unknown environments, has been studied for years.
Some of the works focus on exploring the space quickly\cite{cieslewski2017rapid, dharmadhikari2020motion}, as this paper does.
Meanwhile, other methods place more emphasis on accurate reconstruction\cite{schmid2020efficient, song2017online}.
Among the various proposed methods, the frontier-based approaches are one type of classic approaches.
The methods are first introduced in \cite{yamauchi1997frontier} and evaluated more comprehensively afterwards in \cite{julia2012comparison}.
To detect frontiers in 3D space, a stochastic differential equation-based method is proposed in \cite{shen2012stochastic}.
In the original method\cite{yamauchi1997frontier}, the closest frontier is selected as the next target.
\cite{cieslewski2017rapid} presented a different scheme.
In every decision, it selects the frontier within the FOV that minimizes the velocity change to maintain a consistently high flight speed.
This scheme is shown to outperform the classic method\cite{yamauchi1997frontier}.
In \cite{deng2020robotic}, a differentiable measure of information gain based on frontiers is introduced, allowing paths to be optimized with gradient information.

Sampling-based exploration, as another type of major approaches, generate viewpoints randomly to explore the space.
These methods are closely related to the concept of next best view (NBV) \cite{connolly1985determination}, which computes covering views repeatedly to obtain a complete model of a scene.
\cite{bircher2016receding} first used NBV in 3D exploration, in which it expands RRTs with accessible space and executes the edge with the highest information gain in a receding horizon fashion.
The method was extended to consider uncertainty of localization\cite{papachristos2017uncertainty}, visual importance of different objects \cite{dang2018visual} and inspection tasks \cite{bircher2018receding} later.
To avoid discarding the expanded trees directly, roadmaps are constructed in \cite{witting2018history, wang2019efficient} to reuse previous knowledge.
\cite{schmid2020efficient} maintains and refines a single tree continuously using a rewiring scheme inspired by RRT*.
To achieve faster flight, \cite{dharmadhikari2020motion} samples safe and dynamically feasible motion primitives directly and execute the most informative one.

There are also methods combining the advantages of frontier-based and sampling-based approaches.
\cite{charrow2015information,selin2019efficient} plan global paths toward frontiers and sample paths locally.
\cite{charrow2015information} also presented a gradient-based method to optimize the local path.
\cite{meng2017two} samples viewpoints around frontiers and finds the global shortest tour passing through them.
\cite{song2017online} generates inspection paths that cover the frontier completely using a sampling-based algorithm.

Most of existing methods make decision greedily and does not account for the dynamics of the quadrotor, which leads to inefficient global tours and conservative maneuvers.
In contrast, we plan tours that efficiently cover the entire environment and generate dynamically feasible minimum-time trajectories to enable agile flight.

\subsection{Quadrotor Trajectory Planning}
\label{subs:related_trajectory}

Trajectory planning for quadrotor has been widely studied, which can be categorized into the hard-constrained and soft-constrained approaches in major.
The former is pioneered by minimum-snap trajectory\cite{MelKum1105}, whose closed-form solution was presented \cite{RicBryRoy1312} later.
Based on \cite{MelKum1105}, \cite{CheLiuShe2016, fei2018icra, ding2019efficient} extract convex safe regions for safe trajectory generation.
To obtain a more reasonable time allocation, fast marching\cite{fei2018icra}, kinodynamic search\cite{ding2019efficient} and mixed integer-based methods \cite{tordesillas2019faster} are proposed.
\cite{fei2018icra} also introduced an efficient B\'ezier curve-based method to guarantee feasibility.

Soft-constrained methods typically formulate a non-linear optimization trading off several objectives.
Recently \cite{oleynikova2016continuous, fei2017iros, usenko2017real, zhou2019robust, zhou2020robust, zhou2020raptor} applied them to local replanning, demonstrating their attractive performance.
The methods were revived by \cite{ratliff2009chomp} and extended to continuous-time trajectories\cite{oleynikova2016continuous} later.
To relieve the issue of local minima, \cite{fei2017iros} initializes the optimization with collision-free paths.
\cite{usenko2017real} introduces uniform B-splines for replanning.
More recently, \cite{zhou2019robust} further exploited B-splines and demonstrated fast flight in field tests.
\cite{zhou2019robust} is further improved with topological guiding paths and perception-awareness in \cite{zhou2020robust, zhou2020raptor}.

In this paper, we base our trajectory planning on \cite{zhou2019robust} but extend it to optimize all parameters of B-splines.
In this way, the total trajectory time can be minimized so that the unknown space is explored with a higher navigation speed.

\end{document}
\section{Incremental Frontier Information Structure}
\label{sec:frontier}

As presented in classic frontier-based exploration \cite{yamauchi1997frontier}, frontiers are defined as known-free voxels right adjacent to unknown voxels, which are grouped into clusters to guide the navigation.
Traditionally, the extracted information is too coarse to do fine-grained decision making.
Besides, frontiers are retrieved by processing the entire map, which is not scalable for large scenes and high planning frequencies.
In this work, we extract richer information from frontiers to enable more elaborate planning, and develop an incremental approach to detect frontiers within the locally updated map.

\begin{figure}[t]
	\begin{center}
    {\includegraphics[width=0.99\columnwidth]{frontier}}
      \vspace{-0.8cm}
  \end{center}
   \caption{\label{fig:frontier_detection} Incremental frontier detection and clustering.
   Top: detecting and removing outdated frontiers.
   Bottom: new frontier is detected (left) and PCA is performed, the large cluster is split into two smaller ones (right).
   }
   \vspace{-0.8cm}
\end{figure}

\subsection{Frontier Information Structure}
\label{subs:frontier_structure}

A frontier information structure $ FI_i $ is computed when a new frontier cluster $ F_i $ is created.
It stores all cells $C_{i}$ belonging to the cluster and the average position $ \mathbf{p}_{\text{avg},i} $ of $ C_{i} $.
The axis-aligned bounding box (AABB) $ B_{i} $ of the cluster is also computed, in order to accelerate the detection of frontier changes (Sect.\ref{subs:frontier_detection}).
To serve the exploration planning (Sect.\ref{sec:planning}), candidate viewpoints $\textit{VP}_i$ are generated around the cluster.
Besides, a doubly linked list $ L_{\text{cost},i} $ containing connection costs between $ F_i $ and all other clusters is computed.
Data stored by a FIS is listed in Tab.\ref{tab:structure}.

\subsection{Incremental Frontier Detection and Clustering}
\label{subs:frontier_detection}

As depicted in Fig.\ref{fig:frontier_detection}, every time the map is updated by sensor measurements, the AABB of the updated region $ B_m $ is also recorded, within which outdated frontier clusters are removed and new ones are searched.
It starts with going through all clusters and returning only those whose AABBs ($B_i$) intersect with $ B_m $.
Then, precise checks are conducted for the returned clusters, among which the ones containing cells that are no longer frontier are removed.
These two processes are inspired by the broad/narrow phase collision detection algorithms\cite{ericson2004real}, which eliminate most unaffected clusters in a fast way and significantly reduce the number of expensive precise checks.

After the removal, new frontiers are searched and clustered into groups by the region growing algorithm, similar to the classic frontier-based method.
Among the groups, the ones with a small number of cells typically resulting from noisy sensor observations are ignored.
The remaining groups, however, may contain large-size clusters which are not conducive to distinguishing distinctive unknown regions and making elaborate decisions.
Therefore, we perform Principal Component Analysis (PCA) for each cluster and split it into two uniform ones along the first principal axis, if the largest eigenvalue exceeds a threshold.
The split is conducted recursively so that all large clusters are divided into small ones.

\begin{figure}[t]
	\begin{center}
		{\includegraphics[width=0.9\columnwidth]{viewpoint}}
      \vspace{-0.5cm}
	\end{center}
   \caption{\label{fig:viewpoint} Generating candidate viewpoints for a frontier cluster.
   Within the cylindrical coordinate system centered at the average position of the cluster, points are sampled uniformly.
   }

\end{figure}

\begin{table}[t]
  \centering
  \caption{\label{tab:structure} Data contained by a FIS $ FI_i $ of cluster $F_i$.}
  \begin{tabular}{cc}
  \hline\hline
   \textbf{Data}    &  \textbf{Explanation}    \\
  \hline

    $C_i$    &  Frontier cells that belong to the cluster    \\
    $ \mathbf{p}_{\text{avg},i} $    &  Average position of $C_i$  \\
    $ B_{i} $    &  Axis-aligned bounding box of $C_i$  \\
    $\textit{VP}_{i}$    &  Viewpoints covering the cluster  \\
    $ L_{\text{cost},i} $   &  Doubly linked list of connection costs to all other clusters  \\
  \hline\hline
  \end{tabular}
  \vspace{-1.0cm}
\end{table}

\subsection{Viewpoint Generation and Cost Update}
\label{subs:frontier_viewpoint}

Intuitively, a frontier cluster implies a potential destination to explore the space.
However, unlike previous approaches which simply navigate to the center of a cluster, we desire more elaborate decision making.
To this end, when a cluster $F_i$ is created, we generate a rich set of viewpoints $\textit{VP}_i = \left\{ \mathbf{x}_{i,1}, \mathbf{x}_{i,2}, \cdots, \mathbf{x}_{i,n_i} \right\} $ covering it, where $ \mathbf{x}_{i,j} = \left( \mathbf{p}_{i,j}, \xi_{i,j} \right) $.
$\textit{VP}_i$ are found by uniformly sampling points in the cylindrical coordinate system whose origin locates at the cluster's center, as displayed in Fig.\ref{fig:viewpoint}.
For each of the sampled points $ \mathbf{p} $ lying within the free space, the yaw angle $ \xi $ is determined as the one maximizing sensor coverage to the cluster, by using a yaw optimization method similar to\cite{witting2018history}.
The coverage is evaluated as the number of frontier cells that comply with the sensor model and are not occluded by occupied voxels.
Then, viewpoints with coverage higher than a threshold are reserved and sorted in descending order of coverage.
We reserve at most $ N_{\text{view}} $ viewpoints in $\textit{VP}_i$ ($n_i \le N_{\text{view}} $) to make the local viewpoint refinement (Sect.\ref{sec:planning}) tractable.

To perform global planning of exploration tour (Sect.\ref{sec:planning}), a connection cost between each pair of clusters $ (F_{k_1}, F_{k_2}) $ is required.
Let $ t_{\text{lb}}(\mathbf{x}_{k_1,j_1}, \mathbf{x}_{k_2,j_2}) $ denotes a time lower bound when moving between two viewpoints $ \mathbf{x}_{k_1,j_1} $ and $ \mathbf{x}_{k_2,j_2} $, it is computed by:
\begin{align}
  \label{equ:time_bound}
  \hspace{-0.2cm} t_{\text{lb}}&(\mathbf{x}_{k_1,j_1}, \mathbf{x}_{k_2,j_1}) = \max \bigg\{ \frac{\text{length}\left(P \left(  \mathbf{p}_{k_1,j_1},  \mathbf{p}_{k_2,j_2} \right)\right)}{v_{\text{max}}},  \\ \nonumber
  & \frac{\min\left( |\xi_{k_1,j_1}-\xi_{k_2,j_2}|, 2\pi-|\xi_{k_1,j_1}-\xi_{k_2,j_2}| \right)}{\dot{\xi}_{\text{max}}} \bigg\}
  \end{align}

where $ P \left(  \mathbf{p}_{k_1,j_1},  \mathbf{p}_{k_2,j_2} \right) $ denote a collision-free path between $ \mathbf{p}_{k_1,j_1} $ and $ \mathbf{p}_{k_2,j_2} $ found by a path searching algorithm, $ v_{\text{max}} $ and $ \dot{\xi}_{\text{max}} $ are the limits of velocity and angular rate of yaw.
For each pair $ (F_{k_1}, F_{k_2}) $, we select the viewpoints with highest coverage and estimate the cost as $ t_{\text{lb}}(\mathbf{x}_{k_1,1}, \mathbf{x}_{k_2,1}) $, in which $ P \left(  \mathbf{p}_{k_1,1},  \mathbf{p}_{k_2,1} \right) $ is searched on the voxel grid map using the A* algorithm.

Note that computing connection costs between all pairs of $ N_{\text{cls}} $ clusters from scratch requires $ O(N_{\text{cls}}^2) $ A* searching, which is considerably expensive.
Fortunately, the costs can also be computed in an incremental manner.
When any outdated clusters are removed (Sect.\ref{subs:frontier_detection}), associated cost items in $ L_{\text{cost},i} $ of all remaining FISs are erased.
After that, connection costs from each new cluster to all other clusters are computed and inserted into $ L_{\text{cost},i} $.
Suppose there are $ k_{\text{new}} $ new clusters in each frontier detection, the above update scheme takes $ O(k_{\text{new}} \cdot N_{\text{cls}} ) $ time.
Practically, $ k_{\text{new}} $ is small and can be regarded as a constant factor, resulting in a linear time update of connection costs.

\end{document}
\section{Results}
\label{sec:results}

\subsection{Implementation Details}
\label{subs:results_detail}
We set $ w_{\text{c}} = 1.5 $ in Equ.\ref{equ:start_connect_cost} and \ref{equ:cost_local}.
In global tour planning, the ATSP is solved using a Lin-Kernighan-Helsgaun heuristic solver \cite{helsgaun2000effective}.
In local viewpoint refinement, we reserve at most $ N_{\text{view}} = 15 $ viewpoints in each FIS and truncate the global tour within radius $ R_{\text{rf}}= 5.0 $.
For trajectory optimization, we use $ w_{\text{s,p}}=5.0 $, $ w_{\text{s},\xi}=2.5 $, $ w_{\text{t}}=1.0 $, $ \lambda_{\text{c}} = \lambda_{\text{bs}} = 10.0 $, $\lambda_{\text{d}} = 2.0$ and $ d_{\min} = 0.4 $ and solve the problem with a general non-linear optimization solver NLopt\footnote{https://nlopt.readthedocs.io/en/latest/}.
Cubic B-spline ($ p_b=3 $) is used as the trajectory representation.

To achieve fast exploration, an efficient mapping framework is essential.
In our work, we utilize a volumetric map \cite{han2019fiesta}, which has been successfully applied to fast autonomous flights\cite{zhou2019robust,zhou2020raptor} in complex scenes.
Similar to \cite{wurm2010octomap}, which is widely applied in exploration, \cite{han2019fiesta} builds an occupancy grid representation of the space.
Meanwhile it also maintains an ESDF incrementally to facilitate the trajectory planning.
For brevity we refer interested readers to \cite{han2019fiesta} for more details about our mapping framework.

In all field experiments, we localize the quadrotor by a visual-inertial state estimator\cite{qin2018vins}.
We use a geometric controller\cite{lee2010} for tracking control of the $ (x,y,z,\xi) $ trajectory.
We equipped our customized quadrotor platform with an Intel RealSense Depth Camera D435i.
All the above modules run on an Intel Core i7-8550U CPU.

\begin{figure}[t]
	\begin{center}
		{\includegraphics[width=0.8\columnwidth]{prog_bridge}}
		{\includegraphics[width=0.8\columnwidth]{prog_maze}}
      \vspace{-0.5cm}
	\end{center}
   \caption{\label{fig:progress} The exploration progress of four methods in the bridge (top) and large maze (bottom) scenarios.
   }
   \vspace{-0.3cm}
\end{figure}

\subsection{Benchmark and Analysis}
\label{subs:results_benchmark}

\begin{table}
   \centering
   \caption{\label{tab:computation_time} Average computation time of each proposed component.}
   \begin{tabular}{>{\centering}p{0.7cm}>{\centering}p{0.9cm}>{\centering}p{1.1cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}c}
   \hline\hline
   \textbf{Scene}  & \multicolumn{6}{c}{\textbf{Average computation time (ms)}}                                                                                \\
   \cline{2-7}
                   & \multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Frontier\\Sect.\ref{subs:frontier_detection}\end{tabular}} & \multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}View.+Cost\\Sect.\ref{subs:frontier_viewpoint}\end{tabular}} & \multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Global\\Sect.\ref{subs:planning_global}\end{tabular}} & \multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Local\\Sect.\ref{subs:planning_local}\end{tabular}} & \multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Traj.\\Sect.\ref{subs:planning_traj}\end{tabular}} & \multirow{2}{*}{Total}  \\
                   &                                                                          &                                                                           &                                                                      &                                                                      &                                                                      &                         \\
   Bridge          & 4.69                                                                     & 4.86+5.16                                                                 & 1.12                                                                 & 4.10                                                                 & 4.23                                                                 & 24.17                   \\
   Maze            & 5.21                                                                     & 6.06+10.97                                                                & 3.53                                                                 & 4.98                                                                 & 5.47                                                                 & 36.23                   \\
   \hline\hline
   \end{tabular}
   \vspace{-1.8cm}
\end{table}

We test our proposed framework in simulation.
We benchmark it in a bridge scenario and a large maze scenario.
Three methods are compared: the NBVP \cite{bircher2016receding}, the classic frontier method\cite{yamauchi1997frontier} and the rapid frontier method\cite{cieslewski2017rapid}.
Note that no open source code is available for \cite{cieslewski2017rapid}, so we use our implementation.
In all tests the dynamic limits are set as $ v_{\text{max}}=2.0$ m/s and $ \dot{\xi_{\text{max}}} = 0.9 $ rad/s for all methods.
The FOVs of the sensors are set as $ [80\times60] \ \text{deg} $ with a maximum range of $4.5 $ m.
In both scenarios each method is run for 3 times with the same initial configuration.
Statistics and exploration progresses of the four methods are shown in Tab.\ref{tab:benchmark} and \ref{fig:progress} respectively.
The computation time of each component of our method is listed in Tab.\ref{tab:computation_time}.

\subsubsection{Bridge Scenario}
\label{subsubs:bridge}

Firstly, we compare the four methods in a $ 10 \times 20 \times 5 \ \text{m}^3 $ space containing a bridge, as shown in Fig.\ref{fig:bridge_compare}.
The result indicates that we achieve much shorter exploration time and smaller time variance.
The overall exploration path of our method is significantly shorter, primarily because we plan tours globally.
The executed path is smoother, since we refine motions locally and generate smooth trajectories.
Also, we are able to navigate at a higher flight speed, owing to the minimum-time trajectory planning.

\begin{figure}[t]
	\begin{center}
		{\includegraphics[width=0.75\columnwidth]{maze1}}
		{\includegraphics[width=0.77\columnwidth]{maze2}}
      \vspace{-0.5cm}
	\end{center}
   \caption{\label{fig:four_compare_maze} Path generated by the proposed method (blue), classic frontier method\cite{yamauchi1997frontier}(red), rapid frontier method\cite{cieslewski2017rapid} (green) and NBVP\cite{bircher2016receding}(pink).
   }
   \vspace{-0.25cm}
\end{figure}

\subsubsection{Large Maze Scenario}
\label{subsubs:maze}

We also compare the methods quantitatively in a large maze environment shown in Fig.\ref{fig:four_compare_maze}.
The explored space is $ 20 \times 80 \times 3 \ \text{m}^3 $ large.
In this scenario, all the benchmarked methods take a long time to reach full coverage, due to the complexity of the scene.
In contrast, our method completes the exploration 4+ times faster on average.
Path executed by the four methods after completion are displayed in Fig.\ref{fig:four_compare_maze}.
Noticeably, our method explores the maze in a more sensible order, without revisiting the same place frequently.
In consequence, it produces a much shorter coverage path and an approximately linear exploration rate (Fig.\ref{fig:progress}).
This behavior is owing to the global plan, without which known regions may be revisited many times and slow down the progress, as the benchmarked methods do.
Also note that the computation time in the large maze is longer, mostly due to the larger scene, which naturally involves a greater number of frontier clusters.

\begin{figure}[t]
	\begin{center}
		{\includegraphics[width=0.493\columnwidth]{huge_part2}}
		{\includegraphics[width=0.493\columnwidth]{huge_part1}}
		{\includegraphics[width=0.75\columnwidth]{huge1_rviz}}
    \vspace{-0.5cm}
	\end{center}
  \caption{\label{fig:indoor} Experiments in an indoor scene composed of two room: a small one with tables and chairs (top left), a large room cluttered with obstacles (top right).
  }
  \vspace{-1.4cm}
\end{figure}

\begin{figure}[t]
	\begin{center}
		{\includegraphics[width=0.75\columnwidth]{forest1}}
		{\includegraphics[width=0.75\columnwidth]{forest_rviz1}}
      \vspace{-0.3cm}
	\end{center}
  \caption{\label{fig:forest} Exploration experiment conducted in a forest.
  }
  \vspace{-0.4cm}
\end{figure}

\subsection{Field Exploration Tests}
\label{subs:results_field}

To further validate the proposed method, we conduct extensive field experiments in both indoor and outdoor environments.
In all tests we set the dynamics limits as $ v_{\text{max}}=1.5$ m/s, $ a_{\text{max}}=0.8 $ m/s and $ \dot{\xi}_{\text{max}} = 0.9 $ rad/s.
Note that we do not use any external device for localization and only rely on the onboard state estimator.

First, we present fast exploration tests in two indoor scenes.
The first scene is shown in Fig.\ref{fig:intro}, within which we deploy dozens of obstacles and the quadrotor should perform 3D maneuvers to map the unknown space and avoid obstacles simultaneously.
We bound the space to be explored with a $ 10 \times 6 \times 2  \ \text{m}^3 $ box.
One sample map and the flight trajectory is presented in Fig.\ref{fig:intro}.
The second indoor scene is a larger environment including two rooms, where one room is similar to scene 1 and the other is a part of an office containing tables and chairs.
The space is bounded by a $ 15 \times 11 \times 2 \ \text{m}^3 $ box.
The quadrotor starts by exploring the large room, after which it proceeds to the small one.
The second scene, the online generated map and trajectory are shown in Fig.\ref{fig:indoor}.
Note that in the two scenes the quadrotor starts out at a spot with low visibility, so it only maps a small region of the environment at the beginning.
Finally, to validate our method in natural environments, we conduct exploration tests in a forest.
The size of the area to explore is $ 11 \times 10 \times 2 \ \text{m}^3 $.
The experiment environment and the associated results are displayed Fig.\ref{fig:forest}.

The above experiments demonstrate the capability of our method in complex real-world scenarios.
They also show the merits of our autonomous quadrotor system, among which the state estimation\cite{qin2018vins} and mapping modules\cite{han2019fiesta} are also crucial to fullfil the real-world tasks.
Video demonstration of all experiments is available (Fig.\ref{fig:intro}), we refer the readers to it for more details.

\end{document}
\section{Hierarchical Exploration Planning}
\label{sec:planning}

Instead of adopting greedy exploration strategies or generating conservative maneuvers, we produce global paths to cover the frontiers efficiently and plan safe and agile motions for faster flight.
Our planner takes inspiration from the recent hierarchical quadrotor planning paradigm \cite{CheLiuShe2016, fei2018icra, zhou2019robust}, and makes decisions in three phases, as shown in Fig.\ref{fig:planning}.

\begin{figure}[t]
	\begin{center}
		{\includegraphics[width=0.85\columnwidth]{planning}}
      \vspace{-0.5cm}
	\end{center}
   \caption{\label{fig:planning} Generating exploration motion in three coarse-to-fine steps:
   (1) the shortest tour covering frontier clusters in the entire environment is found.
   (2) a local segment of the global tour is refined.
   (3) a safe minimum-time trajectory is generated to the first viewpoint on the refined tour.
   }
   \vspace{-0.9cm}
\end{figure}

\subsection{Global Exploration Tour Planning}
\label{subs:planning_global}

Our exploration planning begins with finding a global tour to cover existent frontier clusters efficiently.
Inspired by\cite{meng2017two}, we formulate it as a variant of the Traveling Salesman Problem (TSP), which computes an open-loop tour starting from the current viewpoint and passing viewpoints at all clusters.
We reduce this variant to a standard Asymmetric TSP (ATSP) that can be solved quickly by available algorithms, by properly designing the engaged cost matrix $ \mathbf{M}_{\text{tsp}} $.

Assume there are $ N_{\text{cls}} $ clusters totally, $ \mathbf{M}_{\text{tsp}} $ corresponds to a $ N_{\text{cls}}+1$ dimensions square matrix.
The major part is the $ N_{\text{cls}} \times N_{\text{cls}} $ block composed of the connection cost between each pair of frontier clusters, which is computed as:
\begin{align}
  & \mathbf{M}_{\text{tsp}}(k_1,k_2) = \mathbf{M}_{\text{tsp}}(k_2,k_1) \\ \nonumber
  &= t_{\text{lb}}(\mathbf{x}_{k_1,1}, \mathbf{x}_{k_2,1}), \ k_1, k_2 \in \left\{1,2,\cdots,N_{\text{cls}} \right\}
\end{align}
As mentioned in Sect.\ref{subs:frontier_viewpoint}, this information is maintained when frontiers are detected.
Thus, the $ N_{\text{cls}} \times N_{\text{cls}} $ block can be filled without extra overhead.

The first row and column of $\mathbf{M}_{\text{tsp}}$ are associated with the current viewpoint $\mathbf{x}_0 = \left( \mathbf{p}_0, \xi_0 \right)$ and $N_{\text{cls}}$ clusters.
Starting from $\mathbf{x}_0$, the cost to the $k$-th cluster is evaluated by:
\begin{align}
  \label{equ:start_connect_cost}
  \hspace{-0.5cm} \mathbf{M}_{\text{tsp}}(0,k) = t_{\text{lb}}(\mathbf{x}_0,\mathbf{x}_{k,1}) + w_{\text{c}} \cdot c_{\text{c}}(\mathbf{x}_{k,1}),\\ \nonumber
   k \in \left\{1,2, \cdots, N_{\text{cls}} \right\}
\end{align}
here a motion consistency cost $ c_{\text{c}}(\mathbf{x}_{k,1}) $ is introduced, which is generally computed as:
\begin{equation}
   c_{\text{c}}(\mathbf{x}_{k,j}) = \cos^{-1} \frac{(\mathbf{p}_{k,j}-\mathbf{p}_0)\cdot \mathbf{v}_0}{\left\| \mathbf{p}_{k,j}-\mathbf{p}_0 \right\| \left\|\mathbf{v}_0 \right\|}
\end{equation}
where $\mathbf{v}_0$ is the current velocity.
In some cases, multiple tours have comparable time lower bound, so back-and-forth maneuvers may be generated in successive planning steps and slow down the progress.
We eliminate this inconsistency with $ c_{\text{c}}(\mathbf{x}_{k,1}) $, which penalizes large changes in flight direction.

Our problem is different from standard TSP whose solution is a closed-loop tour.
However, we can reduce it to an ATSP by assigning zero connection costs from other clusters to the current viewpoint:
\begin{equation}
   \mathbf{M}_{\text{tsp}}(k,0) = 0, \ k \in \left\{0, 1,2, \cdots, N_{\text{cls}} \right\}
\end{equation}
In this way, going back to the current viewpoint in any closed-loop tours contributes no extra cost, so each closed-loop tour always contains an open-loop one with an identical cost.
As a result, we can obtain the optimal open-loop tour by finding the optimal closed-loop one and retrieving its equal-cost open-loop tour.

\begin{figure}[t]
	\begin{center}
		{\includegraphics[width=0.8\columnwidth]{refine}}
      \vspace{-0.5cm}
	\end{center}
   \caption{\label{fig:planning_local} Refining viewpoints locally using the graph search approach.
   Along a truncated segment of the global tour, multiple viewpoints of each visited cluster are considered to select the optimal set of viewpoints.
   }
   \vspace{-0.5cm}
\end{figure}

\subsection{Local Viewpoint Refinement}
\label{subs:planning_local}

The global tour planning finds a promising order to visit all clusters.
Nonetheless, it involves only a single viewpoint of each cluster, which are not necessarily the best combination among all viewpoints.

To this end, a richer set of viewpoints on a truncated segment of the global tour are considered to further improve the exploration rate, by using a graph search approach, as depicted in Fig.\ref{fig:planning_local}.
We found consecutive clusters whose viewpoints on the global tour are closer than $ R_{\text{rf}} $ to the current position.
To simplify the notation, suppose that $ F_i, 1 \le i \le N_{\text{rf}} $ are the considered clusters.
We create graph nodes for their viewpoints $ \textit{VP}_i $ and the current viewpoint $ \mathbf{x}_0 $.
Then each node is connected to other nodes associated with the next cluster with a directed edge, which compose a directed acyclic graph capturing possible variation of the truncated tour.
We utilize the Dijkstra algorithm to search for the optimal local tour $ \Xi = \left\{ \mathbf{x}_{1,j_1}, \mathbf{x}_{2,j_2}, \cdots, \mathbf{x}_{N_{\text{rf}},j_{N_{\text{rf}}}} \right\} $ that minimizes the cost:
\begin{align}
\label{equ:cost_local}
   c_{\text{rf}}(\Xi) &= t_{\text{lb}}(\mathbf{x}_0,\mathbf{x}_{1,j_1}) + w_{\text{c}} \cdot c_{\text{c}}( \mathbf{x}_{1,j_1} ) \\ \nonumber
   &+ t_{\text{lb}}(\mathbf{x}_{N_{\text{rf}},j_{N_{\text{rf}}}},\mathbf{x}_{N_{\text{rf}}+1,1}) + \sum_{k=1}^{N_{\text{rf}}-1} t_{\text{lb}}(\mathbf{x}_{k,j_k},\mathbf{x}_{k+1,j_{k+1}})
\end{align}
which also consists of time lower bounds and motion consistency.
Note that it is straight forward to incorporate information gain\cite{bircher2016receding,selin2019efficient} to Equ.\ref{equ:cost_local}, however, evaluating information gain for numerous viewpoints is expensive.
Practically, we find that simply adopting viewpoints based on their coverages is much faster and leads to consistently satisfactory results.

\subsection{Minimum-time B-spline Trajectory}
\label{subs:planning_traj}

Given the discrete viewpoints, continuous trajectories are required for smooth navigation.
Our quadrotor trajectory planning is based on a method\cite{zhou2019robust} that generates smooth, safe and dynamically feasible B-spline trajectories.
We go one step further to optimize all parameters of B-splines, so that the total trajectory time is minimized to enable the quadrotor to fully utilize its dynamic capability.

As the quadrotor dynamics are differentially flat\cite{MelKum1105}, we plan trajectories for the flat outputs $ \mathbf{x} \in \left( x,y,z,\xi \right) $.
Let $ \mathbf{X}_{\text{cb}} = \left\{ \mathbf{x}_{\text{c},0}, \mathbf{x}_{\text{c},1}, \cdots, \mathbf{x}_{\text{c},N_b} \right\} $ where $ \mathbf{x}_{\text{c},i} = \left( \mathbf{p}_{\text{c},i}, \xi_{\text{c},i} \right) $ be the $ N_b+1 $ control points of a $ p_b $ degree uniform B-spline, and $ \Delta t_{b} $ be the knot span.
We find the B-spline that trades-off smoothness and total trajectory time, and satisfies safety, dynamic feasibility and boundary state constraints.
It can be formulated as an following optimization problem:
\begin{equation}
   \underset{ \mathbf{X}_{\text{c},b}, \Delta t_b }{\arg \min } \ f_{\text{s}}+ w_{\text{t}} T + \lambda_{\text{c}} f_{c}+ \lambda_{\text{d}}\left( f_{\text{v}} + f_{\text{a}} \right) + \lambda_{\text{bs}} f_{\text{bs}}
\end{equation}
Similar to \cite{zhou2019robust}, $ f_{\text{s}} $ is the elastic band smoothness cost:
\begin{equation}
   f_{\text{s}} = \sum_{i = 0}^{N_b-2} \mathbf{s}_{i}^{\text{T}} \mathbf{R}_{\text{s}} \mathbf{s}_{i}, \ \mathbf{s}_{i} = \mathbf{x}_{\text{c},i+2}-2\mathbf{x}_{\text{c},i+1}+\mathbf{x}_{\text{c},i}
\end{equation}
in which $ \mathbf{R}_{\text{s}} $ is the penalty matrix:
\begin{equation}
   \mathbf{R}_{\text{s}} = \begin{bmatrix}
      w_{\text{s},\text{p}}\mathbf{I}_3 & \mathbf{0} \\
      \mathbf{0}^{\text{T}} & w_{\text{s},\xi}
      \end{bmatrix}
\end{equation}
$ T $ is the total trajectory time depending on $ \Delta t_b $ and the number of B-spline segments:
\begin{equation}
   T = (N_{b}+1-p_b) \cdot \Delta t_b
\end{equation}
$ f_{\text{c}} $, $ f_{\text{v}} $ and $ f_{\text{a}} $ are the penalties to ensure safety and dynamic feasibility.
Given the following function:
\begin{equation}
	\mathcal{P}(\tau_1, \tau_2) = \left\{
	\begin{array}{cl}
	(\tau_1-\tau_2)^{2} & \ \  \tau_1 \le \tau_2 \\
	0 & \ \ \text{else}
	\end{array}
	 \right.
\end{equation}
$ f_{\text{c}} $ is evaluated as:
\begin{equation}
\label{equ:collision_cost}
   f_{\text{c}} = \sum\limits_{i=0}^{N_{b}} \ \mathcal{P}\left( d\left(\mathbf{p}_{\text{c},i}\right), d_{\text{min}} \right)
\end{equation}
where $ d(\mathbf{p}_{\text{c},i}) $ is the distance of point $\mathbf{p}_{\text{c},i}$ to the nearest obstacle, which can be obtained from the Euclidean signed distance field (ESDF) maintained by the mapping module.
Practically a clearance larger than our quadrotor's radius (typically $d_{\min} \ge$ 0.5 m) ensures safety in complex scenes.
$ f_{\text{v}} $ and $ f_{\text{a}} $ penalize infeasible velocity and acceleration:
\begin{equation}
\label{equ:vel_cost}
f_{\text{v}} = \sum_{i=0}^{N_b-1} \left\{\sum\limits_{\mu \in \{x,y,z\}} \mathcal{P}(v_{\text{max}}, |\dot{p}_{\text{c},i,\mu}|) + \mathcal{P}(\dot{\xi}_{\text{max}}, |\dot{\xi}_{\text{c},i}|) \right\}
\end{equation}
\begin{equation}
\label{equ:acc_cost}
f_{\text{a}} = \sum_{i=0}^{N_b-2} \left\{\sum\limits_{\mu \in \{x,y,z\}} \mathcal{P}(a_{\text{max}}, |\ddot{p}_{\text{c},i,\mu}|) + \mathcal{P}(\ddot{\xi}_{\text{max}}, |\ddot{\xi}_{\text{c},i}|) \right\}
\end{equation}
in which the control points of derivatives are utilized:
\begin{align}
    & \dot{\mathbf{x}}_{\text{c},i} = \left[ \dot{p}_{\text{c},i,x}, \dot{p}_{\text{c},i,y}, \dot{p}_{\text{c},i,z}, \dot{\xi}_{\text{c},i} \right]^\text{T} =  \frac{\mathbf{x}_{\text{c},i+1} - \mathbf{x}_{\text{c},i}}{\Delta t_b} \\
    & \ddot{\mathbf{x}}_{\text{c},i} = \left[ \ddot{p}_{\text{c},i,x}, \ddot{p}_{\text{c},i,y}, \ddot{p}_{\text{c},i,z}, \ddot{\xi}_{\text{c},i} \right]^\text{T} =  \frac{\mathbf{x}_{\text{c},i+2} - 2\mathbf{x}_{\text{c},i+1} + \mathbf{x}_{\text{c}, i}}{\Delta t_{b}^{2}}
\end{align}
In Equ.\ref{equ:collision_cost}, \ref{equ:vel_cost} and \ref{equ:acc_cost}, the convex hull property of B-spline is utilized to ensure the feasibility efficiently.
For brevity we refer the reader to \cite{zhou2019robust} for more details.

Lastly, we set the 0th to 2nd order derivatives at the start to the instantaneous state $ ( \mathbf{x}_0, \dot{\mathbf{x}}_0, \ddot{\mathbf{x}}_0) $ for smooth motion.
The 0-th order derivative at the end is also determined by the viewpoint $ \mathbf{x}_{\text{next}} $ to be visited.
In implementation we use cubic B-splines, so the associated cost is:
\begin{align}
\label{equ:cost_boundary}
& f_{\text{bs}} = \left\| \frac{\mathbf{x}_{\text{c},0}+4\mathbf{x}_{\text{c},1}+\mathbf{x}_{\text{c},2}}{6} - \mathbf{x}_{0} \right \|^2 + \left\| \frac{ \dot{\mathbf{x}}_{\text{c},0}+ \dot{\mathbf{x}}_{\text{c},1}}{2} - \dot{\mathbf{x}}_0 \right\|^2 \\ \nonumber
   &+ \left\| \ddot{\mathbf{x}}_{\text{c},0} - \ddot{\mathbf{x}}_0 \right\|^2 + \left\| \frac{\mathbf{x}_{\text{c},N_b-2}+4\mathbf{x}_{\text{c},N_b-1}+\mathbf{x}_{\text{c},N_b}}{6} - \mathbf{x}_{\text{next}} \right \|^2
\end{align}

\end{document}
\section{Conclusions}
\label{sec:conclude}
In this paper, we propose a hierarchical framework for rapid autonomous quadrotor exploration.
An incrementally maintained FIS is introduced to provide the exploration planning with essential information.
Based on FIS, a hierarchical planner plans exploration motions in three sequential steps, which finds efficient global tours, selects a local set of optimal viewpoints, and generates minimum-time local trajectories.
The method makes decisions at high frequency to respond quickly to environmental changes.
Both benchmark and real-world tests show the competence of our method.

One limitation of our method is assuming perfect state estimation, as most methods do.
We evaluate our method in simulation with ground truth localization, while drifts in pose are not considered.
However, error in state estimation exists generally and can not be ignored.
In the future we plan to consider the state estimation uncertainty in our method and evaluate its performance under pose drifts.

\addtolength{\textheight}{0.cm}

\newlength{\bibitemsep}\setlength{\bibitemsep}{0.0\baselineskip}
\newlength{\bibparskip}\setlength{\bibparskip}{0.0pt}
\let\oldthebibliography\thebibliography
\renewcommand\thebibliography[1]
\oldthebibliography{#1
\setlength{\parskip}{\bibitemsep
\setlength{\itemsep}{\bibparskip
}

redpajama_set_name: RedPajamaArXiv