# Femtosecond optical breakdown in silicon

Tzveta Apostolova  Institute for Nuclear Research and Nuclear Energy,
Bulgarian Academy of Sciences, Tsarigradsko chausse 72, 1784 Sofia, Bulgaria
Institute for Advanced Physical Studies, New Bulgarian University, 1618 Sofia,
Bulgaria  Boyan Obreshkov  Institute for Nuclear Research and Nuclear Energy,
Bulgarian Academy of Sciences, Tsarigradsko chausse 72, 1784 Sofia, Bulgaria

###### Abstract

We investigate photoinization, energy deposition, plasma formation and the
ultrafast optical breakdown in crystalline silicon irradiated by intense near-
infrared laser pulses with pulse duration $\tau\leq$ 100 fs. The occurrence of
high-intensity breakdown was established by the sudden increase of the
absorbed laser energy inside the bulk, which corresponds to threshold energy
fluence $\Phi_{th}>$ 1 J/cm2. The optical breakdown is accompanied by severe
spectral broadening of the transmitted pulse. For the studied irradiation
conditions, we find that the threshold fluence increases linearly with the
increase of the pulse duration, while the corresponding laser intensity
threshold decreases. The effect of the high plasma density on the stability of
diamond lattice is also examined. For near threshold fluences, when about 5 %
of valence electrons are promoted into the conduction band, the Si-Si bonds
are softened and large Fermi degeneracy pressure arises (with pressure up to
100 kbar). The mechanical instability of the diamond lattice suggests that the
large number of electron-hole pairs leads directly to ultrafast melting of the
crystal structure.

##  I Introduction

Ultrafast breakdown resulting in damage of dielectrics has been subject of
extensive experimental and theoretical investigations LIB1997 ; Du1994 ;
Stuart1996 ; Kautek1996 ; Varel1996 ; Lenzner1998 ; Pronko1998 ; Li1999 ;
vanDriel1987 ; Apostolova2000 ; Schaffer2001 ; Apostolova2002 ; Gamaly2002 ;
Kudryashov2007 ; Rajeev2009 ; Chimier2011 ; Uteza2011 ; Grojo2013 ;
Lebugle2014 ; Sato2015 ; Zhokhov2018 . Laser-induced breakdown in dielectrics
has been modeled in terms of three basic physical processes, cf. e.g.
Stuart1996 ; Schaffer2001  : i) multiphoton and tunnel ionization generating
free electrons into the conduction band, ii) avalanche ionization involving
free carrier absorption followed by impact ionization and iii) transfer of the
electronic excitation energy to the lattice, eventually resulting in optical
damage after conclusion of the pulse. The relative importance of these
mechanisms have been studied by measuring the dependence of the breakdown
threshold on the duration $\tau$ of the applied laser pulse Stuart1996 ;
Du1994 ; Lenzner1998 ; Chimier2011 ; Uteza2011 . Though experimental results
for the measured damage threshold agree for picosecond to nanosecond pulses, a
controversy exists on the observed scaling trend of the threshold fluence for
femtosecond pulses Du1994 ; Stuart1996 ; Lenzner1998 . Though measurement of
the energy fluence threshold has been extended to sub 10-fs laser pulses
Lenzner1998 ; Chimier2011 , it is difficult to extract information regarding
the time scale for transfer of the electronic excitation energy to the
lattice. Femtosecond time-resolved pump-probe measurements of optical
breakdown thresholds in dielectrics Li1999  have demonstrated that the energy
transfer to the lattice is an ultrafast process that lasts for 100fs.

Femtosecond time-resolved measurements are very effective in probing non-
linear photionization and the temporal behavior of the optical breakdown and
damage thresholds. Non-thermal melting of silicon is tought to be a result of
promotion of a large fraction of the valence electrons into the conduction
band Shank1983 ; Hulin1984 ; Tom1988 ; Sokol1995 ; Rousse2001 ; Sundaram2002 ;
Harb2008 . The high kinetic pressure of the hot electron distribution softens
the transverse acoustic phonon at the Brillouin zone boundary, which results
in instability of the diamond lattice against shear stress Stampfli1990 ;
Stampfli1992 ; Stampfli1999 ; Silvestrelli1996 . Theoretical modeling of the
ultrafast electron dynamics in dielectrics and insulators has been reviewed in
Refs.Gamaly2011 ; Balling2013 . More recently, density-functional-theory based
modeling of ultrafast electron dynamics in silicon were reported Sato2014p1 ;
Sato2014p2 ; Schultze2014 ; Apostolova2020 ; Otobe2020 ; Yamada2021 . The
fundamental physical processes involve coupling of the laser energy into the
electron subsystem on an attosecond time scale as a part of the subcycle
photoionization dynamics and electron-hole pair excitation. When the density
of photoexcited carriers becomes sufficiently large, collective plasma
oscillations in the electron gas emerge. Plasma formation is also accompanied
by highly efficient free-carrier absorption and severe spectral broadening of
the transmitted pulse. The last effect was first observed in Ref.Alfano  and
later it was also observed in gases, solids and liquids Couairon2007 ;
Wang2016 ; Dubietis2019 . Supercontinuum generation exhibits threshold
behavior and occurs when the incident laser power exceeds the critical power
for self-focusing Dubietis2017 , the characteristic broad pedestal on the blue
side of the spectrum of the transmitted pulse is often associated to
occurrence of breakdown Couairon2007 ; Gaeta2000 ; Liu2002 ; Lu2014 ;
Sudrie2002 ; Couairon2005 .

In this paper we report results of theoretical calculations of photoionization
and energy deposition by intense femtosecond laser pulses, which allows us to
determine the optical breakdown threshold in silicon. More specifically we
discuss scaling laws of the absorbed energy and photoelectron density as
functions of the laser fluence, and obtain the dependence of the threshold
fluence on the pulse length and laser polarization direction. We consider
self-broadening and self-phase modulation effects in the laser breakdown
plasma. The photoinduced modification of optical properties are presented in
terms of frequency-dependent reflectivity and absorption coefficients.
Softening of the Si-Si bonds and the emergent mechanical instability of the
diamond lattice are also discussed. Unless otherwise stated, atomic units are
used throughout this paper ($e=\hbar=m_{e}=1$).

##  II Results and discussion

Details of the theoretical approach and methodology can be found in Refs.
Lagomarsino2016 ; Apostolova2020  where we solve numerically the time-
dependent Schrodinger equation for covalent dielectrics interacting with
linearly polarized intense ultrashort laser pulse of near-infrared wavelength.
The time-dependent electric field of the laser pulse transmitted in bulk
silicon is

${\bf E}(t)={\bf E}_{{\rm ext}}(t)+{\bf E}_{{\rm ind}}(t),$ (1)

a superposition of the applied laser electric field ${\bf E}_{{\rm ext}}$ and
an induced electric field ${\bf E}_{{\rm ind}}(t)=-4\pi{\bf P}(t)$ as a result
of the polarization of the photoexcited solid. The applied electric field is
parameterized by a temporary Gaussian function

${\bf E}_{{\rm ext}}(t)={\bf e}e^{-\ln(4)t^{2}/\tau^{2}}F\cos\omega_{L}t$ (2)

where ${\bf e}$ is a unit vector in the direction of the laser polarization,
$\omega_{L}$ is the laser frequency corresponding to photon energy
$\hbar\omega_{L}$ = 1.55 eV, $\tau$ is the pulse length and $F$ is the
electric field strength related to the peak laser intensity by $I=F^{2}$ (in
atomic units). We determine the breakdown threshold from the fluence energy
dependence of the total electronic excitation energy (per unit cell)

$E_{abs}=\int dt{\bf J}(t)\cdot{\bf E}(t)$ (3)

where ${\bf J}(t)=d{\bf P}/dt$ is the photoinduced current.

###  II.1 Optical breakdown threshold of silicon - fluence and pulse length
dependence

As in our earlier investigation Apostolova2020 , the laser wavelength was set
to 800 nm, three different pulse durations $\tau$ = 30, 60 and 90 fs were
considered and two different laser polarization directions - [001] and [111]
were used in our calculation. Figs.1(a-b) show the dependence of the absorbed
energy on the energy fluence $\Phi=I\tau$ of the applied laser. The dependence
on the peak laser intensity is also shown in Figs.1 (c,d). The low fluence
regime is characterized by low conversion efficiency of the laser energy into
electronic excitation energy (less than 0.01 eV/atom). Dielectric breakdown
threshold is characterized by an abrupt step-like increase of the absorbed
energy by nearly 1 eV/atom. For above threshold fluences, energy absorption at
near infrared wavelengths increases slowly following the scaling law
$E_{abs}\sim\Phi^{1/2}$. The energy absorption curves shown in Fig.1 (c-d) are
parallel, electrons gain more energy for longer pulse duration. That is
because the probability for photoionization increases linearly with $\tau$, as
a result the absorbed energy also increases linearly with the increase of the
pulse length, such that electrons gain more energy by interacting with longer
pulses.

