# Universality in odd-even harmonic generation
and application in terahertz waveform sampling

Doan-An Trieu  Computational Physics Key Laboratory K002, Department of
Physics,
Ho Chi Minh City University of Education, Ho Chi Minh City 72711, Vietnam
Ngoc-Loan Phan<EMAIL_ADDRESS>Computational Physics Key Laboratory
K002, Department of Physics,
Ho Chi Minh City University of Education, Ho Chi Minh City 72711, Vietnam
Quan-Hao Truong  Computational Physics Key Laboratory K002, Department of
Physics,
Ho Chi Minh City University of Education, Ho Chi Minh City 72711, Vietnam
Hien T. Nguyen  Tay Nguyen University, Buon Ma Thuot City 63161, Vietnam  Cam-
Tu Le  Atomic Molecular and Optical Physics Research Group, Advanced Institute
of Materials Science, Ton Duc Thang University, Ho Chi Minh City 72915,
Vietnam  Faculty of Applied Sciences, Ton Duc Thang University, Ho Chi Minh
City 72915, Vietnam  DinhDuy Vu  Computational Physics Key Laboratory K002,
Department of Physics,
Ho Chi Minh City University of Education, Ho Chi Minh City 72711, Vietnam
Van-Hoang Le<EMAIL_ADDRESS>Computational Physics Key Laboratory K002,
Department of Physics,
Ho Chi Minh City University of Education, Ho Chi Minh City 72711, Vietnam

###### Abstract

Odd-even harmonic generation contains the system’s dynamical symmetry
breaking, and its decoding is desirable not only in understanding physics but
also for applications of dynamics extraction. In this work, we discover a
simple universal relation between the odd-even harmonics and the asymmetry of
the THz-assisted laser-atomic system – atoms in a primary midinfrared laser
pulse combined with a THz laser. We demonstrate numerically and then
analytically derive the dependence of the harmonic even-to-odd ratio on the
THz electric field. Notably, the functional form of this dependency reveals a
universal scaling independent of the parameters of both the primary pulse and
atomic target. This universality inspires us to propose a pump-probe scheme
for THz waveform sampling from the even-to-odd ratio, accessible from
conventional compact setups.