Fluence thresholds increase by approximately 1 J/cm2 for increasing pulse
duration by 30fs. Fig.2 (a,b) show the threshold fluence $\Phi_{th}$ and
threshold intensity $I_{th}$ as a function of $\tau$, respectively. We find
linear dependence $\Phi_{th}=\Phi_{0}+I_{\infty}\tau$, where $\Phi_{0}$ is the
threshold fluence extrapolated to infinitely short pulse duration,
$I_{\infty}$ is the asymptotic threshold intensity for continuous wave laser
irradiation ($\tau\rightarrow\infty$). Consequently the threshold intensity
decreases with the increase of the pulse length according to
$I_{th}=\Phi_{0}/\tau+I_{\infty}$. These two parameters depend on the laser
polarization direction: $I_{\infty}=3\times 10^{13}$ W/cm2, $\Phi_{0}=0.95$
J/cm2 for [001] polarization direction. $I_{\infty}=2.8\times 10^{13}$ W/cm2,
$\Phi_{0}=0.82$ J/cm2 for [111] polarization direction. The result agrees
qualitatively with the experimental results for damage thresholds in fused
silica reported in Refs.Lenzner1998 ; Chimier2011 .

Figs.3 (a-d) presents the photo-electron density and the absorbed laser energy
per Si atom as a function of the laser fluence. The low fluence regime is
characterized by perturbative three-photon absorption with super-linear
scaling trend of the photoelectron density $n_{c}\sim\Phi^{3}$. For higher
laser fluences $\Phi\geq 1$ J/cm2, a qualitative change of electron dynamics
occurs with sub-linear scaling law $n_{c}\sim\Phi^{1/3}$. Similar result is
found in the dependence of the absorbed energy on the laser energy fluence,
shown in Fig.3 (b,d). At low fluences, perturbative scaling trend is exhibited
with $E_{abs}\sim\Phi^{4}$ due to four-photon absorption. The sudden step-like
increase of laser energy deposition is associated with the onset of free-
carrier absorption and dielectric breakdown of Si. Free carriers promoted into
the conduction band behave like those in simple metal described by a Drude
model. Because their effective mass is large compared to the electron mass,
their screened plasma frequency
$\omega_{p}^{2}=4\pi(32/a_{0}^{3})e^{2}/(m^{\ast}\epsilon_{b}$) is in the
infrared range and nearly matches the laser oscillation frequency (here
$\epsilon_{b}\approx\epsilon_{0}$ is the background dielectric constant of the
photoexcited Si). At this point, resonant energy transfer occurs from the
laser pulse to the electrons and an optical breakdown occurs. The deposited
energy at threshold $E_{break}\approx$ 0.2 eV/atom is independent on the
parameters of the laser pulse and slightly exceeds the corresponding melting
temperature of silicon (1700 Kelvin). For above threshold fluences, the
electronic excitation energy $E_{abs}>$ 1 eV/atom exceeds a few times the
thermal melting threshold and a characteristic square-root scaling law
$E_{abs}\sim\Phi^{1/2}$ is exhibited. This sub-linear scaling trend of the
absorbed energy in the breakdown plasma is a result of balance between kinetic
and potential energy gain by electrons.

The carrier density and absorbed energy depend sensitively on the laser
polarization direction as shown in Fig.4 (a)-(e). The energy fluence threshold
is lower for a laser pulse linearly polarized along the [111] direction as
compared to the [001] direction because in the latter case the laser is
aligned with Si-Si bonds, which increases the probability for photoionization.
Fig.4 (a-e) demonstrate this laser polarization dependent modification of the
energy fluence threshold.

The density of conduction electron states (per unit cell) is plotted in Fig.5
for a laser linearly polarized along the [111] direction. For each pulse
duration the distribution is structured into two broad absorption peaks with
energies 2.5 and 4 eV above the valence band maximum, which are associated
with perturbative 3 and 4-photon transitions. The spectral weight of the
absorption peaks increases with the increase of the laser intensity. For near
threshold fluences, when substantial part of valence electrons are promoted
into the conduction band, the role of Pauli’s exclusion principle becomes
prominent: photoinization saturates. Depending on the pulse duration,
saturation is reached for laser intensity $I>7\times 10^{13}$ W/cm2 for 30fs
pulse, and $I>6\times 10^{13}$ W/cm2 for $\tau$ = 60 and 90fs, cf. also Fig.5.

###  II.2 Supercontinuum generation

The time evolution of the electric field of the transmitted pulse is shown in
Fig.6. In the low intensity regime dielectric response is exhibited 6(a-c):
the pulse exhibits temporary Gaussian profile with electric field ${\bf
E}_{{\rm ext}}/\epsilon_{0}$, where $\epsilon_{0}=12$ is the static bulk
dielectric constant of silicon. The field strength does not exceed 0.1 V/Å. In
the high intensity regime 6(d-f), and above threshold fluences, the field
strength reaches 0.2 V/Åon the rising edge of the pulse and exhibits an
exponential growth in time, the peak field strength is substantially delayed
relative to the applied laser and reaches 0.4 V/Å, next a steep edge is formed
in the back of the pulse and frequency up-chirp is exhibited. This frequency
upshift causes severe spectral broadening of the transmitted pulse, which
extends into the visible region and results in supercontinuum generation (cf.
below).

The time evolution of the electronic excitation energy is shown in Fig.7
(a-c). In the high fluence regime, the electronic excitation energy increases
in cumulative manner following the optical cycles of the driving pulse. The
laser energy irreversibly deposited to electrons exceeds the thermal melting
threshold after conclusion of the pulse. The energy transfer rate increases
rapidly with the increased peak intensity, in the above threshold regime more
than 1 eV per Si-atom is deposited over few-cycles prior to the pulse peak. At
that point metallic response is exhibited: the electric field is screened out
by free carriers, such that their number and excitation energy quickly
saturate.

The spectral profile of the transmitted pulse is also shown in Fig.8. For
below threshold fluences, the spectrum is centered around the laser wavelength
and exhibits Gaussian profile. For above threshold fluences in the breakdown
regime, a broad shoulder appears on the blue side of the central wavelength
and extends into the visible region. This effect is consistent with plasma
model of superbroadening proposed in Ref. Bloembergen1973 ; Yablonovitch1974 :
the generation of free carriers creates dense plasma in the wake of the pulse,
which decreases the refractive index and produces self-phase modulation in the
breakdown regime. As a result the laser frequency becomes time-dependent and
upshifts in the back of the pulse, which causes severe spectral broadening and
manifests in broad supercontinuum generation.

###  II.3 Optical properties of laser-excited silicon

The real and imaginary parts of the dielectric function of the photoexcited Si
are shown in Fig.9(a-b). The main feature in the response of the photoexcited
solid is the negative divergence of ${\rm Re}\epsilon$ in the infrared region
(cf. Fig.9(a)), as a result of plasma formation of free carriers. Fig.9(b)
shows that the imaginary part of the dielectric function becomes appreciable
below the direct absorption edge. The diminution of the main absorption peak
is accompanied by transfer of spectral weight below the absorption edge and
into the Drude peak near $\omega=0$. Similar results were found in
Ref.Sato2014p1 . An interesting feature of the response is that the dielectric
tensor has off-diagonal components in the laser-excited state, even though the
crystal exhibits cubic symmetry, cf. Fig.10 (a-b).

In Fig.11 the frequency-dependent absorption coefficient $\alpha(\omega))$
(left panel) and the normal-incidence reflectivity $R(\omega)$ (right panel)
of the photo-excited Si are also shown. Here the absorption coefficient is
defined by $\alpha=4\pi\kappa/\lambda$, where $\kappa$ is the extinction
coefficient of the laser-exicted Si and $\lambda$ is the wavelength in vacuum.

For below threshold fluences, the dielectric response of Si is exhibited. The
main absorption peak gradually decreases with increasing the laser intensity,
the position of the peak undergoes a progressive blue-shift. A strongly
absorbing state of silicon for near infrared laser wavelengths is produced in
the breakdown plasma, when the pumped solid becomes opaque. For such large
absorption coefficients $\alpha\sim 10^{5}$ cm-1 in the NIR region, the light
probes only a thin surface layer less than 1 $\mu$m in thickness. As a result,
the reflectivity of the excited Si becomes very sensitive to the surface
quality of the crystal (contamination and oxidation of the surface layer).

###  II.4 Bonding properties of photo-excited Si

Fig. 12(b-c) display the change of electron density distribution in the (110)
plane in Si after photo-excitation by intense 30fs laser pulse. For
comparison, Fig.12(a) shows the ground state valence charge density map
characterized by bond charges located in the middle of each Si-Si bond. The
local pseudo-potential model produces closed elliptic density contour lines
with long axis perpendicular to the bonding direction. Photo-ionization causes
redistribution of charge density, electrons which occupy bonding sites move
into the interstitial anti-binding regions, thereby reducing the amplitude of
bond charges, cf. Fig. 12(b-c). More electrons are displaced from bonds
parallel to the laser polarization direction ([111] direction) as compared to
bonds that are oriented perpendicularly ([1 $\bar{1}$ 1] direction). The bond
charges remain as centers of inversion of the diamond lattice, their amplitude
decreases and shape changes.

There is simple criteria for stability of the diamond lattice, cf. Heine1976 ;
Kudryashov2002 . The transverse acoustic (TA) phonons in Si are unstable when
the bond charge vanishes Martin1969 , the crystal no longer resists to shear
stress and melts. Because the bond charge in Si is $Z_{b}=2/\epsilon_{0}=1/6$
Phillips1968 , the creation of a single electron-hole pair destroys the bond
charge $Z_{b}$ in six Si-Si bonds, thus the excitation of electron-hole pairs
of number density $n_{c}$ softens the transverse acoustic phonon modes by

$\omega_{{\rm TA}}^{2}\rightarrow\omega_{{\rm
TA}}^{2}\left(1-\epsilon_{0}\frac{n_{c}}{n_{0}}\right)$ (4)

where $n_{0}=32/a_{0}^{3}$ is the average bulk electron density. The simple
estimation predicts that when 8 % of the valence electrons are promoted into
the conduction band, the diamond lattice becomes unstable and melts on
ultrashort time scale. This is consistent with both our result and more
detailed calculation based on phenomenological bond-charge model in Ref.
Biswas1982 .

Charge transfer from bonding sites results in filling of vacuum and
interstitial regions, which partially destroys the covalent bonding
interaction between the Si atoms. The freed electrons create an uncompensated
outward kinetic pressure, such that crystal is no longer in mechanical
equilibrium and expands uniformly. However this lattice expansion should be
slow on a ultrashort time scale (because the speed of sound in Si
$c_{s}=10^{4}$ m/s). Ref.Stampfli1990  modeled the elasticity of the photo-
excited Si using a phenomenological valence force field model with three
adjustable parameters - central force constant, bond-bending constant and a
non-equilibrium photo-induced pressure constant. They considered shear
distortions of the diamond lattice and showed that large internal pressure
makes negative the elastic shear modulus, softens the transverse acoustic
phonons and results in ultra-fast melting of silicon. Our independent
calculation of the change of macroscopic stress tensor of Si for fixed-ion
positions Nielsen1985 , shows that large uncompensated positive pressure
arises after photo-excitation and that more than 80 % of the photoinduced
pressure is due increase of the electronic kinetic energy. The electronic
kinetic pressure is evaluated from

$p=\frac{1}{3}\sum_{v,{\bf G}}\int_{{\rm BZ}}\frac{d^{3}{\bf
k}}{4\pi^{3}}\left(|c_{v,{\bf G}+{\bf k}}|^{2}-|c^{0}_{v,{\bf G}+{\bf
k}}|^{2}\right)({\bf G}+{\bf k})^{2}$ (5)

where ${\bf G}$ are the reciprocal-lattice wave-vectors, the quasi-momentum
${\bf k}$ extends over the first Brillouin zone, $v$ labels initially occupied
valence band states, $c^{0}_{v,{\bf G}+{\bf k}}$ and $c_{v,{\bf G}+{\bf k}}$
are the electronic momentum-space wave-functions before and after the
excitation. Fig.13 shows the dependence of the photo-induced pressure on the
laser intensity. For high level of electronic excitation $n_{c}\geq 10^{22}$
cm-3, the pressure $p\sim$ 10 GPa exceeds the elastic strength of Si. More
importantly, $p$ exhibits a characteristic sub-linear scaling trend with
$p\sim I^{5/9}$ and because $n_{c}\sim I^{1/3}$, the famous relation between
pressure and density $p\sim n_{c}^{5/3}$ for degenerate Fermi gas is obtained.
Thus the Pauli’s exclusion principle is prominent in the breakdown regime.
This result differs from the prediction of the phenomenological model in
Ref.Stampfli1990 , which obtains pressure increasing linearly with the
increase of conduction electron density. The breakdown plasma gives rise to
large Fermi pressure, that can not be compensated by electron-ion attraction.
Thus electron emission from the silicon surface occurs, leaving behind
positive holes. The charged Si surface may become electrostatically unstable,
and emit massive positively charged projectiles on a longer time scale.

##  III Conclusion

The optical breakdown threshold of crystalline silicon is obtained as a
function of the laser energy fluence, pulse duration and polarization
direction. The occurrence of breakdown is associated with plasma formation and
highly efficient free-carrier absorption. The transmitted pulse is distorted
and becomes subject to self-steepening and self-phase modulation, the
corresponding breakdown field strength is self-limited and does not exceed 0.4
V/Å. The severe spectral broadening of transmitted pulse results in generation
of supercontinuum extending into the visible spectral region. The optical
properties of Si undergo irreversible modification in the infrared spectral
region: both absorption and reflectivity of near infrared wavelengths are
greatly enhanced in the breakdown plasma. In this regime of high electronic
excitation, an unbalanced photo-induced degeneracy pressure results in
mechanical instability of the lattice, which may be considered as a precursor
to ultrafast melting and ablation of Si. We plan more detailed investigation
of the elasticity of the photoexcited solid for a follow-up paper.

## Acknowledgements

Support from the Bulgarian National Science Fund under Contract No. KP-06-KOST
is acknowledged. This research is based on work supported by the Air Force
Office of Scientific Research under award number FA9550-19-1-7003.

Figure 1: Fluence and intensity dependence of the the absorbed energy (in
eV/atom) in bulk silicon irradiated by pulsed laser of wavelength 800 nm.
Different symbols in (a-d) designate the laser pulse duration: 30fs (squares),
60 fs (upper triangles) and 90 fs (circles). The laser is linearly polarized
along the [001] crystal direction in Figs. (a),(c), and is polarized along the
[111] direction in Figs.(b) and (d).  Figure 2: Pulse length dependence of the
threshold fluence (a) and threshold intensity (b) after irradiation of silicon
by a pulsed laser with wavelength 800 nm. Different symbols in (a-b) designate
the laser polarization direction: [001] crystal direction (squares) and [111]
direction (upper triangles). Figure 3: Fluence dependence of the absorbed
energy per Si atom (in eV) (a),(c) and the conduction electron density (in
cm-3) (b),(d) in bulk silicon subjected to pulsed laser irradiation with near
infrared wavelength 800 nm. Different symbols in (a-d) designate the laser
pulse duration: 30fs (black squares), 60 fs (upper triangles) and 90 fs
(circles). The laser is linearly polarized along the [001] direction in
Fig.(a,b) and is polarized along the [111] direction in Fig. (c,d). Figure 4:
Fluence dependence of the conduction electron density (in cm-3) (a),(c),(e)
and the absorbed energy(in eV/atom) (b),(d),(f) in bulk silicon irradiated by
pulsed laser with near-infrared wavelength 800 nm. The pulse duration is: 30fs
in Figs.(a),(b), 60fs in Figs.(c),(d) and 90fs in Fgis.(e), (f). Different
symbols in (a-b) designate the laser polarization direction: [001] crystal
direction (squares) and [111] direction (upper triangles).  Figure 5: Density
of conduction electron states per unit cell (in eV-1) in bulk silicon after
irradiation by a pulsed laser with near-infrared wavelength 800 nm. Different
line types in Fig.(a-c) designate different laser intensities: In Fig. (a) the
pulse length is 30 fs and the intensity is $I=3.\times 10^{13}$ W/cm2 (solid
line) , 5.$\times 10^{13}$W/cm2 (dashed line), 6.$\times 10^{13}$W/cm2 (dotted
line), 7. $\times 10^{14}$W/cm2 (short-dashed line) and 8.0 $\times 10^{14}$
W/cm2 (short dotted line); Fig. (b) the pulse length is 60 fs and the
intensity is 3.$\times 10^{13}$W/cm2 (solid line), 4.5$\times 10^{13}$W/cm2
(dashed line), 5.$\times 10^{13}$ W/cm2 (dotted line), 6.$\times 10^{14}$W/cm2
(short-dashed line), and 7.$\times 10^{14}$W/cm2 (short dotted line); Fig. (c)
the pulse length is 90 fs and the intensity is 3.$\times 10^{13}$ W/cm2 (solid
line), 3.8$\times 10^{13}$W/cm2 (dashed line), 5.$\times 10^{13}$W/cm2 (dotted
line), 6.$\times 10^{14}$W/cm2 (short-dashed line), and 7.$\times
10^{14}$W/cm2 (short dotted line). In Fig. (a-c) the laser is linearly
polarized along the [111] crystal direction. Figure 6: Time evolution of the
electric field (in V/Å) in bulk silicon irradiated by pulsed laser with near-
infrared wavelength 800 nm. In Fig. (a) and Fig.(d) the pulse length is 30fs,
and the peak laser intensity is 2 $\times 10^{13}$ W/cm2 and $10^{14}$ W/cm2,
respectively. In Fig. (b) and Fig.(e) the pulse length is 60 fs and the laser
intensity is $3\times 10^{13}$ W/cm2 and 8 $\times 10^{13}$ W/cm2,
respectively. In Fig. (c) and Fig.(e) the pulse duration is 90fs and the peak
laser intensity is $2\times 10^{13}$ W/cm2 and $5\times 10^{13}$ W/cm2,
respectively. In Fig. (a-f), the laser is linearly polarized along the [001]
crystal direction.