Introduction — Recently, high-order harmonic generation (HHG) from an
asymmetric laser-target system has gained great attention since it reveals
deeper structures in HHG spectra [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. The most
prominent feature in such symmetry-breaking systems is the emergence of even
harmonics, which apparently separate them from symmetric counterparts, which
emit only odd harmonics [12]. Therefore, the odd-even pattern in HHG spectra
must encode the dynamical asymmetry information of the laser-target system [1,
2, 3, 7, 8, 4, 5, 6, 9, 10], and thus a universal relation between them can be
helpful in many extraction applications.

One route to break the spatial-temporal symmetry is adding a weaker external
static electric field to the primary multi-cycle laser pulse [13, 14, 15, 16,
17, 18]. However, to yield a visible effect, the static field needs to be
impractically strong ($\sim$MV/cm). Nevertheless, with recent developments of
powerful terahertz (THz) sources, a quasi-static strong electric field can be
engineered into the HHG process [19, 20, 21, 22, 23, 24, 25, 18, 26, 27]. This
THz-assisted HHG has been studied mostly in the harmonic conversion efficiency
or the plateau structure [19, 20, 21, 22, 23] but less in other notable
aspects such as odd-even harmonic spectra. Some initial estimates have been
made [14, 26], but a direct relation between the odd-even pattern and the
asymmetricity of the THz-assisted laser-atomic system has not yet been
established. This relation is not only of intellectual interest but can be
critical to many applications such as extracting quantum dynamics inside
atoms/molecules, manipulating electron trajectories on the attosecond time
scale, or in THz waveform sampling, which has been gaining attention recently
[3, 22, 28, 29, 18].

Figure 1: (a) Pump-probe scheme for sampling THz waveforms with sequential
steps: (i) experimental setup with the THz pulse (pump) and a primary mid-IR
laser (probe) both shinning on a gas jet; (ii) measurements of intensities of
even and two adjacent odd harmonics at the cutoff region with pump-probe time
delays; (iii) extracting step with the first computing the even-to-odd ratio
from HHG signals and then extracting THz electric field using the universal
rule. (b) Comparison between the extracted (dotted curve) and benchmark (solid
curve) THz pulses. Five-cycle trapezoidal primary pulses with the intensity of
$2.5\times 10^{14}$ W/cm2 and wavelength of 2000 nm are used as probe lasers.

In this work, we accomplish two goals. First, we thoroughly study the response
of the odd-even HHG to the variation of the THz electric field in a THz-
assisted laser-atomic system. The purpose is to achieve a quantitative
connection between the measurable harmonic even-to-odd ratio and the symmetry-
breaking THz electric field. Interestingly, we discover a universal and simple
rule for the even-to-odd ratio as a function of the scaled THz electric field.
This finding stimulates us to proceed with the second goal of proposing a
general method for sampling the THz pulse waveform. Our method can retrieve
the time-resolved waveform with high accuracy from the harmonic even-to-odd
ratio by an algorithm shown in Fig. 1. In short, we collect the HHG signals
emitted from the gas jet inside the intersection between a primary midinfrared
(mid-IR) linearly polarized laser and the THz pulse. Our algorithm allows a
direct reconstruction of the instantaneous THz strength through the even-to-
odd ratio of harmonics of the HHG readout. By tuning the pump-probe time
delay, we can profile the entire THz pulse. In this vein, we note that
waveform detection is another important area of THz science that is
intensively investigated in addition to its generation [28, 29, 18]. Within
currently available methods, active matter has to be carefully chosen; thus,
finding a new temporal detection scheme free from external parameters is
experimentally meaningful.

Universal response of even-to-odd ratio to THz electric field — We first look
thoroughly into the response of resolved odd-even harmonics emitted from a
hydrogen atom in the combined linearly polarized primary mid-IR laser pulse
and the slow-varying THz field, as shown in Fig. 2(a). Here, the HHG data are
calculated by numerically solving the time-dependent Schrödinger equation
(TDSE method) [30]. Since electron wave packets are mostly spread along the
laser polarization direction, the one-dimensional model is sufficient for the
HHG study as it requires low computational overhead but is quantitatively
comparable to the three-dimensional model [31, 32]. We consider an example of
a mid-IR pulse with a 10-cycle trapezoidal envelope, the intensity of
$2.5\times 10^{14}$ W/cm2, and the wavelength of 2000 nm. It is clear that
when the added THz field $E_{T}$ is much weaker than the laser field
[$<10^{-5}$ a.u. (51 kV/cm)], the spectrum is pure-odd harmonics [12]. Upon
increasing the THz strength, the space-time symmetry (half-period-time
translation combined with spatial inversion) of the original system is
violated, leading to the emergence of even harmonics.

For a more focused picture, Fig. 2(b) shows the response of harmonic intensity
with the THz strength for the two selected even (500${\mathrm{th}}$) and odd
(501${\mathrm{st}}$) harmonics at the cutoff, denoted respectively as H500 and
H501. Their behavior can be visually partitioned into two regions. (i) For the
first region, the intensity of the even harmonics gradually increases, while
that of the odd harmonics is almost unchanged. (ii) After the first
intersection at $E_{T}\approx 4\times 10^{-5}$ a.u. ($\sim 0.2$ MV/cm) that
indicates the equal presence of the even and odd harmonics, the second region
begins with the intensity fluctuation of the odd and even harmonics around
their average values. Notably, the odd and even harmonics at some specific THz
field values alternatively undergo deep minimums corresponding to pure-even or
pure-odd spectra, as illustrated in Fig. 2(a).

Figure 2: Response to changing of THz electric field for (a) Intensity of
resolved odd-even harmonic spectra, (b) Selected odd (H501) and even (H500)
harmonics at the cutoff, and (c) Harmonic even-to-odd ratio. The primary mid-
IR laser pulse with the same parameters as in Fig. 1 is used for the HHG
process; the THz field with a frequency of 1.3 THz (231 $\mu$m) is used; the
color bar in Panel (a) decodes the HHG intensity in arbitrary units; the
dotted horizontal line in Panel (c) shows the unity.  Figure 3: Dependence of
harmonic even-to-odd ratio on pure (left panels) and scaled (right panels) THz
electric field for a hydrogen atom in various THz-assisted primary laser
pulses with (a) different intensities (fixed wavelength of 2000 nm) and (b)
different wavelengths (fixed intensity of 2.5$I_{0}$), where
$I_{0}=1\times$1014 W/cm2. The grey-shaded areas in Panels (c) and (d)
highlight the region of stable even-to-odd ratio.

Since the absolute harmonic efficiency is strongly affected by the laser
field, we introduce a dimensionless quantity by taking the intensity ratio
between the even harmonic and the average of the two adjacent odd ones, which
is referred to as the even-to-odd ratio. The variation of this quantity with
respect to the THz field is shown in Fig. 2(c). However, taking the intensity
ratio does not entirely cancel the effect of the primary laser still on the
THz-dependent even-to-odd ratio, as can be seen in the sensitivity to
different laser intensities and wavelengths shown in Figs. 3(a) and 3(b).
Specifically, the sensitivity to laser wavelength [Fig. 3(b)] is much more
severe than laser intensity [Fig. 3(a)]. Remarkably, we find that by rescaling
the THz electric field $E_{T}$ by the factor $E_{0}/\omega_{0}^{3}$ as

$\gamma=\dfrac{E_{0}}{\omega_{0}^{3}}E_{T},$ (1)

where $E_{0}$ and $\omega_{0}$ are the primary laser’s peak amplitude and
carrier frequency, the data collapse quantitatively within the range
$\gamma\gtrsim 0.1$ (see Fig. 3). Here, the region $\gamma$ below $0.1$
suffers from significant numerical noise. We also numerically verify the data
collapse with respect to the duration (5, 10, 15, 20 cycles) of the primary
laser pulse and, more interestingly, with different active atomic targets (H,
He, Ne, Ar) [30]. Furthermore, the universal response as a function $\gamma$
is observed not only for harmonics at the cutoff, but also for those below the
cutoff with an additional condition of good phase matching in HHG experiments
[30].

Analytical formula — The universal relation between the harmonic even-to-odd
ratio and scaled THz electric field observed numerically above motivates us to
uncover its underlying physics. We show that it is not accidental but can be
proven rigorously by the quantum-orbit theory [33]. Within this framework, the
generation of harmonic order $N$ is the coherent interference of attosecond
bursts emitted each half of an optical cycle:
$H(N\omega_{0})\approx|D_{1}e^{-i\Phi_{1}}+D_{2}e^{-i\Phi_{2}}|^{2}$, where
$D_{1(2)}$ and $\Phi_{1(2)}$ are their amplitude and phase. Therefore, the
even-to-odd ratio is mostly controlled by the relative intensity
$|D_{2}|^{2}/|D_{1}|^{2}$ and phase difference $\Delta\Phi$ between attosecond
bursts, which in turn is associated with the difference of quasi-classical
action $\Delta S$ of free electrons in the laser electric field as
$\Delta\Phi=N\pi-\mathrm{Re}(\Delta S)$. Adding a perturbative THz field does
not change significantly the intensity of attosecond bursts, i.e.,
$|D_{2}|^{2}/|D_{1}|^{2}\approx 1$. However, it distorts the quasi-classical
action by $S_{T}=C\gamma$ to the first order, leading to the change in the
phase of attosecond bursts, which in turn modifies the even-to-odd ratio as

$\eta=\tan^{2}\left(C\gamma\right).$ (2)

Here, $C=2\sin\theta~{}(\Delta\theta\cos\Delta\theta-\sin\Delta\theta)$ with
$\theta=\omega_{0}(t_{r}+t_{i})/2$ and
$\Delta\theta=\omega_{0}(t_{r}-t_{i})/2$ is a real dimensionless coefficient
depending on the instants of harmonic ionization $t_{i}$ and recombination
$t_{r}$. We emphasize that $C$ is a constant for each harmonic energy, and
does not contain any parameters of the laser or target, either explicitly or
implicitly. Specifically for harmonics at the cutoff, $C\approx 2.558$. See
more detailed calculations in Supplementary [30].

The analytical formula (2) demonstrates a direct connection between the even-
to-odd ratio $\eta$ (a normalized quantity characterizing the asymmetry of
measurable output) and the dimensionless scaled THz electric field $\gamma$
(the normalized symmetry-breaking factor). Most importantly, Eq. (2) is free
from parameters of the primary laser pulse parameters (except implicitly in
$\gamma$) and the atomic target, implying the universality of the response of
harmonic even-to-odd ratio to the scaled THz field. Besides, this analytical
expression also shows a periodic modulation of the harmonic even-to-odd ratio
with the period of $\pi/C$. It alternatively undergoes maxima and minima,
generating instants of the pure-even and pure-odd harmonic spectra. Based on
the relation (2), we further refer to the dimensionless quantity $\gamma$ as a
parameter describing the asymmetric degree of the laser-target system and call
it the asymmetry parameter.

Figure 4: Response of harmonic even-to-odd ratio to a wide range of scaled THz
fields calculated by the analytical formula $\eta=\tan^{2}(2.558\,\gamma)$
(black solid curve) and numerical TDSE method (dashed and dotted color curves)
for harmonics at the cutoff using different primary lasers enclosed in the
legend. The grey-, cream-, and mauve-shaded areas cover three regions with
different underlying physics mechanisms.

Figure 4 shows an overall visualization of the even-to-odd ratio versus the
asymmetry parameter $\gamma$ for the harmonics at the cutoff. This
visualization is based on the data obtained in two different ways, using
analytical relations and direct numerical calculations by the TDSE method.
Comparing the results of the two methods reveals three regimes with different
levels of analytic-numeric agreement, implying various physical mechanisms.
(i) For a small asymmetry parameter $0.1\lesssim\gamma\lesssim 0.6$ (gray-
shaded area), the analytical formula quantitatively predicts the even-to-odd
ratio obtained from numerically solving the TDSE with various primary laser
parameters and atomic targets. The agreement is expected within perturbative
strengths of the THz field because it only modifies the electron quasi-
classical motion in the continuum energy region but does not affect the
ionization and recombination steps. (ii) For higher asymmetric parameters
$0.6\lesssim\gamma\lesssim 3.5\pi/C\equiv 4.3$ (cream-shaded area), the
relation (2) fails to match the magnitude. Moreover, the disparity caused by
the laser’s parameters becomes prominent, unlike the perturbative regime.
However, the reversal points (from $>1$ to $<1$, and vice versa) still fall
into the sequence of $(k+0.5)\pi/2C$ with $k\in\mathbb{N}$. The main reason is
moderately strong THz field participates in the ionization step beyond
distorting the electron quasi-classical motion. This induces the imbalance
between adjacent attosecond bursts, thus reducing the interference contrast
and causing irregular modulation of the even-to-odd ratio. (iii) With
$\gamma\gtrsim 4.3$ (mauve-shaded area), the numerical even-to-odd ratio
becomes highly disordered, and the classical description (2) fails to predict
both the magnitude and spacing between the reversal points. Here, the intense
THz field modifies ionized electrons’ travel time compared to the field-free
case. Additionally, intense THz fields dominate over the primary laser field
in driving electron trajectories, thus altering the plateau structure of HHG
spectra.

In short, universality manifests itself in two different aspects: the even-to-
odd magnitude itself in the first region and its main oscillation frequency in
the second region. We note that the predictive power of Eq. (2) is optimized
[the (i) regime] if the primary laser’s parameters vary in an appropriate
working range (intensity within $[1.0-4.0]\times 10^{14}$ W/cm2 and wavelength
longer than 1200 nm (mid-IR laser)) [30]. In fact, the primary laser field
$E_{0}$ must be high enough not only to ensure a low controlled ratio
$E_{T}/E_{0}$ to avoid deforming the HHG plateau structure, but also not
exceed the intensity saturation. Also, the wavelength should be long so that
the asymmetry parameter $\gamma$ can fall into the stable range.

Application in THz waveform sampling — With laser parameters that optimize the
validity of Eq. (2), we propose a pump-probe THz waveform sampling scheme,
illustrated in Fig. 1(a). Here, the THz pulse becomes the pump pulse, being
probed by a delayed mid-IR laser. The two pulses meet inside a gas jet,
triggering the HHG process. The odd-even HHG [Fig. 1(ii)] and, consequently,
the even-to-odd ratio [Fig. 1(iii)] are recorded at each time delay. With the
analytical formula (2), the waveform of the THz pulse can be easily
reconstructed [Fig. 1(iii)].

Figure 1 shows an example of the waveform measurement for the THz pulse
$E_{T}(t)=E_{T0}\exp\left(-\omega_{T}^{2}t^{2}/36\pi^{2}\right)\sin\omega_{T}t$
with $E_{T}=257$ kV/cm and frequency $\omega_{T}=1.3$ THz, inspired by the
pulse recently reported in [34]. The five-cycle trapezoidal primary pulses
with $2.5\times 10^{14}$ W/cm2 intensity and 2000 nm wavelength are used as
probe lasers. The time resolution of the THz waveform sampling is related to
the probe pulse duration, which is about 20 fs for the five-cycle laser pulse
with three cycles in the flat part. Figure 1(b) demonstrates the validity of
our proposed procedure. It indicates a good consistency between the extracted
THz waveform and the input data. We have also examined and affirmed the
validity of the proposed method in sampling THz pulses with complicated
waveforms implying a broad frequency band (not shown). We note that the
carrier-envelope phase of the detected THz pulse might be flipped by $\pi$
since the proposed method can extract the magnitude only but not its sign.

Using the optimized universal rule can detect the THz electric field within a
wide THz field range of about $[20,2000]$ kV/cm imposed by the working range
of the probing laser pulses [30]. We emphasize that the detectable range can
be expanded if we use the suboptimal regime [regime (ii) of Fig. 4] of
$\gamma$, looking only at its stable spacing between reversal points. We leave
the details in future work.

In conclusion, we have both numerically and analytically demonstrated the
universal dependence of the even-to-odd ratio on the scaled THz electric field
(asymmetry parameter of the system). The approach to derive the universal rule
in this Letter can be generalized to other asymmetric laser-target systems
that may be meaningful in probing atomic/molecular quantum dynamics,
controlling electron dynamics within attosecond time scale, or extracting
asymmetric factors of laser-target systems.

Based on this universal rule, we have proposed a pump-probe method for THz
waveform sampling using the even-to-odd ratio, which is measurable within
current compact laser setups. Unlike the previous methods for THz detection
involving electronic or optical excitation of targets under THz pulse, our
proposed method works on the perturbative regime in which the THz field only
affects the dynamics of quasi-free electrons in the continuum energy region
leading to modulation of the even-to-odd ratio but does not directly interact
with the targets. This independence of the targets and probe laser parameters
makes the method feasible for detecting a wide range of THz electric fields.

###### Acknowledgements.

Acknowledgments -  This work was funded by Vingroup and supported by Vingroup
Innovation Foundation (VINIF) under project code VINIF.2021.DA00031. The
calculations were executed by the high-performance computing cluster at Ho Chi
Minh City University of Education, Vietnam.

## References

  * Frumker  _et al._ [2012] E. Frumker, N. Kajumba, J. B. Bertrand, H. J. Wörner, C. T. Hebeisen, P. Hockett, M. Spanner, S. Patchkovskii, G. G. Paulus, D. M. Villeneuve, A. Naumov, and P. B. Corkum, Phys. Rev. Lett. 109, 233904 (2012).
  * Chen  _et al._ [2013] Y. J. Chen, L. B. Fu, and J. Liu, Phys. Rev. Lett. 111, 073902 (2013).
  * Pedatzur  _et al._ [2015] O. Pedatzur, G. Orenstein, V. Serbinenko, H. Soifer, B. Bruner, A. Uzan, D. Brambila, A. Harvey, L. Torlina, F. Morales, O. Smirnova, and N. Dudovich, Nat. Phys. 11, 815 (2015).
  * Kraus  _et al._ [2015] P. Kraus, O. I. Tolstikhin, D. Baykusheva, A. Rupenyan, J. Schneider, C. Z. Bisgaard, T. Morishita, F. Jensen, L. B. Madsen, and H. J. Wörner, Nat. Commun. 6, 7039 (2015).
  * Chen  _et al._ [2019] C. Chen, F.-M. Guo, Y.-J. Yang, Y.-J. Chen, and S.-P. Yang, EPL (Europhysics Letters) 127, 34004 (2019).
  * Li  _et al._ [2019] W. Y. Li, R. H. Xu, X. J. Xie, and Y. J. Chen, Phys. Rev. A 100, 043421 (2019).
  * Yin  _et al._ [2021] H. Yin, D. Liu, and F. Zeng, Appl. Phys. Lett. 119, 151105 (2021).
  * Nguyen  _et al._ [2022a] H. T. Nguyen, K.-N. H. Nguyen, N.-L. Phan, C.-T. Le, D. Vu, L.-P. Tran, and V.-H. Le, Phys. Rev. A 105, 023106 (2022a).
  * He  _et al._ [2022] L. He, S. Sun, P. Lan, Y. He, B. Wang, P. Wang, X. Zhu, L. Li, W. Cao, P. Lu, and C. D. Lin, Nat. Commun. 13, 4595 (2022).
  * Nguyen  _et al._ [2022b] K.-N. H. Nguyen, N.-L. Phan, C.-T. Le, D. Vu, and V.-H. Le, Phys. Rev. A 106, 063108 (2022b).
  * Wang  _et al._ [2017] X. Wang, A.-T. Le, Z. Zhou, H. Wei, and C. D. Lin, Phys. Rev. A 96, 023424 (2017).
  * Ben-Tal  _et al._ [1993] N. Ben-Tal, N. Moiseyev, and A. Beswick, J. Phys. B: At. Mol. Opt. Phys. 26, 3017 (1993).
  * Bao and Starace [1996] M.-Q. Bao and A. F. Starace, Phys. Rev. A 53, R3723 (1996).
  * Wang  _et al._ [1998] B. Wang, X. Li, and P. Fu, J. Phys. B: At. Mol. Opt. 31, 1961 (1998).
  * Hong  _et al._ [2007] W. Hong, P. Lu, W. Cao, P. Lan, and X. Wang, J. Phys. B: At. Mol. Opt. 40, 2321 (2007).
  * Zhao  _et al._ [2011] G. Zhao, X. Guo, T. Shao, and K. Xue, New J. Phys. 13, 093035 (2011).
  * Le  _et al._ [2022] C.-T. Le, N.-L. Phan, D. D. Vu, C. Ngo, and V.-H. Le, Phys. Chem. Chem. Phys. 24, 6053 (2022).
  * Silaev  _et al._ [2022a] A. A. Silaev, A. A. Romanov, and N. V. Vvedenskii, Opt. Lett. 47, 4664 (2022a).
  * Hong  _et al._ [2009] W. Hong, P. Lu, P. Lan, Q. Zhang, and X. Wang, Opt. Express 17, 5139 (2009).
  * Balogh  _et al._ [2011] E. Balogh, K. Kovacs, P. Dombi, J. A. Fulop, G. Farkas, J. Hebling, V. Tosa, and K. Varju, Phys. Rev. A 84, 023806 (2011).
  * Jia  _et al._ [2015] C. Jia, J. Wang, Q.-Y. Li, F.-M. Guo, J.-G. Chen, S.-L. Zeng, and Y.-J. Yang, Opt. Express 23, 32222 (2015).
  * Ge  _et al._ [2015] X.-L. Ge, H. Du, J. Guo, and X.-S. Liu, Opt. Express 23, 8837 (2015).
  * Babushkin  _et al._ [2022] I. Babushkin, A. Demircan, U. Morgner, and A. Savel’ev, Phys. Rev. A 106, 013115 (2022).
  * Milošević [2022] D. B. Milošević, Opt. Lett. 47, 1669 (2022).
  * Brennecke  _et al._ [2022] S. Brennecke, M. Ranke, A. Dimitriou, S. Walther, M. J. Prandolini, M. Lein, and U. Frühling, Phys. Rev. Lett. 129, 213202 (2022).
  * Silaev  _et al._ [2022b] A. A. Silaev, A. A. Romanov, and N. V. Vvedenskii, J. Phys. Conf. Ser. 2249, 012004 (2022b).
  * Mao  _et al._ [2022] W. Mao, A. Rubio, and S. A. Sato, Phys. Rev. B 106, 024313 (2022).
  * Zhang  _et al._ [2017] X. C. Zhang, A. Shkurinov, and Y. Zhang, Nat. Photonics 11, 16 (2017).
  * Zhang  _et al._ [2021] Y. Zhang, K. Li, and H. Zhao, Front. Optoelectron. 14, 4 (2021).
  * [30] See Supplemental Materials for Theoretical model and numerical methods, Independent of even-to-odd ratio on driving laser pulse parameters and atomic targets, Analytical relation between even-to-odd ratio and THz electric field, Universality of even-to-odd ratio for harmonics below cutoff, and Working range.
  * Chirilă  _et al._ [2010] C. C. Chirilă, I. Dreissigacker, E. V. van der Zwan, and M. Lein, Phys. Rev. A 81, 033412 (2010).
  * Majorosi  _et al._ [2018] S. Majorosi, M. G. Benedict, and A. Czirják, Phys. Rev. A 98, 023401 (2018).
  * Salières  _et al._ [2001] P. Salières, B. Carré, L. Le Déroff, F. Grasbon, G. G. Paulus, H. Walther, R. Kopold, W. Becker, D. B. Milošević, A. Sanpera, and M. Lewenstein, Science 292, 902 (2001).
  * Shalaby  _et al._ [2016] M. Shalaby, C. Vicario, K. Thirupugalmani, S. Brahadeeswaran, and C. P. Hauri, Opt. Lett. 41, 1777 (2016).


# Harmonious Sampling for Mobile Manipulation Planning

Mincheul Kang1, Donghyuk Kim1 and Sung-Eui Yoon2 1Mincheul Kang
<EMAIL_ADDRESS>and 1Donghyuk Kim<EMAIL_ADDRESS>are
with the School of Computing, and 2Sung-Eui Yoon (Corresponding author,
<EMAIL_ADDRESS>is with the Faculty of School of Computing, KAIST at
Daejeon, Korea 34141

###### Abstract

Mobile manipulation planning commonly adopts a decoupled approach that
performs planning separately on the base and the manipulator. While this
approach is fast, it can generate sub-optimal paths. Another direction is a
coupled approach jointly adjusting the base and manipulator in a high-
dimensional configuration space. This coupled approach addresses sub-
optimality and incompleteness of the decoupled approach, but has not been
widely used due to its excessive computational overhead. Given this trade-off
space, we present a simple, yet effective mobile manipulation sampling method,
harmonious sampling, to perform the coupled approach mainly in difficult
regions, where we need to simultaneously maneuver the base and the
manipulator. Our method identifies such difficult regions through a low-
dimensional base space by utilizing a reachability map given the target end-
effector pose and narrow passage detected by generalized Voronoi diagram. For
the rest of simple regions, we perform sampling mainly on the base
configurations with a predefined joint configuration, accelerating the
planning process. We compare our method with the decoupled and coupled
approaches in six different problems with varying difficulty. Our method shows
meaningful improvements experimentally in terms of time to find an initial
solution (up to 5.6 times faster) and final solution cost (up to 17$\%$ lower)
over the decoupled approach, especially in difficult scenes with narrow space.
We also demonstrate these benefits with a real, mobile Hubo robot.

##  I INTRODUCTION

Mobile manipulation planning has been attracted much attention in recent years
since the rise of demands for autonomous indoor services, e.g., households
tasks or delivery services. Mobile manipulators usually equipped with wheels
on a base for the mobility perform various manipulation tasks, e.g., a pick-
and-place problem, with robotic arms.

For the mobile manipulation planning problem, a target end-effector pose is
given. Given the pose, we first compute its base configuration by utilizing a
reachability inversion [1] and then set joints of the manipulator based on
inverse kinematics for computing a goal configuration in the joint
configuration space (C-space) consisting of the base and manipulator.
Unfortunately, this joint C-space can be high-dimensional, and thus tends to
lead a slow planning process.

Therefore, a decoupled approach that plans separately for the base and the
manipulator has been commonly used in practice, since this decoupled approach
works relatively fast thanks to planning in separated and thus lower-
dimensional spaces [2]. Nonetheless, this decoupled approach may find sub-
optimal results or let the manipulator to place in a narrow space due to the
fixed base used for manipulation planning. This is mainly caused by its
decoupled planning and sampling that does not consider the joint planning
between the base and manipulator.

Figure 1:  This figure shows a sequence of the mobile Hubo robot grasping the
yellow beverage. The heat map on the bottom right represents a sample density
in terms of the 2D floor projected from samples generated from our harmonious
sampling in the joint C-space consisting of the base and the manipulator.
Samples are distributed intensively near the goal configuration, resulting in
efficient and effective exploration of solution paths.

Main contributions. In this work, we present a simple, yet effective mobile
manipulation sampling method, harmonious sampling, for optimal mobile
manipulation planning, to efficiently explore a high-dimensional C-space
consisting of both base and manipulator (Fig. 1). Our harmonious sampling
adjusts the sampling space for the base and the manipulator, and guides to
sample more on difficult regions where we need to simultaneously change the
configurations of both base and manipulator (Sec. IV-B). We identify such
regions based on a low-dimensional base space, especially, manipulation
regions that are computed by considering the target end-effector pose and
obstacles creating narrow passages (Sec. IV-A). For the rest of regions that
are guided from the base regions, we sample configurations of the base body,
while the manipulator has a predefined configuration. We also suggest a
region-specific k-nearest neighbor search for effectively connecting samples
generated differently by these two different regions (Sec. IV-C).

To validate our method, we compare our method with the coupled and decoupled
approaches for the base and manipulator in six different problems with
different characteristics (Sec. V-A). Overall, we are able to observe that our
method works robustly across different scenes, while the other tested methods
work well in some cases, but not in other cases (Sec. V-B). Furthermore, we
test our approach with the mobile Hubo in a real environment mimicking a
convenient store and observe similar results (Sec. V-C).

##  II RELATED WORK

In this section, we discuss prior work on planning for mobile manipulators.

###  II-A Sampling-based planning

The sampling-based approach [3, 4] is one of the most prominent strategies,
which has successfully tackled down the high-dimensional problems thanks to
its scalable space representation with the probabilistic completeness. In
addition, to get efficient movement, many manipulation approaches have used to
optimal motion planning algorithms such as Informed RRT∗ [5], BIT∗ [6], FMT∗
[7], and LazyPRM∗ [8], since Karaman et al. [9] introduced RRT∗ and PRM∗.
Burget et al. [10] suggested BI2RRT∗ that extends the Informed RRT∗ towards a
bidirectional search for a task-constrained mobile manipulator. Schmitt et al.
[11] proposed an optimal manipulation planner for continuous grasps,
placements, and actions by connecting the configurations. For achieving the
almost-sure asymptotic optimality of mobile manipulation planning, our method
is based on LazyPRM∗, while its sampling is guided by our proposed harmonious
sampling.

Guided or biased sampling with various heuristic has been widely used and
improved the performance of the optimal planners [12, 13, 14, 15]. While these
methods guided sampling for a given joint C-space, our harmonious sampling
adaptively adjusts the sampling dimensionality for the base and manipulator
configurations in mobile manipulation. Interestingly, Gochev et al. [16, 17]
suggested a dimensionality-reduction algorithm to efficiently handle high-
dimensional problems with adaptive dimensionality, which is, however,
applicable to only grid-based search algorithms with resolution-completeness.
On the other hand, our method enables the adaptive dimensionality for sampling
to the sampling-based approach for designing more effective exploration to the
mobile manipulation problem.

###  II-B Mobile manipulation planning

In order to plan an optimal trajectory for a mobile manipulation planning
problem, we have to consider the degrees-of-freedom (DoFs) of both robotic
manipulator and its base body. Unfortunately, constructing a random geometric
graph in such high-dimensional C-spaces with conventional sampling-based
planners can be a significant burden due to the curse of dimensionality.

To alleviate the complexity, a divide and conquer strategy has been widely
used [2, 18, 19]. By partitioning the entire planning problem into a set of
sub-problems (e.g., goal configuration generation and path planning separately
for base body and manipulator), it can efficiently subdivide the search space
into a set of sub-problem, while achieving a reasonable quality of the robot
trajectory.

For goal configuration generation, there have been several works using
reachability inversion [1, 20, 21]. These approaches compute a set of suitable
base positions reaching the target end-effector pose in consideration of
reachability for the manipulator. The goal configurations can be then computed
by inverse kinematics (IK) with the constructed reachability map to identify
better robot placements.

Given a target end-effector pose in the workspace, various goal configurations
can be found during the planning process by adding and considering multiple
goal configurations [22, 23]. We adopt this approach of considering multiple
goal configurations to find better ones among possible goal configurations.
Since our harmonious sampling effectively explores the joint C-space, we can
naturally consider these multiple goal configurations for the base and the
manipulator in the planning process.

##  III Overview

In this section, we give an overview of our approach. We present our adaptive
sampling algorithm to effectively explore the whole configurations of a given
mobile manipulation problem.

The main task of mobile manipulation is to reach a given target end-effector
pose, $p_{goal}$, in the task space. Due to the high dimensionality of the
mobile manipulation problem, separate planning, i.e., decoupled planning, of
the base body and manipulator is commonly adopted. Unfortunately, the
manipulator not only has high DoFs, but also is easily placed in a narrow
space, because the configurations computed from the given target pose
$p_{goal}$ are typically located very close to obstacles such as shelves and
objects. Manipulation in this narrow space with the decoupled planning may
lead to sub-optimality and increase the overall planning time due to the fixed
base chosen from the separately-run base planning for manipulation planning.
To alleviate these problems, we explore the whole C-space effectively through
our harmonious sampler guided by a low-dimensional base space.

In this paper, we assume an omnidirectional mobile manipulator robot, whose
base configuration can be represented by three parameters: 2D position and its
orientation ($q_{x},q_{y},q_{\theta}$). We define its space to be a base
space, $R$, and partition the 3D base space $R$ into a set of manipulation
regions $R_{m}$ and the rest $R_{b}(=R\setminus R_{m})$. We represent these
regions in a grid form for efficient processing. Intuitively, $R_{m}$
corresponds to difficult regions such that simultaneously adjusting the base,
$b$, and manipulator, $m$, is required. As a result, in these manipulation
regions, high-dimensional planning in the joint C-space is necessary to get a
feasible solution. On the other hand, in the base regions $R_{b}$, mainly
considering the DoFs of the base body is performed to reduce the
dimensionality in less cluttered regions. Note that our sampling process is
guided by the 3D base space $R$, but all of generated samples reside in the
joint C-space consisting of DoFs of both the base and manipulator; the samples
generated with the base regions $R_{b}$ have a predefined configuration for
joints.

Input: $p_{goal}$: target end-effector pose, $q_{pre}^{m}$: predefined
configuration of joints. $q_{init}$: initial configuration

1 $V\leftarrow\\{q_{init}\\},E\leftarrow\emptyset,Q_{goal}\leftarrow\emptyset$

2 $R\leftarrow IdentifyManipulationRegions(p_{goal})$

3 $\mathnormal{H}\leftarrow SetHarmoniousSampler(R)$

4 while  _$Termination\ condition\ is\ not\ satisfied$_ do

5 if  _$Q_{goal}=\emptyset$ or $rand(0,1)<\rho_{goal}$_ then

6 $q\leftarrow AddGoalConfiguration(p_{goal})$

7 if  _$IsCollisionFree(q)$_ then

8 Insert $q$ to $V$ and $Q_{goal}$

9

10 else

11 continue

12

13

14 else

15 $q\leftarrow HarmoniousSampling(\mathnormal{H},q_{pre}^{m})$

16 if  _$IsCollisionFree(q)$_ then

17 Insert $q$ to $V$

18

19 else

20 continue

21

22

23 $Q_{near}\leftarrow RegionSpecificNearestNeighbor(q,R)$

24 foreach  _$q_{near}\in Q_{near}$_ do

25 Insert $(q_{near},q)$ to $E$

26

27 UpdateSolutionPath($G,Q_{goal}$)

28

29

30return _$SolutionPath(G)$_

31

Algorithm 1 The proposed algorithm

Alg. 1 shows our overall approach based on LazyPRM∗ [8]. We briefly explain
how LazyPRM∗ constructs a graph $G=(V,E)$ consisting of a vertex set, $V$, and
a edge set, $E$, with our proposed methods highlighted in the yellow color. We
first identify $R_{m}$ using the reachability map and generalized Voronoi
graph (GVG), and then construct a sampling distributor to harmonize sampling
of two identified regions, $R_{m}$ and $R_{b}$ (line 2-3 of Alg. 1). We design
such a sampler, a harmonious sampler $\mathnormal{H}$, to consider the ratio
of hyper-volumes of these two different sampling spaces to draw more samples
on a larger space (line 12 of Alg. 1).

By drawing more samples in a larger space, there exists a connection issue
between samples generated differently from those two different regions. To
ameliorate this issue, we propose a region-specific k-nearest neighbor (k-NN)
search that effectively links between samples of different regions (line 17 of
Alg. 1). In addition, to select the best one among various goal
configurations, $Q_{goal}$, from a target end-effector pose $p_{goal}$, we
dynamically add goal configurations in the planning process (line 6 of Alg.
1). Lastly, we compute the shortest path from $q_{init}$ to $q_{goal}\in
Q_{goal}$ on the constructed graph $G$, and then check collisions for edges of
the path (line 20 of Alg. 1). This is because LazyPRM∗ delays edge collision
checking after finding a better path to reduce the overhead of collision
checking (line 18-19 of Alg. 1).

##  IV Main approach

(a) Reachability map.

(b) A narrow passage.

Figure 2:  These figures show our manipulation regions $R_{m}$ (visualized in
2D) computed by the target end-effector pose $p_{goal}$ and the narrow
passage. $(a)$ The blue cells indicate reachability maps computed from the
target end-effector pose (blue box). $(b)$ The green and purple spheres
represent the geometric volumes of the base $S_{b}$ and the manipulator
$S_{m}$, respectively. The red cells on the 2D floor visualize manipulation
regions caused by obstacles such that $S_{b}$ has no collision, but $S_{m}$
has collision. Black lines on the floor show the edges of the GVG.

In this section, we explain how to identify manipulation regions and how to
sample them harmoniously across different regions. We also introduce a region-
specific k-nearest neighbor (k-NN) search for effective planning with our
harmonious sampling.

###  IV-A Manipulation Region Identification

We compute the manipulation regions $R_{m}$ in the base space, to efficiently
identify difficult regions in the joint C-space where simultaneously adjusting
the base and the manipulator. To identify such regions, we adopt the
reachability map [1] from the target end-effector pose $p_{goal}$, and utilize
the GVG to deal with narrow passage problems caused by various objects. Note
that while we compute the manipulation regions in the base space, their
sampling is done in the joint C-space consisting of the base and manipulator.

Reaching the target end-effector pose. The manipulator is required to move to
reach the given target end-effector pose $p_{goal}$. For realizing the task,
we first aim to locate our robot in an appropriate base location in the base
space $R$. For computing those base locations, we use a reachability map from
the end-effector pose.

The reachability map defined in the 3D base space represents a set of suitable
base configurations computed for realizing the target end-effector pose
$p_{goal}$. The blue cells in Fig. 2(a) are 2D example cells of the
reachability map computed from the target end-effector pose shown in the
figure. When the base is located in the reachability map, it is guaranteed to
reach the target end-effector pose [1]. We set our 3D manipulation regions to
cover those reachability map. Once a random sample for the base configuration
is generated in the manipulation region, we also sample a random configuration
for the manipulator, to explore the whole high-dimensional space including the
base and manipulator configurations. Details on sample generation is in Sec.
IV-B.

Avoiding obstacles in narrow passage. When we have many obstacles, it can
potentially create narrow passages. In these narrow passages, it is required
to simultaneously adjust the base and manipulator to identify a collision-free
path. As a result, we expand the manipulation regions $R_{m}$ to include such
narrow passages as well.

To effectively identify such narrow passages, we construct GVG for the
workspace analysis. Each vertex and edge of GVG is associated with its closest
obstacles, and thus accessing those edges provides useful information on
identifying narrow passages [24]. We also adopt a _sphere expansion_ algorithm
[25], by simplifying the geometric description of a robot as a set of spheres.
Specifically, we first represent the geometric volume of the base and
manipulator by spheres, which are represented by $S_{b}$ and $S_{m}$,
respectively. We traverse edges of GVG in 2D workspace ($x$, $y$) and perform
collision detection between our sphere representations with a random
orientation and the environment. We then identify regions where $S_{b}$ has no
collision, but $S_{m}$ has collision to be difficult regions for manipulation,
and thus add them to our 3D manipulation regions $R_{m}$. Fig. 2(b) shows an
example of GVG in the tested environment and our sphere representations of the
base and the manipulator of the predefined pose.

Input: $\mathnormal{H}$: harmonious sampler, $q_{pre}^{m}$: predefined
configuration of joints.

1 if  _$rand(0,1) <\rho_{sample}$_ then

2 $q_{rand}^{b}\leftarrow GetRandomSampleForBase(\mathnormal{H})$

3 if  _$InManipulationRegions(q_{rand}^{b})$_ then

4 $q_{rand}^{m}\leftarrow GetRandomSampleForManipulator()$

5

6 else

7 $q_{rand}^{m}\leftarrow q_{pre}^{m}$

8

9

10else

11 $q_{rand}\leftarrow UniformSampler()$

12

13return _$q_{rand}$_

14

Algorithm 2 Harmonious sampling

###  IV-B Harmonious Sampling

The purpose of our harmonious sampler $\mathnormal{H}$ is to balance the
sampling distribution between two different regions, i.e., $R_{m}$ and
$R_{b}$. To achieve our goal, we design $\mathnormal{H}$ to adjust the
sampling probability adaptively by considering the ratio of hyper-volumes of
these two different sampling spaces to draw more samples on a larger space.
Note that the manipulation regions $R_{m}$ is defined in the base space $R$,
but is associated with the larger sampling space consisting of the base and
manipulator, i.e., the joint high-dimensional C-space.

Our harmonious sampler $\mathnormal{H}$ first generates a random sample for
the base configuration $q_{rand}^{b}$ (line 2 of Alg. 2) and then we check
whether the sample is in either $R_{m}$ or $R_{b}$. If $q_{rand}^{b}\in
R_{m}$, we also sample a random configuration of the manipulator,
$q_{rand}^{m}$ (line 3-4 of Alg. 2). Otherwise, $q_{rand}^{m}$ is initialized
with predefined values of joints, $q_{pre}^{m}$, which can be initial values
of joints (line 5-6 of Alg. 2). Note that $q_{rand}$ is composed of
$q_{rand}^{b}$ and $q_{rand}^{m}$(=$q_{rand}\setminus q_{rand}^{b}$).

To generate a random sample for the base configuration $q_{rand}^{b}$
according to our harmonious sampler $\mathnormal{H}$, we use a probability
mass function (PMF) defined in a grid form over the 3D base space. Each cell
of the grid over the 3D base space has only two types of indicating whether
the cell is in $R_{m}$ or $R_{b}$. Our harmonious sampler $\mathnormal{H}$
then uses the inverse transform sampling method that works for an arbitrary
sampling PMF.

To compute a sampling probability for two regions $R_{m}$ and $R_{b}$, we
utilize the notion of the sampling hyper-volume associated to each cell,
$\gamma_{c}$. Specifically, when a cell $\gamma_{c}$ is on the manipulation
region $R_{m}$, its sampling probability is computed based on the hyper-volume
of its sampling space, $\nu(\gamma_{c})$, which is simply calculated by
multiplying the intervals between the lower bound, $q_{d}^{min}$, and upper
bound, $q_{d}^{max}$, for each parameter space $d$ of its sampling space:

$\nu(\gamma_{c}|\gamma_{c}\in
R_{m})=\prod_{d=1}^{D(R_{m})}w_{d}*(q_{d}^{max}-q_{d}^{min}),$ (1)

where $D(R_{m})$ is the DoF associated with the sampling space of $R_{m}$, and
$w_{d}$ is the weight to match the units for each dimension (Sec. IV-D). This
approach is also used for computing a sampling probability for cells on the
base regions $R_{b}$.

The PMF for our harmonious sampler $\mathnormal{H}$ of cells $\gamma_{c}$ in
the base space $R$ can be defined as:

$\mathbb{P}(\gamma_{c}|R)=\nu(\gamma_{c})/\sum\limits_{i=1}^{N}\nu(\gamma_{i}),$
(2)

where $N$ is the number of cells in the base space $R$, whose cell index is
$i$. Consequently, the harmonious sampler $\mathnormal{H}$ varies the sampling
probability according to sampling DoFs associated with each cell, i.e. three
sampling DoFs for $R_{b}$ and ten sampling DoFs for $R_{m}$.

Note that our harmonious sampler $\mathnormal{H}$ is simple yet useful, but
cannot guarantee the optimality. To ensure the probabilistic completeness and
the asymptotic optimality of our harmonious sampling, we use both our
harmonious sampler $\mathnormal{H}$ at a probability of $\rho_{sample}$ and
the uniform sampler over the entire sampling space at a probability of
$(1-\rho_{sample})$ [13, 14] (Alg. 2).

###  IV-C Linking Different Regions

Our harmonious sampling adaptively adjusts its sampling space depending on
whether a sample is generated on the base regions $R_{b}$ or manipulation
regions $R_{m}$. Note that we generate more samples on the joint C-space
associated with $R_{m}$, compared to that associated with $R_{b}$, thanks to
the wider sampling DoF associated with $R_{m}$. Fig. 3 illustrates the
different sampling densities on the joint C-space; in this example, $R_{m}$
and $R_{b}$ are illustrated in the bottom one-dimensional space.

Figure 3:  Our region-specific k-NN links samples generated through different
regions. Blue and orange dotted lines, i.e., connections, are made by our
approach for the red dot.

We found that this varying sampling probability causes a sub-optimal
performance, especially due to sub-optimal linking along the boundary of two
different regions. For example, when we naively perform k-NN search for
samples (e.g., the red dot in Fig. 3) generated from the manipulation regions
$R_{m}$, identified samples are likely to be ones from $R_{m}$ due to its high
sample probability causing small distances between samples generated from
$R_{m}$ (Fig 3).

To alleviate this connectivity issue between two different regions, we propose
to use a simple region-specific k-NN search. For samples associated with the
manipulation regions, we perform two independent k-NN search: one is performed
only with samples associated with $R_{m}$, resulting in blue dotted lines in
Fig. 3, and the other k-NN is performed only with samples associated with
$R_{b}$, resulting in the orange dotted lines in Fig. 3, for improving the
connectivity between two different regions. For samples associated with the
base regions, we simply use the original k-NN search that is performed with
all the available samples.

We found that our region-specific k-NN search may generate unnecessary edges
between samples that are located in different regions, causing a more memory
overhead, approximately 25% over using the original k-NN search. Fortunately,
our method is based on LazyPRM∗ [8], and thus there is not much difference in
terms of the number of edge collision checking, since we perform edge
collision checking on the solution paths. As a result, using region-specific
k-NN search improves the path costs by 7% in our tested environments without a
significant runtime performance overhead, thanks to better connectivity
between two different regions (Table V-B).

###  IV-D Distance Metric

We perform planning in the joint C-space consisting of the base and
manipulator. Since the base and joints have different kinds of quantities,
i.e., linear and angular quantities, we have the incompatibility issue, which
can lower the efficacy of performing nearest neighbor search during planning
with respect to the asymptotic optimality [9].

In case of the incompatibility, LaValle et al. [26] considered a robot
displacement by using coefficients of different quantities as weights, $w$, to
match the units, as the following:

$dist=\sqrt{\sum_{d}^{D}(w_{d}*\rVert q_{d}-q^{\prime}_{d}\rVert)^{2}},$ (3)

where $D$ is the DoF of the joint C-space; for the mobile Hubo robot,
$q_{x},q_{y},q_{\theta}$ correspond to the base configuration, and the rest to
the manipulator.

We explain how to determine weights $w_{i}$ to match the units using the robot
displacement. For two configurations $q$ and $q^{\prime}$, the robot
displacement metric, $disp(\cdot,\cdot)$, can be defined as:

$disp(q,q^{\prime})=\max_{a\in\mathnormal{A}}{\rVert\tau(q,a)-\tau(q^{\prime},a)\rVert},$
(4)

where $\tau(q,a)$ is the position of point $a$ for the robot $\mathnormal{A}$
with the configuration $q$ in the workspace.

Since the robot displacement metric $disp(\cdot,\cdot)$ yields the maximum
amount in the workspace, the first joint with the largest change of the
manipulator should have a larger weight. While fixing weights of $q_{x},q_{y}$
to $1$ for the simplicity, the maximum displacement for each angular quantity
is treated as arc lengths in the workspace that can be made given its joint.
One can then show that the weight for a joint of a manipulator can be computed
by the distance from the joint to the end-effector; the weight intuitively
corresponds to the maximum arc length per radian.

#### Supporting multi-goal configurations.

Given the target end-effector pose $p_{goal}$, our work supports multi-goal
configurations in the joint C-space by planning the base and manipulator
simultaneously with our harmonious sampling. Since possible goal
configurations from $p_{goal}$ can be infinite, our method dynamically adds
goal configurations by injecting the reachability map and IK solver [23]
during planning according to a probability of $\rho_{goal}$ (line 5-6 of Alg.
1). Our harmonious sampler $\mathnormal{H}$ then generate samples in the joint
C-space according to $R_{m}$, naturally supporting dynamically added goal
configurations.

##  V Experiments and Analysis

(a) Grasping an object on the table (Problem 1).

(b) Grasping an object on the table with ground obstacles (Problem 2).

(c) Grasping an object on the shelf (Problem 3).

(d) Grasping an object on the shelf with ground obstacles (Problem 4).

(e) The worst case of the predefined posture (Problem 5).

(f) The narrow passage (Problem 6).

Figure 4:  This shows our test scenes with the start and goal configurations,
$q_{init}$ and $q_{goal}$. $(a)$ and $(b)$ are to grasp an object on the
table. $(c)$ and $(d)$ are to grasp an object on the shelf. $(b)$ and $(d)$
consist of additional objects, obstructing the robot. $(e)$ shows a difficult
case of the predefined posture in $(b)$. The manipulator is located under the
table. $(f)$ shows the narrow passage where adjusting a manipulator is needed
to move a base.

We test our method on a machine that has 3.40GHz Intel i7-6700 CPU and 16GB
RAM. We use the simulated mobile Hubo robots for testing our method within the
OpenRAVE simulator (Fig. 4). We also use a real mobile Hubo robot [27] for
testing the usability in a real environment (Fig. 1). A manipulator of the
robot has seven joints, and its base consists of $x$, $y$ and $\theta$. We
partition the base space $R$ into a grid form using $0.1m$ resolution for $x$,
$y$ and $30$ degree resolution for $\theta$. We use LazyPRM∗ [8] as our base
optimal planner, and set $\rho_{sample}$ $0.8$ and $\rho_{goal}$ $0.01$ for
our tests.

To validate our method, we compare our method against two different
approaches: the decoupled and coupled approaches. The decoupled approach [2]
divides the planning process into two separate planning for the base and
manipulator, while the coupled approach [10] performs the planning within the
joint C-space representing both base and manipulators and thus having 10 DoFs.

For the decoupled planning, we need to allocate separate time budget for base
and manipulator planning parts for the same-time comparison against other
approaches. Unfortunately, this issue was not explicitly discussed in prior
techniques, and left for user setting. We test different time budgets for the
decoupled method, and found $1.5:8.5$ time budget ratio for the base and the
manipulator shows the best performance on average. Fig. 5(c) shows results
with varying ratios in one of difficult testing environments.

We experiment for various planning with a single goal configuration except for
the multiple goal approach generating multiple goal configurations given the
target end-effector pose. For all the experiments, we perform $50$ tests given
a time budget of 100 seconds, and report their average performance.

(a) Result of Problem 3 w/ different methods

(b) Result of Problem 4 w/ different methods

(c) Results w/ different time budget ratios for the decoupled approach in
Problem 4

Figure 5:  This figure shows the cost over planning time of different methods.
$(a)$ and $(b)$ show that our method (Ours + Multi) achieves the best
performance. $(c)$ shows that $1.5:8.5$ time budget ratio for the base and the
manipulator used in the decoupled approach shows the best performance among
the tested budget ratios.We visualize graphs from the time that an initial
solution is computed.

###  V-A Test problems

We prepare six different problems to validate our approach, as shown in Fig.
4. Our test scene is $5m\times 4m$ wide. The first and second problems (table
scenes) are designed to grasp an object on the table as relatively simpler
problems with a few obstacles for the manipulator. The third and fourth
problems (shelf scenes) have a narrow shelf with obstacles on the bookshelf.
Moreover, to show the effect of base movement, the second and fourth problems
additionally place ground obstacles that interfere with the movement of a
robot.

To intentionally create a difficult case that we may encounter in practice, we
use predefined postures for the fifth and sixth that are difficult for
manipulating objects or passing narrow passages without dynamically modifying
the predefined postures. For other scenes, we use the basic posture of the
mobile Hubo robot.

TABLE I:  Various statistics, averaged out from 50 independent tests. For the decoupled approach, we also show separate times (s) spent on finding an initial solution for each of the base and the manipulator.  | Decoupled | Coupled | Ours + | Ours +
---|---|---|---|---
| Approach | Approach | Single | Multi
| # of time-outs | 0 | 0 | 0 | 0
Problem 1 | Initial solution time | 5.70 | 11.52 | 2.23 | 2.45
Fig. 4(a) | (Base/Manipulator) | (0.57/5.13)
| Final distance cost | 4.61 | 8.10 | 4.61 | 4.03
| # of time-outs | 0 | 0 | 0 | 0
Problem 2 | Initial solution time | 6.55 | 27.72 | 2.79 | 3.08
Fig. 4(b) | (Base/Manipulator) | (0.93/5.62)
| Final distance cost | 5.40 | 9.65 | 5.18 | 4.81
| # of time-outs | 11 | 0 | 0 | 0
Problem 3 | Initial solution time | 39.15 | 9.25 | 7.00 | 4.14
Fig. 4(c) | (Base/Manipulator) | (0.75/38.40)
| Final distance cost | 6.92 | 8.73 | 5.71 | 5.62
| # of time-outs | 12 | 21 | 0 | 0
Problem 4 | Initial solution time | 41.50 | 52.74 | 24.34 | 30.26
Fig. 4(d) | (Base/Manipulator) | (2.05/39.45)
| Final distance cost | 7.24 | 12.76 | 6.37 | 6.23
| # of time-outs | 0 | 0 | 0 | 0
Problem 5 | Initial solution time | 9.55 | 31.96 | 3.94 | 3.65
Fig. 4(e) | (Base/Manipulator) | (2.07/7.48)
| Final distance cost | 6.97 | 9.84 | 5.79 | 5.08
| # of time-outs | - | 28 | 0 | 0
Problem 6 | Initial solution time | - | 70.90 | 36.91 | 33.13
Fig. 4(f) | (Base/Manipulator)
| Final distance cost | - | 16.45 | 12.16 | 11.73

###  V-B Result Analysis

We compare our method with the coupled and decoupled approaches across six
different problems. Initially, our method identifies the manipulation regions
$R_{m}$ and constructs the harmonious sampler $\mathnormal{H}$ ($\approx$0.5s
at the most); these times are included in the reported planning time. Table I
shows various run-time results of the tested methods. Problem 1 and 2 are
simple with a few obstacles, so all the methods find a solution in reasonable
planning time (Table I). In these simple cases, the simplest, decoupled
approach finds initial solutions faster than the coupled approach that
performs sampling in the joint C-space. For Problem 3, on the other hand, the
coupled approach finds solutions faster than the decoupled approach.
Furthermore, the decoupled approach did not find solutions in $11$ of $50$
tests. This clearly indicates the problem of the decoupled approach that does
not work well at the narrow space. Overall, the decoupled and coupled
approaches have the weakness in terms of showing the robust performance across
different scenes.

Fortunately, our methods, especially, Ours + Multi considering multiple goal
configurations, find solutions faster than both coupled and decoupled
approaches, and achieve the shortest path in all of Problem 1 to 4. Fig. 5(a)
and Fig. 5(b) show how different methods behave as a function of planning
time. Our method (Ours + Multi) computes shorter paths given any fixed
planning time. These results clearly show benefits of our approach that
adaptively performs sampling for base and manipulators in their joint C-space,
even in simple and middle-level problems.

Let us now consider most difficult problems among the tested ones. In Problem
5, our method finds the initial solution much faster, more than one order of
magnitude for Ours + Multi, than the coupled and decoupled approaches. In the
case of Problem 6, also, ours achieves better solutions than the coupled
approach, even though the decoupled cannot solve the problem. These results
demonstrate the efficiency and robustness of our methods.

Single vs. multiple goal configurations. Considering multiple goal
configurations shows to produce smaller costs across all the tests given the
time budget (100s), thanks to considering various candidate configurations to
reach the target end-effector pose (Table I). Nonetheless, ours with multiple
goal configurations shows similar performance in terms of finding an initial
solution over considering only a single goal configuration. This is because
the overhead of considering multiple goal configurations is about 2.8% of the
overall planning time, resulting in slower performance for identifying the
initial solution. However, this overhead is small, and thus is paid well as we
have more planning time for achieving shorter paths.

TABLE II:  The results of the simple problem 2 and the difficult problem 4
with varying volume ratios for our harmonious sampler and region-specific NN
search.

k-NN type | Original | Region-specific
---|---|---
Ratio of $\nu(\gamma_{c}|\gamma_{c}\in R_{m})$ | Ours | 0.2$\times$Ours | Ours (Eq. 1) | 5$\times$Ours
Problem 2  Fig. 4(b) | # of time-outs | 0 | 0 | 0 | 0
Initial solution time | 2.99 | 5.26 | 2.79 | 3.08
Final distance cost | 5.18 | 5.45 | 5.18 | 5.33
Problem 4  Fig. 4(d) | # of time-outs | 4 | 1 | 0 | 1
Initial solution time | 31.93 | 37.19 | 24.34 | 29.61
Final distance cost | 6.85 | 6.79 | 6.37 | 6.61

Analysis of the harmonious sampler. To see benefits of our harmonious sampler
$\mathnormal{H}$, we check how the performance of our method behaves according
to $\nu(\gamma_{c}|\gamma_{c}\in R_{m})$ of Eq. 1 for Problem 2 and 4, given
the region-specific NN search. Specifically, we compare the performance of our
approach (Ours) against decreasing and increasing
$\nu(\gamma_{c}|\gamma_{c}\in R_{m})$ by five times; they are denoted as
0.2$\times$Ours and 5$\times$Ours, respectively. Table V-B shows that the
number of time-outs, cost, and planning time are increased when
$\nu(\gamma_{c}|\gamma_{c}\in R_{m})$ is increased or decreased. These results
support the benefit of our harmonious sampler $\mathnormal{H}$ taking into
account the hyper-volumes of different regions.

Analysis of the region-specific k-NN search. To see the benefits of our
region-specific k-NN search, we compare its performance against the original
k-NN search that performs k-NN search with all the samples in the joint
C-space. Table V-B shows that in the simple problem (Problem 2), there is not
much difference in terms of the final distance cost between two methods. On
the other hand, in a difficult problem (Problem 4), our region-specific k-NN
search (Ours) reduces all of three measures (including the final distance
cost) over using the original k-NN search. Specifically, using the region-
specific k-NN search shows $1.3$ times performance improvement in terms of
finding an initial solution time, and $7\%$ cost reduction, compared to using
the original k-NN search.

TABLE III:  Results with a predefined posture, which is chosen from five random postures. These results are averaged out from 50 independent tests.  5 Random predefined postures | Decoupled | Coupled | Ours + | Ours +
---|---|---|---|---
Approach | Approach | Single | Multi
| # of time-outs | 0 | 0 | 0 | 0
Problem 1 | Initial solution time | 4.79 | 8.88 | 2.91 | 4.18
Fig. 4(a) | (Base/Manipulator) | (0.54/4.25)
| Final distance cost | 6.75 | 9.04 | 6.76 | 5.95
| # of time-outs | 13 | 21 | 1 | 1
Problem 4 | Initial solution time | 41.78 | 55.19 | 17.04 | 12.64
Fig. 4(d) | (Base/Manipulator) | (1.60/40.18)
| Final distance cost | 9.28 | 12.79 | 8.27 | 7.72
| # of time-outs | - | 29 | 0 | 0
Problem 6 | Initial solution time | - | 63.39 | 42.73 | 45.57
Fig. 4(f) | (Base/Manipulator)
| Final distance cost | - | 16.13 | 14.43 | 14.65

Analysis with different predefined postures. We also test different methods
with five random postures as the predefined posture. Table III shows results
in Problem 1, 4, and 6; other scenes show similar tendencies. Our methods find
solutions faster than both decoupled and coupled approaches and achieve the
shortest paths in all of tested problems even with different initial postures,
demonstrating the robustness of our approach. These results show that our
method robustly achieves a better performance than both decoupled and coupled
approaches.

Figure 6:  This figure shows our experimental environment mimicking a
convenient store with the mobile Hubo robot. There are various kinds of
beverages and snacks on the shelf, and our goal is to grasp the yellow
beverage. Our approach computes a shorter path over the decoupled approach in
the same time budget.

###  V-C Real Robot Test

We integrated our approach with the real, mobile Hubo robot. We tested in an
experimental environment designed for mimicking a convenience store with
various kinds of beverages and snacks on the shelf (Fig. 6). We defined a
scenario that the robot moves to the shelf and grasps the yellow beverage on
the shelf (Fig. 1). In this scenario, we compared our method with the
decoupled approach. Our method stretched out the manipulator while moving the
base body. On the other hand, the decoupled approach moves to the shelf
without adjusting the manipulator, and then adjusts the manipulator to grasp
the yellow beverage, resulting in an inferior path. As a result, our method
grasped the target object with a better solution than the decoupled approach;
i.e., the cost of our approach is 3.962, while the decoupled one has 4.318.

##  VI Conclusion and future work

In this paper, we have introduced a simple, yet effective mobile manipulation
sampling method named harmonious sampling that adaptively adjusts the sampling
space for the base and the manipulator. Our harmonious sampling works upon the
low-dimensional base space to rapidly explore a wide-open space, while
exploring in a high-dimensional space where it is necessary to consider full
DoFs of the mobile manipulator. Through the experiments, we have shown the
improvement of performance in a variety of indoor environments with practical
tasks. Furthermore, we have tested our approach by integrating with the mobile
Hubo in the real environment.

In general, it is critical to consider the uncertainty of sensing and robot
control. Since our method simultaneously operates for the base and
manipulator, it is important to handle the uncertainties, which are under
active research in robotics. In future work, we would like to efficiently
handle them, and support anytime planning.

## Acknowledgment

We would like to thank anonymous reviewers for constructive comments. This
work was supported by MI/KEIT (10070171) and NRF/MSIT (No.2019R1A2C3002833).

## References

  * [1] Nikolaus Vahrenkamp, Tamim Asfour, and Rüdiger Dillmann,  “Robot placement based on reachability inversion”,  in ICRA. IEEE, 2013, pp. 1970–1975.
  * [2] Sachin Chitta, E Gil Jones, Matei Ciocarlie, and Kaijen Hsiao,  “Perception, planning, and execution for mobile manipulation in unstructured environments”,  Robotics and Automation Magazine, vol. 19, no. 2, pp. 58–71, 2012\.
  * [3] Lydia E Kavraki, Petr Svestka, J-C Latombe, and Mark H Overmars,  “Probabilistic roadmaps for path planning in high-dimensional configuration spaces”,  IEEE Transactions on Robotics and Automation, vol. 12, no. 4, pp. 566–580, 1996.
  * [4] Steven M LaValle,  “Rapidly-exploring random trees a new tool for path planning”,  Tech. Rep. 98-11, Iowa State University, 1998.
  * [5] Jonathan D Gammell, Siddhartha S Srinivasa, and Timothy D Barfoot,  “Informed RRT*: Optimal sampling-based path planning focused via direct sampling of an admissible ellipsoidal heuristic”,  in IROS. IEEE, 2014, pp. 2997–3004.
  * [6] Jonathan D Gammell, Siddhartha S Srinivasa, and Timothy D Barfoot,  “Batch informed trees (BIT*): Sampling-based optimal planning via the heuristically guided search of implicit random geometric graphs”,  in ICRA. IEEE, 2015, pp. 3067–3074.
  * [7] Lucas Janson, Edward Schmerling, Ashley Clark, and Marco Pavone,  “Fast Marching Tree: A fast marching sampling-based method for optimal motion planning in many dimensions”,  The International Journal of Robotics Research, vol. 34, no. 7, pp. 883–921, 2015.
  * [8] Kris Hauser,  “Lazy collision checking in asymptotically-optimal motion planning.”,  in ICRA. IEEE, 2015, pp. 2951–2957.
  * [9] Sertac Karaman and Emilio Frazzoli,  “Sampling-based algorithms for optimal motion planning”,  The International Journal of Robotics Research, vol. 30, no. 7, pp. 846–894, 2011.
  * [10] Felix Burget, Maren Bennewitz, and Wolfram Burgard,  “BI2RRT*: An efficient sampling-based path planning framework for task-constrained mobile manipulation”,  in IROS. IEEE, 2016, pp. 3714–3721.
  * [11] Philipp S Schmitt, Werner Neubauer, Wendelin Feiten, Kai M Wurm, Georg V Wichert, and Wolfram Burgard,  “Optimal, sampling-based manipulation planning”,  in ICRA. IEEE, 2017, pp. 3426–3432.
  * [12] Brendan Burns and Oliver Brock,  “Toward optimal configuration space sampling.”,  in RSS, 2005, pp. 105–112.
  * [13] Baris Akgun and Mike Stilman,  “Sampling heuristics for optimal motion planning in high dimensions”,  in IROS. IEEE, 2011, pp. 2640–2645.
  * [14] Jauwairia Nasir, Fahad Islam, Usman Malik, Yasar Ayaz, Osman Hasan, Mushtaq Khan, and Mannan Saeed Muhammad,  “RRT*-SMART: A rapid convergence implementation of rrt”,  International Journal of Advanced Robotic Systems, vol. 10, no. 7, pp. 299, 2013.
  * [15] Scott Kiesel, Tianyi Gu, and Wheeler Ruml,  “An effort bias for sampling-based motion planning”,  in IROS. IEEE, 2017, pp. 2864–2871.
  * [16] Kalin Gochev, Alla Safonova, and Maxim Likhachev,  “Planning with adaptive dimensionality for mobile manipulation”,  in ICRA. IEEE, 2012, pp. 2944–2951.
  * [17] Kalin Gochev, Alla Safonova, and Maxim Likhachev,  “Incremental planning with adaptive dimensionality.”,  in ICAPS, 2013.
  * [18] Siddharth Srivastava, Eugene Fang, Lorenzo Riano, Rohan Chitnis, Stuart Russell, and Pieter Abbeel,  “Combined task and motion planning through an extensible planner-independent interface layer”,  in ICRA. IEEE, 2014, pp. 639–646.
  * [19] Caelan Reed Garrett, Tomás Lozano-Pérez, and Leslie Pack Kaelbling,  “Sample-based methods for factored task and motion planning”,  in RSS, 2017.
  * [20] Jun Dong and Jeffrey C Trinkle,  “Orientation-based reachability map for robot base placement”,  in IROS. IEEE, 2015, pp. 1488–1493.
  * [21] Andreas Hertle and Bernhard Nebel,  “Identifying good poses when doing your household chores: Creation and exploitation of inverse surface reachability maps”,  in IROS. IEEE, 2017, pp. 6053–6058.
  * [22] Mike Stilman, Jan-Ullrich Schamburek, James Kuffner, and Tamin Asfour,  “Manipulation planning among movable obstacles”,  in ICRA. IEEE, 2007, pp. 3327–3332.
  * [23] Dmitry Berenson, Siddhartha S Srinivasa, Dave Ferguson, Alvaro Collet, and James J Kuffner,  “Manipulation planning with workspace goal regions”,  in ICRA. IEEE, 2009, pp. 618–624.
  * [24] Donghyuk Kim, Junghwan Lee, and Sung-eui Yoon,  “Cloud RRT∗: Sampling cloud based rrt∗”,  in ICRA. IEEE, 2014, pp. 2519–2526.
  * [25] Yuandong Yang and Oliver Brock,  “Efficient motion planning based on disassembly”,  in RSS, 2005, pp. 97–104.
  * [26] Steven M LaValle,  Planning algorithms,  Cambridge university press, 2006.
  * [27] M. Lee, Y. Heo, J. Park, H. Yang, P. Benz, H. Jang, H. Park, I. Kweon, and J. Oh,  “Fast perception, planning, and execution for a robotic butler: Wheeled humanoid m-hubo”,  in IROS. IEEE, 2019.