Figure 7: Time evolution of the absorbed energy (in eV/atom) in bulk silicon
irradiated by pulsed laser with near-infrared wavelength 800 nm. The laser
energy fluence (in J/cm2)is indicated above each curve. The pulse duration is
30 fs in Fig.(a), 60 fs in Fig.(b) and 90 fs in Fig.(c). The laser is linearly
polarized along the [001] crystal direction.

Figure 8: Spectrum of the transmitted pulse in bulk silicon irradiated by
pulsed laser with near-infrared wavelength 800 nm. Different line types in
(a-b) designate different laser intensities: In Fig.(a) the pulse duration is
30, the peak laser intensity is 5.$\times 10^{13}$W/cm2 (solid line) and
$10^{14}$W/cm2 (dashed line). In Fig. (b) the pulse duration is 60 fs and the
peak laser intensity is 3$\times 10^{13}$W/cm2 (solid line) and 9$\times
10^{13}$W/cm2 (dashed line) and in Fgi.(c)the pulse length is 90 fs, the peak
laser intensity is 3$\times 10^{13}$ W/cm2 (solid line) and 9$\times
10^{13}$W/cm2 (dashed line). In Fgi.(a-c) the applied laser is linearly
polarized along the [001] crystal direction. Figure 9: (a) Real and (b)
imaginary part of the dielectric function of photoexcited Si irradiated by
90fs laser pulse with wavelength 800 nm, linearly polarized along the [111]
axis. In Fig. (a,b) the different line types designate different laser
intensities: 2$\times 10^{13}$ W/cm2 (solid line), 3.6$\times 10^{13}$W/cm2
(dashed line), 4$\times 10^{13}$W/cm2 (dotted line) and 5$\times 10^{13}$W/cm2
(dashed-dot line). Figure 10: Real (a) and imaginary part (b) of the
dielectric function of photoexcited Si irradiated by 90fs laser pulse with
wavelength 800 nm and peak laser intensity 4$\times 10^{13}$ W/cm2. The laser
is linearly polarized along the [001] axis (solid line) and along the [111]
axis (dotted line). Figure 11: Left panel- frequency dependent absorption
coefficient of Si induced by 90fs laser pulse with wavelength 800 nm, linearly
polarized along the [111] crystal direction. The peak intensity of the applied
laser is 2$\times 10^{13}$ W/cm2 in (a), 3$\times 10^{13}$W/cm2 in Fig.(b),
3.4$\times 10^{13}$W/cm2 in Fig.(c) and 4$\times 10^{13}$W/cm2 Fig.(d). Right
panel- the corresponding frequency-dependent reflectivity of photoexcited
silicon. The laser intensity is 2$\times 10^{13}$W/cm2 in Fgi.(e), 3$\times
10^{13}$ W/cm2 in Fig.(f), 3.4$\times 10^{13}$W/cm2 in Fig.(g) and 4$\times
10^{13}$W/cm2 in Fig.(f).  Figure 12: Contour plots of the electron number
density in the (110) plane in silicon. Fig.(a) - ground state electron
density. Fig.(b) and Fig.(c) display the local density difference map $n({\bf
r})-n_{0}({\bf r})$ after irradiation of Si by 30fs laser pulse with
wavelength 800 nm and peak laser intensity 1014 W/cm2. In Fig.(b) and Fig.(c)
the laser is linearly polarized along the [111] and [001] axis, respectively.
In (b,c), negative contours (charge depletion) are indicated by the dotted
lines and positive contours (charge accumulation) by the solid lines. The
contour interval is 0.01 a.u. in (a), and it is 0.0005 a.u. in Fig.(b) and
Fig.(c).  Figure 13: Laser intensity dependence of the photoinduced pressure
(in kbar) in crystalline Si after the irradiation by pulsed laser of near-
infrared wavelength 800 nm. The pulse duration is 30 fs in Fig.(a), 60 fs in
Fig.(b) and 90fs in Fig.(c). The laser is linearly polarized along the [001]
axis (squares) and along the [111] axis (circles).

## References

  * (1) N. Bloembergen, Journal of Nonlinear Optical Physics and Materials 6, 377 (1997).
  * (2) D. Du , X. Liu, G. Korn, J. Squier and G. Mourou, Appl. Phys. Lett. 64 3071 (1994).
  * (3) B. C. Stuart, M. D. Veit, S. Herman, A. M. Rubenchik, B. W. Shore, M. D. Perry, Phys. Rev. B 53 , 1749 (1996).
  * (4) W. Kautek, J. Krüger, M. Lenzner, S. Sartania, Ch. Spielmann, F. Krausz, Appl. Phys. Lett. 69, 3146 (1996).
  * (5) H. Varel, D. Ashkenasi, A. Rosenfeld, R. Herrmann, F. Noack and E. E. B. Campbell, Appl. Phys. A 62, 293 (1996).
  * (6) M. Lenzner, J. Kruger, S. Sartania, Z. Cheng, C. Spielmann, G. Mourou, W. Kautek, and F. Krausz, Phys. Rev. Lett. 80 , 4076 (1998).
  * (7) P. P. Pronko, P. A. VanRompay, C. Horvath, F. Loesel, T. Juhasz, X. Liu, and G. Mourou, Phys. Rev. B 58, 2387 (1998).
  * (8) M. Li, S. Menon, J. P. Nibarger and G. N. Gibson, Phys. Rev. Lett. 82, 2394 (1999).
  * (9) H. M. van Driel, Phys. Rev. B 35, 8166 (1987).
  * (10) T. Apostolova, Y. Hahn, J. Appl. Phys. 88, 1024 (2000).
  * (11) C. B. Schaffer, A. Brodeur and E. Mazur, Meas. Sci. Technol. 12, pp. 1784-1794 (2001).
  * (12) T. Apostolova, D. H. Huang, P. M. Alsing, J. McIver, D. A. Cardimona, Phys. Rev. B. 66, 075208 (2002).
  * (13) E. G. Gamaly, A. V. Rode,B. Luther-Davies and V. T. Tikhonchuk, Physiscs of Plasmas 9, 949 (2002).
  * (14) S. I. Kudryashov, A. Joglekar, G. Mourou, A. A. Ionin V. D. Zvorykin and A. J. Hunt, in Proc. of SPIE Vol. 6733, p.67332H (2007).
  * (15) P. P. Rajeev, M. Gertsvolf, P. B. Corkum, D. M. Rayner, Phys. Rev. Lett. 102, 083001 (2009).
  * (16) B. Chimier, O. Uteza, N. Sanner, M. Sentis, T. Itina, P. Lassonde, F. Legare, F. Vidal, and J. C. Kieffer, Phys. Rev. B 84, 094104 (2011).
  * (17) O. Utéza, N. Sanner, B. Chimier, A. Brocas, N. Varkentina, M. Sentis, P. Lassonde, F. Legare, and J. C. Kieffer, Appl. Phys. A 105, 131 (2011).
  * (18) D. Grojo, S. Leyder, P. Delaporte, W. Marine, M. Sentis, and O. Uteza, Phys. Rev. B 88, 195135 (2013).
  * (19) M. Lebugle, N. Sanner, O. Utéza, and M. Sentis, Appl. Phys. A 114, 129 (2014).
  * (20) S. A. Sato, K. Yabana, Y. Shinohara, T. Otobe, K.-M. Lee and G. F. Bertsch, Phys. Rev. B 92, 205413 (2015).
  * (21) P. A. Zhokhov and A. M. Zheltikov, Scientific reports 8, 1824 (2018).
  * (22) C. V. Shank, R. Yen, and C. Hirlimann, Phys. Rev. Lett. 51, 900, (1983).
  * (23) D. Hulin, M. Combescot, J. Bok, A. Migus, J. Y. Vinet, A. Antonetti, Phys. Rev. Lett. 52, 1998 (1984).
  * (24) H. W. K. Tom, G. D. Aumiller, C. H. Brito-Cruz, Phys. Rev. Lett. 60, 1438 (1988).
  * (25) K. Sokolowski-Tinten, J. Bialkowski, D. von der Linde, Phys. Rev. B 51, 14186 (1995).
  * (26) A. Rousse, C. Rischel, S. Fourmaux, I. Uschmann, S. Sebban, G. Grillon, P. Balcou, E. Förster, J. P. Geindre, P. Audebert, and J. C. Gauthier, Nature, 410, 65 (2001).
  * (27) S. K. Sundaram, and E. Mazur, Nat. Mater. 1, 217, (2002).
  * (28) M. Harb, R. Ernstorfer, C. T. Hebeisen, G. Sciaini, W. Peng, T. Dartigalongue, M. A. Eriksson, M. G. Lagally, S. G. Kruglik, R. J. Dwayne Miller, Phys. Rev. Lett. 100, 155504 (2008)
  * (29) P. Stampfli and K. H. Bennemann, Phys. Rev. B 42, 7163 (1990)
  * (30) P. Stampfli and K. H. Bennemann, Phys. Rev. B 46, 10686 (1992)
  * (31) P. Stampfli and K. H. Bennemann, Phys. Rev. B 49, 7299 (1999).
  * (32) P. L. Silvestrelli, A. Alavi, M. Parrinello and D. Frenkel Phys. Rev. Lett. 77, 3149 (1996).
  * (33) E. G. Gamaly, Phys. Rep. 508, 91 (2011).
  * (34) P. Balling and J. Schou, Rep. Prog. Phys. 76, 036502 (2013).
  * (35) S. A. Sato, K. Yabana, Y. Shinohara, T. Otobe, and G. F. Bertsch, Phys. Rev. B 89, 064304 (2014).
  * (36) S. A. Sato, Y. Shinohara, T. Otobe, and K. Yabana, Phys. Rev. B 90, 174303 (2014).
  * (37) M. Schultze, K. Ramasesha, C. D. Pemmaraju, S. Sato, D. Whitmore, A. Gandman, J. Prell, L. Borja, D. Prendergast, K. Yabana, D. Neumark, S. Leone, Science, 346, 1348 (2014).
  * (38) T. Apostolova, B. Obreshkov, I. Gnilitskyi, Appl. Surf. Sci. 519, 146087 (2020).
  * (39) T. Otobe, Phys. Rev. Applied 13, 024062 (2020).
  * (40) S. Yamada and K. Yabana, Phys. Rev. B 103, 155426 (2021).
  * (41) R. R. Alfano and S. L. Shapiro, Phys. Rev. Lett. 24, 584 (1970).
  * (42) A. Couairon, A. Mysyrowicz, Phys. Rep. 441, 47 (2007).
  * (43) Q. Z. Wang, P. P. Ho, and R. R. Alfano _I_ n The Supercontinuum Laser Source, pp. 33-100. Springer, New York, NY, (2016).
  * (44) A. Dubietis and A. Couairon. _U_ ltrafast supercontinuum generation in transparent solid-state media, Springer (2019).
  * (45) A. Dubietis, G. Tamošauskas, R. Šuminas, V. Jukna, and A. Couairon, Lith. J. Phys. 57, no. 3 (2017).
  * (46) A. L. Gaeta, Phys. Rev. Lett. 84, 3582 (2000).
  * (47) W. Liu, S. Petit, A. Becker, N. Akozbek, C. M. Bowden, and S. L. Chin, Opt. Commun. 202, 189 (2002).
  * (48) C. H. Lu, Y. J. Tsou, H. Y. Chen, B. H. Chen, Y. C. Cheng, S. D. Yang, M. C. Chen, C. C. Hsu, and A. H. Kung, Optica 1, 400 (2014)
  * (49) L. Sudrie, A. Couairon, M. Franco, B. Lamouroux, B. Prade, S. Tzortzakis, and A. Mysyrowicz, Phys. Rev. Lett. 89, 186601 (2002).
  * (50) A. Couairon, L. Sudrie, M. Franco, B. Prade, and A. Mysyrowicz, Phys. Rev. B 71, 125435 (2005).
  * (51) S. Lagomarsino, S. Sciortino, B. Obreshkov, T. Apostolova, C. Corsi, M. Bellini, E. Berdermann and C. J. Schmidt, Phys. Rev. B 93, 085128 (2016).
  * (52) N. Bloembergen, Opt. Comm. 8, 285 (1973).
  * (53) E. Yablonovitch, Self-phase modulation and short-pulse generation from laser-breakdown plasmas, Phys. Rev. A 10 , pp. 1888-1895 (1974).
  * (54) V. Heine, J. A. Van Vechten, Phys. Rev. B 13, 1622 (1976).
  * (55) S. I. Kudryashov and V. I. Emelyanov, J. Expt. Th. Phys. 94, pp. 94-107 (2002).
  * (56) R. M. Martin, Phys. Rev. 186, 871 (1969).
  * (57) J. C. Phillips, Phys. Rev. 166, 832 (1968).
  * (58) R. Biswas, V. Ambegaokar, Phys. Rev. B 26, 1980 (1982).
  * (59) O. H. Nielsen, R. M. Martin, Phys. Rev. B 32, 3792 (1985).


#  Cell A* for Navigation of Unmanned Aerial Vehicles in Partially-known
Environments

Wenjian Hao, Rongyao Wang, Alexander Krolicki, Yiqiang Han
Department of Mechanical Engineering
Clemson University
Clemson, SC 29634 USA
Graduate Researcher<EMAIL_ADDRESS>Researcher,
<EMAIL_ADDRESS>Researcher<EMAIL_ADDRESS>Assistant Professor<EMAIL_ADDRESS>Author

###### Abstract

Proper path planning is the first step of robust and efficient autonomous
navigation for mobile robots. Meanwhile, it is still challenging for robots to
work in a complex environment without complete prior information. This paper
presents an extension to the A* search algorithm and its variants to make the
path planning stable with less computational burden while handling long-
distance tasks. The implemented algorithm is capable of online searching for a
collision-free and smooth path when heading to the defined goal position. This
paper deploys the algorithm on the autonomous drone platform and implements it
on a remote control car for algorithm efficiency validation.

##  1 Introduction

Unmanned Aerial Vehicle (UAV) is a critical component in the ecosystem of
autonomous systems designs. UAVs have unparalleled flexibility in the three-
dimensional space, creating novel applications that are impossible for
autonomous ground vehicles (AGV) to accomplish. These extra degrees of freedom
are coupled with added complexity to the dynamics of the system. UAVs navigate
primarily in unstructured environments but in some cases they must navigate
structured environments where a clear path is not always present. As a vital
step for autonomous navigation, path planning must schedule a smooth path and
avoid obstacles in real-time for UAVs and AGVs alike. In order to solve the
UAVs path planning problem, various algorithms have been proposed to obtain an
energy-efficient and obstacle-free route. The existing algorithms can be
classified into the following categories: sample-based algorithms, heuristic-
based algorithms, and intelligent bionic algorithms. For more algorithms,
details refer to [1, 2, 3].

This study focuses on the application of an extended A* algorithm which is
also a heuristic-based algorithm [4]. Since the sampling-based approach cannot
guarantee an optimal solution, we do not want to put our mission-critical path
planning tasks up to chance or luck of the draw. The heuristic-based
algorithms have been widely implemented in various fields such as autonomous
driving [5], autonomous parking [6], etc. The standard A* algorithm’s basic
idea depends on a global or local heuristics table to track the best candidate
node with the lowest computed heuristic cost. The algorithm works efficiently
in a reasonably scaled map with known obstacles. To achieve dynamic re-
planning in the vicinity of dynamic obstacles in partially known environments,
several A* variants have been employed such as D*, Anytime A* (ARA) [7, 8],
etc.

Although A* and its variants have been successfully applied to various mobile
robots in two-dimensional structured environments, practical use in UAVs path
planning still suffers from the increased degrees of freedom as well as a lack
of explicit prior knowledge of the environment. For example, micro-
scale/small-scale UAVs in urban environments bring about challenges associated
with complexity from the highly structured environment, noisy signals due to
EMI, and a need for energy efficiency due to battery constraints for flights
covering long distances of unexplored areas. In many practical use cases of
UAVs, GPS signals are limited, especially for indoor or rural areas.
Therefore, we need to mainly rely on data collected from distance sensors to
obtain a local map, rather than an accurate plan available for AGVs on roads
with high-definition maps. Besides, 3-D environments are often too large to
achieve accurate global planning and inevitable bias due to position
measurement error, execution bias, etc. and adds complexity.

To overcome the aforementioned problems, this study proposes two extensions to
the A* and its variants deployed in three-dimensional space. First, in real-
time, we prioritize a local cubic search cell as the drone’s sampling space
and pick nodes with the lowest cost at each layer to get a ’cell path.’ In
this way, we can re-plan a cell path at each time step to reduce the execution
bias without accumulating error. Secondly, to avoid aggressive path smoothing,
we introduce an extra cost for the distance between the current position and a
3-D straight line connecting the starting point and goal point. In this way,
we can drive the drone close to the theoretic shortest path and make long-
distance path planning in an energy- and time-efficient manner.

The experimentation of this algorithm is divided into two parts. The first
part is conducted within a customized 2-D simulation environment. For this
simple 2-D example, the dynamical system model is defined as a bicycle model
for an AGV, which is essentially a simplified version of the Ackermann
steering model for low-speed application [9]. The 2-D simulation was
implemented in MATLAB[10] for demonstration purposes. In this case, 2D testing
is implemented based on our modification of hybrid A* for a non-holonomic
vehicle. The second example is based on a more complex UAV model simulated in
a 3-D GAZEBO virtual environment [11]. It utilizes a simultaneous localization
and mapping (SLAM) algorithm for UAVs localization [12] and a stereo camera to
create a point cloud map of unknown obstacles. A schematic of the environment
is shown in Figure 1.

This paper is organized as follows: Section 2 briefly introduces the A*
algorithm and its major variants and presents details of the extensions we
propose. Section 3 shows the algorithm verification in the 2-D world by
simulating the algorithm on an RC car with constraints in a 2-D Matlab
simulator. Section 4 exhibits the experiment execution on the drone platform
in both simulator and the real world with the results. Finally, section 5
concludes the paper.

Figure 1: Path planning in 2-D and 3-D scenarios

##  2 Extension of Hybrid A* Algorithm

In this section, we briefly present the A* algorithm with its variants, and
our extended version of the A* algorithm: Cell A*. We will also discuss the
different treatments when applying the proposed algorithm to AGV with a
bicycle model on a 2-D plane and a UAV with holonomic behavior in 3-D space.

###  2.1 Hybrid A*

In real-world applications, standard A* has a disadvantage that the planned
path could be infeasible for robots with differing kinematic constraints. In
practical application, Hybrid A* is proposed to include non-holonomic
constraints during expanding path search, as shown in unstructured outdoor
environments [5] and unknown semi-structured environments [13]. A standard
hybrid A* flow chart is shown in Figure 2.

Figure 2: Schematics of hybrid A* in steps

In Figure 2 the heuristic function $H$ in Equation 1 maintains a tracking
record of two distances. The first part is $f_{t}$, the euclidean distance
from the current state to the goal state. $g_{t}$ is the euclidean distance
from the current state to the start state. $w_{t}$ is the hyperparameter with
a value between 0 and 1, which balances the priority between chasing the final
goal and finding an optimal path.

$\begin{gathered}f_{t}=||P_{t}(x_{t},y_{t},z_{t})-P_{goal}(x_{g},y_{g},z_{g})||_{2}\\\
g_{t}=||P_{t}(x_{t},y_{t},z_{t})-P_{0}(x_{0},y_{0},z_{0})||_{2}\\\
H_{t}=f_{t}+w_{g}\cdot g_{t}\end{gathered}$ (1)

Hybrid A* generates a path based on the global map as well as a predefined
non-holonomic sampler with a sampling direction aligned with the goal
position. Before implementing hybrid A* planning, all environment information
is assumed to be available to maintain a rapid search on a small scale map.
Yet in the real-word application, planning across the entire map requires a
large amount of memory storage and high latency cost. Furthermore, the
environment information may not be fully known. Hence, for practical use
purposes, when the robot is in a semi-known environment, it is necessary to
navigate and re-plan its path by constructing a map and generating a path
simultaneously [13], which can, in turn, lead to tens of thousands of nodes of
information.

With sufficient computing power, the Hybrid A* is guaranteed to achieve
optimal path planning for a given problem. However, in applications for
micro/small scale autonomous vehicles, the computing power is limited, and
real-time online searching capability is more of a priority. Therefore, we
seek an extension to the current algorithm for an efficient algorithm suitable
for edge computing on a mobile robot.

###  2.2 Extension of A* algorithm

To maintain hybrid A* exploratory capability while shrink the minimum
necessary size of the map and number of nodes stored, this paper proposes an
extension of the hybrid A* algorithm that can dissect the global hybrid A*
into local hybrid A* within a cell that has an adaptive size and ensures the
smoothness of the path.

The proposed algorithm shares the same backbone of Hybrid A* while expanding
the adaptive searching capability into the 3-D neighboring nodes. The details
of the algorithm consist of the following steps: Considering the focus of this
paper is the drones navigation, the algorithm should sample its neighbor nodes
in all the directions.

  1. 1.

Generate Search Nodes:
For predictive nodes sampling, we choose a cubic cell for the drone:
$Cell_{t}(cellsize,gridsize)$, where the $cellsize$ is the number of nodes on
each cell side, the $gridsize$ is the length of each grid in the cell. In this
way, the cell can cover all the directions and expand the range of searching
if necessary. The search cell ($Cell^{search}_{t}$) of position $P_{t}$ at
time $t$ would be $Cell_{t}(cellsize,gridsize)+P_{t}$

  2. 2.

Evaluate the costs of all the valid nodes in $Cell^{search}_{t}$:
To reduce the computation of the evaluation, we compare the adjacency cells
with the real-time updated local occupancy map. We only take the valid cell
points, which is not in the obstacles point cloud set. The cost contains two
parts. One is an A*-like heuristic cost $H_{t}$, which tracks the distance
between the current position ($P_{t}$) and goal position ($P_{goal}$); the
other cost $L_{t}$ is the distance between the current position ($P_{goal}$)
and the Euclidean shortest path which connects the initial position and the
goal position. While $H_{t}$ pushes the mobile robots close to the goal
position, $L_{t}$ ensures a theoretic shortest path of the robots and
compensates the bias generated by the path following, especially for drones in
complex outdoor environments. The total cost ($J$) is shown in Equation 2,
$w1$ and $w2$ is weights of $H_{t}$ and $L_{t}$.

$\begin{gathered}J=w_{1}H_{t}+w_{2}L_{t}\end{gathered}$ (2)

Where $H_{t}$ is in Equation 3

$\begin{gathered}H_{t}=||P_{t}(x_{t},y_{t},z_{t})-P_{goal}(x_{g},y_{g},z_{g})||_{2}\end{gathered}$
(3)

To track the cost term $L_{t}$, we require a partially known map which enables
us initialize the starting point $P_{0}(x_{0},y_{0},z_{0})$, goal position
$P_{goal}(x_{g},y_{g},z_{g})$ and obtain current position
$P_{t}(x_{t},y_{t},z_{t})$ with Simultaneous Localization and Mapping (SLAM)
or GPS. The cost term $L_{t}$ is defined in Equation 4.

$\begin{gathered}3D\
line:\frac{x-x_{0}}{x_{g}-x_{0}}=\frac{y-y_{0}}{y_{g}-y_{0}}=\frac{z-z_{0}}{z_{g}-z_{0}}\\\
\\\ L_{t}=|distance(P_{t},line)|\end{gathered}$ (4)

  3. 3.

Verify Obstacles:

    1. (a)

If front obstacles on path are out of avoidance range: Set
$Cell^{search}_{t}(3,bigstep)$, where $cellsize=3$, $gridsize=bigstep$ means
it is a one layer cell with 3 nodes on each edge so as to enable a fast path
planning. Meanwhile, $J=w_{1}H_{t}+w_{2}L_{t}$

    2. (b)

Else, if it is able to pass the obstacle: $J=w_{1}H_{t}+w_{2}L_{t}$
if not: $J=w_{1}H_{t}-w_{2}L_{t}$

    3. (c)

We allow the sign of $L_{t}$ term to change to balance the need to find the
shortest path based on global heuristics and encourage local replanning to
avoid obstacles along the shortest path. To avoid unforeseen obstacles or
dynamic obstacles along the global path, the robot may need to deviate from
the shortest path and pivot to exploration in neighboring cells. We allow the
sign of $L_{t}$ in heuristic function to reverse to encourage robot navigating
away from the global planned path.

  4. 4.

Search an optimal path in the search cell $Cell^{search}_{t}$:
Set $cellpath$ list, in each layer of the search cell, select the node with
the lowest cost ($N_{t}=argmin_{node}J(Cell^{search}_{t}$), then pass the
lowest cost node to the cell path set ($cellpath\leftarrow N_{t}$), finally
choose the first $cellsize-1$ nodes as the cell path for safety sake at time
step $t$. The planned path ($path_{t}$) at time step $t$ would be:
$path_{t}\leftarrow cellpath_{1...(cellsize-1)}$.

  5. 5.

Check for Constraints:
Integrate constraints into the cost function for some mobile robots with
kinematic limitations like RC cars steering.

  6. 6.

Check for Arrival: If the distance between the current position and goal
position is acceptable, STOP.

The summary of the algorithm is shown below:

Input: Current state $P_{t}(x_{t},y_{t},z_{t})$ $\leftarrow$ SLAM algorithm

Obstacles set $\leftarrow$ points cloud from stereo camera

Output: Planned trajectory in the search cell at each time step $t$

Initialization

  1. 1.

Set start state $P_{0}(x_{0},y_{0},z_{0})$ and goal state
$P_{goal}(x_{g},y_{g},z_{g})$

  2. 2.

Set starting node = $node_{0}$

  3. 3.

Heuristic function H($P_{t},P_{goal}$), L($P_{goal},P_{0}$)

  4. 4.

Define nodal cost function $J(Cell_{t}^{search},Obstacles)$

for  _node in $Cell_{t}^{search}$_ do

SafeToPass = node $\notin$ $Obstacles\ set$? 1:-1

$J_{node}=w_{1}H_{t}+w_{2}L_{t}\times SafeToPass$

end for

while  _goal reached == false_ do

if  _obstacles on path $\notin$ search range_ then

$Cell_{t}^{search}$ = $genSearchCell_{node}(3,gridsize)$ \+ $P_{t}$

$Cell\ Path\leftarrow node=argmin_{node}J(Cell_{t}^{search},Obstacles)$

else

for  _i in range ( $cellsize$)_ do

$Cell_{t}^{search}=genSearchCell_{node}(i,gridsize)+P_{t}$

$Cell\ Path\leftarrow node=argmin_{node}J(Cell_{t}^{search},Obstacles)$

end for

end if

end while

return _Cell Path_

Algorithm 1 Cell A* for UAVs in 3D environment

##  3 Algorithm Planning Efficiency Validation in 2-D Environment

In this section, we validate the Cell A* in a MATLAB 2-D simulation
environment [10]. The entire validation consists of three parts. Firstly we
implement it with non-holonomic constraints to mimic the real car driving
scenario. We then compare the Cell A* with constraints with the local hybrid
A*, finally, we compare the proposed Cell A* algorithm with the D* lite
algorithm, which can ensure the shortest planned path in terms of the planned
path length and planning time. To get a fair benchmark, we do not optimize any
of these comparison algorithms.

###  3.1 Cell A* Validation with kinematic constraints

Unlike a drone’s motion being holonomic, a car-like robot’s motion is
restricted by Ackermann steering kinematic model. Considering a minimum
functional model [9], the vehicle states can be represented by
$(x,y,\theta,d)$, where $x,y$ is coordinate position on a 2D plane, $\theta$
is the current yaw angle of the car with respect to x coordinate axis, and $d$
represents the direction of robot’s motion. For simplicity, The rear axle is
considered as the center of rotation. A schematic of the model setup is shown
in Figure 3.

Figure 3: Ackermann Sampler Demonstration

The expansion of the search cell is governed by two parameters, $l$, and $d$.
$l$ control the distance range of the sampling cell edge in the 2D plane. A
larger sampling distance $L$ means sampler can reach further on a 2D plane,
thus increase exploratory performance. Change of $d$ can change the growth
direction of sampling cells. Similar to the purpose of reversing the sign of
$L_{t}$ in nodal cost function, changing the sign of $d$ also serves the
purpose of encouraging the exploration behavior of the agent when an obstacle
is detected along the global planned path. During the simulation, to further
promote exploration behavior, the sampling steps within one search are also
increased if local/dynamic obstacles are detected.

To accommodate the need for local replanning, which comes with the constrain
of nonholonomic behavior, $g_{t}$ in the nodal cost function $J$ no longer
tracks distance to the global start state, but rather the distance to the
start state of the current sampling start position. In this case, the global
start position got abandoned after first sampling, hence reducing the
possibility of circling back to the original state when encountering
obstacles.

A simple 2-D simulation demonstration is shown in Figure 4, where black blocks
represent different obstacles, two red points are the starting point and the
endpoint. As exhibited in the figure, when the robot encounters an obstacle
that leads to zero new neighbors from the current sampler setting, sampler
changes searching direction to move away from the direct path between the
start state and goal state. Once front obstacles are absent, the planner
return to searching the shortest path to the final goal from current state as
expected.

Figure 4: Cell A* real-time planning process

###  3.2 Experiment for Comparison with Local Hybrid A*

To compare the extension algorithm with the standard hybrid A* algorithm in
the same benchmark, we controlled the searching steps of both cell A* and
hybrid A* to compare the path difference between them under similar nodes
visited. The experiment is implemented in MATLAB.

As shown in Figure 5, once hybrid A* was constrained within a small step of
search, due to the limit of sampling depth, it failed to escape from obstacles
in certain occasions like scenario with goal state $(15,50)$. The robot
started to move back to the trapped corner if the local nodes sampler fail to
expand itself into a safe-to-pass area. Since cell A* adjusts sampling
distance according to the current sampling situation, it can increase sampling
distance to search deeper into the environment, hence increase the possibility
of finding a path away from obstacles.

Figure 5: Cell A* (Blue) vs. Local Hybrid A* (Green) for four different
destinations

Table 1 details how these two methods perform when they both head to the same
goal position. "Failed" in Local Hybrid A* means robot ended up in a dead
corner and failed to escape. From table, it can be inferred that cell A* needs
less computation, and it is more stable compared to the local hybrid A*
algorithm.

Table 1: Cell A* vs. Local hybrid A* Algorithm | Goal state | Time (sec) | Path Length (m) | Average Nodes
---|---|---|---|---
| (6, 60) | 40.53 | 80.01 | 30.5
Cell A* | (15, 50) | 46.33 | 91.93 | 27.79
| (30, 43) | 43.67 | 87.08 | 27.96
| (6, 60) | 47.43 | 86.5 | 31.1
Local Hybrid A* | (15, 50) | Failed | Failed | Failed
| (30, 43) | 53.293 | 91.5 | 29.56

###  3.3 Experiment for Comparison with D* lite

To further measure the navigation performance of Cell A*, a comparison between
D* Lite and Cell A* in discrete state-space is implemented in MATLAB. In this
test, the robot has 8 degrees of freedom of motion in a 2-D plane. D* Lite is
built based on dynamic programming; therefore, it can return the shortest path
with any given 2-D map information at the cost of high time complexity. This
comparison focuses on the overall traveling distance and planning time. The
result is shown in Figure 6.

As shown in the picture, Cell A* has a different path pattern compared to the
path generated from D* Lite. To further analyze the difference, a detailed
statistic is shown in Table 2.

Figure 6: Cell A* (Blue) vs. D* Lite (Green) in discrete map Table 2: Cell A* vs. D* lite Algorithm | Goal state | Total Search Time (sec) | Path Length (number of nodes)
---|---|---|---
| (15, 25) | 0.12 | 115
Cell A* | (41, 25) | 0.13 | 230
| (55, 65) | 0.25 | 403
| (15, 25) | 1.45 | 110
D* lite | (41, 25) | 6.12 | 206
| (55, 65) | 2.07 | 331

As shown in Table 2, Cell’s A* has much less computational time due to the
strictly constrained search depth. However, the path nodes are not
significantly larger than optimal solution from D* Lite, which indicates
cell’s A* can avoid some redundant searches from D* Lite.

##  4 Algorithm Implement for 3D Case and Results

This section exhibits the results of the algorithm implemented on a drone
platform in a 3D environment.

###  4.1 Drones Execution Set-up

Before we deployed the algorithm in the real world, we tested it in a gazebo
virtual environment built by GAAS group [14]. The simulation is implemented in
the robotics operating system (ROS) [15]. Localization is enabled by SLAM via
a stereo camera sensor. For the real-world implementation, we modified a TBS
Discovery drone frame to accommodate the on-board sensors, Pixracer PX4
autopilot, and NVIDIA TX2 single board computer. The TX2 build environment is
configured using jetpack 4.2, which runs on ubuntu 18.04. The ROS build
version is melodic. The communication between the high-level (e.g. Local pose
updates from the navigation algorithm) and low-level control (flight
controller) is done over MAVROS. The PX4 controller [16] takes care of the
path following part in both simulation and real-world experiments in real-
time. The above codes are all written with python language, and the physical
drone we use is exhibited in Figure 7.

Figure 7: Modified TBS Discovery development platform.

On-board sensors shown in Figure 7 include a downward facing LiDAR for
altitude measurement. The introduction of LiDAR provides accurate altitude
reading by fusing the measurements with the on-board barometer readings. The
on-board sensing module accurately tracks the drone’s planar velocity using an
optical flow sensor and is fused with IMU data to more accurately measure the
speed of the drone. The forward-facing stereo camera (Zed 2 RGBD camera) is
used to generate the local map data.

###  4.2 Experiment and Results

####  4.2.1 Experiment for Obstacles Avoidance

In this experiment, the algorithm is tested on its ability to keep away from a
static wall obstacle. The drone starts form position $(0,0,2.5)$ and heads to
the goal position point $(7,0,2.5)$. The results is shown in Figure 8, where
the red points are the recorded trajectory.

Figure 8: Obstacles avoidance test

As shown in Figure 8, the drone correctly avoid the wall while it is within
its defined obstacles avoidance range.

####  4.2.2 Experiment for Long Distance Task

In this experiment, the algorithm is tested on its ability for a long-distance
path planning task. The drone takes off at position $(0,0,6)$ and the goal
position is point $(45,-6,8)$, the euclidean distance between these two points
is $45.44m$. A planned route scheduled with the algorithm is exhibited in
Figure 9.

Figure 9: Long distance task test

A clean path is shown in Figure 9 that the drones avoid some trees at first
and then flies straightly to the defined position, which is the balcony of the
building. The top subplot in Figure 9 shows the detected occupancy grids in a
3-D space. The 3D representation of the world is mostly sparse since we design
Cell A* with efficiency as the top priority. We demonstrate that online
trajectory planning can be achieved for an unknown environment.

###  4.3 Experiment Results

Table 3 shows the results achieved by the two drone deployment experiments. It
can be inferred that the drone simulation with the proposed algorithm can
effectively handle unknown obstacles like a wall or foliage on trees during
navigating in a semi-unknown environment.

Table 3: Experiment Data Experiment | Task’s Euclidean Distance (m) | Time (sec) | Exploit Nodes
---|---|---|---
Obstacle avoidance | 7 | 22.2319 | 13000
Long Distance | 45.44 | 98.6619 | 91000

##  5 Conclusion and Future Work

This paper presents an extension work of the Hybrid A* algorithm to enable
UAVs run in a 3-D world or a car-like robot in a 2-D world without any prior
information of the environment. By adjusting the nodal cost function, we allow
the cost function to be a linear combination of the euclidean distance between
the current location to the goal location and the distance between the current
and the shortest euclidean path between the initial position and the goal
position. In contrast, the sign of the second term is adaptive, depending on
the obstacle in front of the planned path. The proposed Cell A* requires much
less computing costs for the on-board computer. The algorithm requires at
least one order of magnitude less sampling nodes than the Hybrid A* algorithm
at each time step during online replanning. By taking advantage of the
computational saving of the Cell A*, the drone in the simulation performs
navigation with high efficiency not only in terms of obtaining a collision-
free and smooth trajectory path but also saving much less time compared to the
standard A* algorithm and its variants.

## References

  * [1] Zhangjie Fu, Jingnan Yu, Guowu Xie, Yiming Chen, and Yuanhang Mao.  A heuristic evolutionary algorithm of uav path planning.  Wireless Communications and Mobile Computing, 2018:2851964:1–2851964:11, 2018.
  * [2] Tanaka Kitampura.  3-d path planning in a dynamic environment using an octree and an artificial potential field.  IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 1–7, 1995.
  * [3] Liang Yang, Juntong Qi, J. Xiao, and Xia Yong.  A literature review of uav 3d path planning.  In Proceeding of the 11th World Congress on Intelligent Control and Automation, pages 2376–2381, 2014.
  * [4] David Ferguson, Maxim Likhachev, and Anthony (Tony) Stentz.  A guide to heuristic-based path planning.  In Proceedings of the International Workshop on Planning under Uncertainty for Autonomous Systems, International Conference on Automated Planning and Scheduling (ICAPS), June 2005.
  * [5] J. Petereit, T. Emter, C. W. Frey, T. Kopfstedt, and A. Beutel.  Application of hybrid a* to an autonomous mobile robot for path planning in unstructured outdoor environments.  In ROBOTIK 2012; 7th German Conference on Robotics, pages 1–6, 2012\.
  * [6] S. Sedighi, D. Nguyen, and K. Kuhnert.  Guided hybrid a-star path planning algorithm for valet parking applications.  In 2019 5th International Conference on Control, Automation and Robotics (ICCAR), pages 570–575, 2019.
  * [7] Anthony Stentz.  The focussed d* algorithm for real-time replanning.  In In Proceedings of the International Joint Conference on Artificial Intelligence, pages 1652–1659, 1995.
  * [8] Maxim Likhachev, David I. Ferguson, Geoffrey J. Gordon, Anthony Stentz, and Sebastian Thrun.  Anytime dynamic a*: An anytime, replanning algorithm.  In ICAPS, 2005.
  * [9] W. Wang, T. Hsu, and T. Wu.  The improved pure pursuit algorithm for autonomous driving advanced system.  In 2017 IEEE 10th International Workshop on Computational Intelligence and Applications (IWCIA), pages 33–38, 2017.
  * [10] The Mathworks, Inc., Natick, Massachusetts.  MATLAB version 9.7.0.1247435 (R2019b), 2019.
  * [11] N. Koenig and A. Howard.  Design and use paradigms for gazebo, an open-source multi-robot simulator.  In 2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (IEEE Cat. No.04CH37566), volume 3, pages 2149–2154 vol.3, 2004.
  * [12] R. Latif and A. Saddik.  Slam algorithms implementation in a uav, based on a heterogeneous system: A survey.  In 2019 4th World Conference on Complex Systems (WCCS), pages 1–6, 2019.
  * [13] Dmitri Dolgov, Sebastian Thrun, Michael Montemerlo, and James Diebel.  Path planning for autonomous vehicles in unknown semi-structured environments.  I. J. Robotic Res., 29:485–501, 04 2010.
  * [14] Tao Wu Cyanine gi, Hanyang Wang.  Generalized autonomy aviation system.  Github, 2019.
  * [15] Shaun Edwards and Chris Lewis.  Ros-industrial: applying the robot operating system (ros) to industrial applications.  In IEEE Int. Conference on Robotics and Automation, ECHORD Workshop, 2012.
  * [16] E. A. Niit and W. J. Smit.  Integration of model reference adaptive control (mrac) with px4 firmware for quadcopters.  In 2017 24th International Conference on Mechatronics and Machine Vision in Practice (M2VIP), pages 1–6, 2017.