# Supplementary material for ’Interplay of superexchange and vibronic effects
in the hidden order of Ba2MgReO6 unravelled from first principles’

Dario Fiore Mosca  CPHT, CNRS, École polytechnique, Institut Polytechnique de
Paris, 91120 Palaiseau, France  Collège de France, Université PSL, 11 place
Marcelin Berthelot, 75005 Paris, France  Cesare Franchini  University of
Vienna, Faculty of Physics and Center for Computational Materials Science,
Vienna, Austria  Department of Physics and Astronomy ”Augusto Righi”, Alma
Mater Studiorum - Università di Bologna, Bologna, 40127 Italy  Leonid V.
Pourovskii  CPHT, CNRS, École polytechnique, Institut Polytechnique de Paris,
91120 Palaiseau, France  Collège de France, Université PSL, 11 place Marcelin
Berthelot, 75005 Paris, France

##  I First principles methods

In the following, we first describe our electronic structure calculations for
paramagnetic BMRO. We then proceed with the description of how the intersite
exchange interaction (IEI) and Jahn-Teller (JT) terms of the full effective
Hamiltonian (5) are calculated on the basis of this electronic structure.
Lastly, we describe the intersite elastic coupling (EC) calculations by the
Density Functional Perturbation Theory (DFPT) approach.

###  I.1 Correlated electronic structure calculations

The electronic structure of Ba2MgReO6 in the paramagnetic phase is computed
with the DFT+dynamical mean-field theory(DMFT) method. The quantum impurity
problem for the $t_{2g}$ shell of the Re ion is solved within the quasi-atomic
Hubbard-I (HI) approximation Hubbard (1963); the method is abbreviated below
as DFT+HI. We employ a full charge self-consistent DFT+DMFT
implementationAichhorn  _et al._ (2009, 2011, 2016) based on the full-
potential LAPW code Wien2kBlaha  _et al._ (2018) that includes the spin-orbit
coupling effect with the standard variational treatment.

The Wannier orbitals representing Re $t_{2g}$ orbitals are constructed from
the manifold of $t_{2g}$ Kohn-Sham (KS) bands enclosed by the energy window
[-1.36:2.04] eV relative to the KS Fermi level. The Kanamori rotationally-
invariant on-site Coulomb repulsion between Re $t_{2g}$ is parameterized by
$U_{\mathrm{K}}=3.77$ eV and $J_{\mathrm{H}}^{\mathrm{K}}=0.39$ eV, which
correspond to the full $d$-shell parameters $F^{0}=U=$3.2 eV and
$J_{\mathrm{H}}$=0.5 eV. Those values of $U$ and $J_{\mathrm{H}}$ are in
agreement with our previous works on $d^{1}$ and $d^{2}$ DP Fiore Mosca  _et
al._ (2021); Pourovskii  _et al._ (2021).

The calculations of Figure 1 in the main are carried out for the experimental
cubic lattice structures of Ba2MgReO6 with lattice parameter $a=$ 8.0802 Å .
The local density approximation of the DFT exchange-correlation potential is
employed, together with 400 ${\bf k}$-point mesh in the full Brillouin zone,
and the Wien2k basis cutoff $R_{\mathrm{mt}}K_{\mathrm{max}}=$7\. The double-
counting correction is evaluated using the fully-localized limit with the
nominal 5$d$ shell occupancy of 1.

###  I.2 Calculation HIEI

We evaluate the IEI from the converged electronic paramagnetic electronic
structure of BMRO using the force-theorem in Hubbard-I (FT-HI)
approachPourovskii (2016). This method considers small symmetry-breaking
fluctuations in the density matrix of the ground state (GS)
$J_{\mathrm{eff}}=3/2$ multiplet that occur simultaneously on two neighboring
magnetic (Re) sites $i$ and $j$. The IEI $V_{KK^{\prime}}^{QQ^{\prime}}(ij)$
is then evaluated by considering the response of DFT+DMFT grand potential to
these two-site fluctuations. The method is similar in spirit to the force-
theorem methods for the symmetry-broken magnetic state of Refs. Liechtenstein
_et al._ , 1987; Katsnelson and Lichtenstein, 2000, but is formulated,
instead, for the symmetry-unbroken paramagnetic state. The application of FT-
HI in the case of BMRO is fully analogous to its previous applications to
$J_{\mathrm{eff}}=3/2$ DP Fiore Mosca  _et al._ (2021); Pourovskii (2023). The
reader can thus find more detailed description in Appendix of Ref. Pourovskii,
2023 and Suppelementary of Ref. Fiore Mosca  _et al._ , 2021; the full
derivation of the method is given in Ref. Pourovskii, 2016.

###  I.3 Calculation $H_{JT}$

Since the distortions of ReO6 octahedra couple to quadrupoles and we are
interested in the competition between the F$xy$ and AF$x^{2}-y^{2}$
quadrupolar orders, we consider the corresponding oxygen sublattice
distortions belonging to the Brillouin zone $\Gamma$ and $X$-points
irreducible representations that induce $T_{2g}$ and $E_{g}$ local distortions
of the ReO6 octahedra. The corresponding $\Gamma_{5}^{+}$ and $X_{2}^{+}$
structures (see Suppl. Fig. 1) were generated with the help of ISODISTORT tool
Campbell  _et al._ (2006); Stokes  _et al._ (2022) for a set of small
distortions with the amplitude in the range from about 0.01 Å to 0.03 Å.

Figure 1: Distorted unit cells of BMRO implementing the $\Gamma_{5}^{+}$ and
$X_{2}^{+}$ distortions. The displayed distortion amplitude is strongly
exaggerated for better visibility as compared to the distortions amplitude
actually used in calculations. The plots are generated by VESTA Momma and
Izumi (2011).

Subsequently we carry out DFT+HI calculations using the same setup as
described in Sec. I.1. In these self-consistent calculations, we follow the
approach of Ref. Delange  _et al._ , 2017 by averaging the Boltzmann weights
within the $J_{\mathrm{eff}}$=3/2 quadruplet in order to suppress the
unphysical self-interaction contributionBrooks  _et al._ (1997) to the
quadruplet’s splitting.

Due to the JT coupling $-g_{\Gamma}O^{Q}_{2}q_{\Gamma Q}$, octahedral
distortions split the $J_{\mathrm{eff}}$=3/2 quadruplet into two Kramers
doublets. The eigenvalues of normalized $J_{\mathrm{eff}}$=3/2 quadrupolar
operators $O_{2}^{Q}$ are $\pm 1/2$. Hence, the JT splitting of the quadruplet
is simply given by $g_{\Gamma}q_{\Gamma Q}$. As shown in Suppl. Fig. 2a, the
calculated splitting indeed shows a perfect linear scaling vs. the distortion
amplitude $q_{\Gamma Q}$ within the studied range, hence, no quadratic JT
couplingIwahara and Chibotaru (2023) needs to be introduced. The larger slope
in the case of $X_{2}^{+}$ distortion is due to a stronger JT coupling to the
$E_{g}$ octahedral distortions.

Figure 2: a. Calculated $J_{\mathrm{eff}}$=3/2 quadruplet splitting vs.
distortion amplitude for BMRO at ambient pressure is shown by symbols. The
lines are linear fits to the data. b. Calculated JT coupling vs. pressure
(symbols). The lines are linear (for $X_{2}^{+}$) and quadratic (for
$\Gamma_{5}^{+}$) regression fits to the data.

We repeated the same calculations for a set of volumes corresponding to
compressed and expanded BMRO. The resulting JT coupling are depicted in Suppl.
Fig. 2b vs. pressure (calculated from volume using our theoretical equation of
states, see Sec. III); the JT coupling shows a smooth increase vs. pressure
allowing to fit it with a simple linear or quadratic regression. As one sees,
the $E_{g}$ JT coupling is almost 3 times stronger than the $T_{2g}$ one. Our
values for the JT coupling calculated by DFT+HI agree well with the DFT
results of Ref. Soh  _et al._ , 2023.

###  I.4 Calculation $H_{EC}$

To determine the effective Hamiltonian $H_{EC}$, we utilize the projective
phonon approach explained in detail in the subsequent section. For the phonon
calculations in Ba2MgReO6 we adopted the DFPT approach as incorporated in the
Vienna Ab initio Simulation Package (VASP) Kresse and Hafner (1993); Kresse
and Furthmüller (1996). These calculations involve the use of DFT with LDA
exchange-correlation potential, employing a $2\times 2\times 2$ k-mesh and a
high energy cutoff of 600 eV together with a total energy convergence of 10-6
eV. The choice of an LDA potential has been based on difficulties in
converging accurate phonons within the Perdew-Burke-Ernzerhof exchange
correlation functional. We note, however, that using LDA in this context has
proven to be accurate in determining both volumetric, bond-length and pressure
properties (See Suppl. Table 4 and Table 5). The calculations are based on a
$2\times 2\times 2$ supercell structure containing 32 f.u. This supercell was
generated through the phonopy package Togo  _et al._ (2023); Togo (2023)
starting from a fully relaxed structure with lattice parameter of 7.9992 Å.

For the structural relaxation of the smaller 8 f.u. supercell we used the LDA
potential together with a $12\times 12\times 12$ k-mesh, an energy cutoff of
500 eV, and a convergence criterion of 10-5 eV on the electronic self-
consistent cycle. The volume minimization and atomic positions were optimized
using a quasi-Newton algorithm with a step width of 0.2 Å and a convergence of
ionic relaxation of 10-3 eV. Detailed information about the optimized
structure is available in the cif file, which is included at the conclusion of
the Supplementary Material.

##  II Projective Phonon Formalism

The convention used in phonopy for the dynamical matrix calculation is the
following:

$D_{\alpha\beta}(ij,\mathbf{q})=\frac{1}{\sqrt{m_{i}m_{j}}}\sum_{l}\Phi_{\alpha\beta}(i0,jl)exp(i\mathbf{q}\cdot[\mathbf{r}(jl)-\mathbf{r}(i0)])\
,$ (1)

where

$\Phi_{\alpha\beta}(i0,jl)=\frac{\partial^{2}V}{\partial r_{\alpha
0}(i)r_{\beta l}(j)}$ (2)

is the corresponding force constant matrix. Here ($\alpha$, $\beta$) represent
the Cartesian coordinates, ($i$,$j$) the atomic sites and ($0$,$l$) the unit-
cell index. From here on, for simplicity, we will drop the unit-cell indexes.
In order to project onto the ligand ion components with specific irreducible
representation coordinates, i.e.
$\Phi_{\alpha\beta}(ij)\rightarrow\Phi_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(ij)$
we apply a projection operator $C_{\Gamma Q}^{\alpha}$ such that $r_{\Gamma
Q}=\sum_{\alpha}C_{\Gamma Q}^{\alpha}r_{\alpha}$.

$\Phi_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(ij)=\frac{\partial^{2}V}{\partial
r_{\Gamma Q}(i)r_{\Gamma Q^{\prime}}(j)}=\sum_{\alpha,\beta}(C_{\Gamma
Q}^{\alpha})^{t}\Phi_{\alpha\beta}(ij)(C_{\Gamma^{\prime}Q^{\prime}}^{\beta})\
.$ (3)

Here $\Gamma\in\\{T_{2g};E_{g}\\}$ while $Q$ runs over the corresponding
distortions $\\{xy,xz,yz;z^{2}x^{2}-y^{2}\\}$. The Jahn-Teller irreducible
representations are constituted from linear combinations of the Cartesian
displacements $X_{i}$, $Y_{i}$, and $Z_{i}$ of the ligand ions situated at the
vertices of the distorted octahedron. These displacements are measured
relative to their positions in an undistorted reference structure, indicated
by $x_{i}^{(0)},y_{i}^{(0)},z_{i}^{(0)}$. Employing $X_{i}$, $Y_{i}$, and
$Z_{i}$, one can construct JT normal coordinates $Q_{i}$, which adhere to the
irreducible representations of the octahedral symmetry group $O_{h}$ Bersuker
(2006).

$A_{1g}$ | $Q_{1}$ | $(X_{2}-X_{5}+Y_{3}-Y_{6}+Z_{1}-Z_{4})/\sqrt{6}$
---|---|---
$E_{g}$ | $Q_{2}$ | $(X_{2}-X_{5}-Y_{3}+Y_{6})/2$
| $Q_{3}$ | $(2Z_{1}-2Z_{4}-X_{2}+X_{5}-Y_{3}+Y_{6})/2\sqrt{3}$
$T_{2g}$ | $Q_{yz}$ | $(Z_{3}-Z_{6}+Y_{1}-Y_{4})/2$
| $Q_{xz}$ | $(X_{1}-X_{4}+Z_{2}-Z_{5})/2$
| $Q_{xy}$ | $(Y_{2}-Y_{5}+X_{3}-X_{6})/2$
Table 1: The octahedral JT coordinates $Q_{i}$ are defined based on the
Cartesian distortions $X_{i}$, $Y_{i}$, and $Z_{i}$ of the system.

We test the consistency of the method in two cases:

  1. 1.

By comparing the $\Gamma-$point projected frequencies against the experimental
and computed ones of a non-magnetic double perovskite containing Re, i.e.
Ba2WReO6.

  2. 2.

By comparing the projected-phonon spectra to the fully computed one for
Ba2MgReO6 .

The implementation of the projective phonon formalism for calculation of EC is
available on GitHub at the following link
https://github.com/dariofiosca/ElastiCouplings.git

###  II.1 Comparison with Ba2WReO6

We start by addressing point (1). Ba$2$WReO$6$ is a non-magnetic double
perovskite whose $\Gamma$-point irreducible phonon frequencies were determined
experimentally through Raman spectroscopy Pásztorová  _et al._ (2023) and via
theoretical first-principles involving both DFT and lattice dynamic
calculations (LDC) Diao  _et al._ (2013); Rodrigues  _et al._ (2018). The
$\Gamma$-point frequencies can be directly derived from equation 1 where
$D_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(0)=1/\sqrt{m_{1}m_{j}}\sum_{l}\Phi_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(0,l)$,
and their comparison is detailed in Suppl. Table 3.

The phonon spectra alongside on-site and EC were calculated subsequent to a
thorough DFT-based relaxation of the volume and internal atomic coordinates,
using the LDA potential. For the relaxation process, a $10\times 10\times 10$
k-mesh was utilized, a 400 eV cutoff for the plane wave basis, and an
electronic self-consistent cycle convergence of 10-4 eV. The final lattice
parameters are listed in Table 3.

For phonon calculations, a $2\times 2\times 2$ supercell was constructed using
the phonopy tool, and the DFPT routine in VASP was employed. The calculations
were performed with an enhanced energy convergence of 10-6 eV and an increased
energy cutoff of 600 eV. Calculations included both on-site, nearest-neighbor
EC, and next-nearest-neighbor EC.

Table 2:  Comparison between lattice constants for Ba2WReO6. The $\dagger$ is used to mark the result of this work.  | DFT $\dagger$ | Exp. 100 K Pásztorová  _et al._ (2023) | Exp. 300 K Pásztorová  _et al._ (2023)
---|---|---|---
$a$ (Å) | 8.1046 | 8.1010 | 8.1120
Table 3:  The frist column represents the computed irreducible frequencies as $\Gamma$ obtained from phonopy, while the second are the ones obtained from the projective method including onsite JT and both nearest and next nearest neighbors EC. The $\dagger$ represents the results of this work. The third colum includes the experimentally defined Raman modes at 300 K and the comparison with experimental data at 300 K on the ceramic material, together with two theoretical calculations of phonon spectra.  R. mod. | Phonopy$\dagger$ | Projected$\dagger$ | Exp. Ref Pásztorová  _et al._ (2023) | Exp Ref. Diao  _et al._ (2013) | DFT Ref. Diao  _et al._ (2013) | LDC Ref. Rodrigues  _et al._ (2018)
---|---|---|---|---|---|---
$E_{g}$ | 534 | 492 | 536 | 538 | 614 | 522
$T_{2g}$ | 421 | 418 | 439 | 441 | 405 | 432
$A_{1g}$ | 780 | 784 | 807 | 812 | 833 | 776

###  II.2 Comparison phonons in BMRO

From the projected force constant matrix
$\Phi_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(ij)$ it is possible to recalculate
the Dynamical Matrix

$D_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(\mathbf{q})=\frac{1}{\sqrt{m_{Q}m_{Q^{\prime}}}}\sum_{ij}\Phi_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(ij)exp(i\mathbf{q}\cdot[\mathbf{r}_{\Gamma}^{Q}(i)-\mathbf{r}_{\Gamma^{\prime}}^{Q^{\prime}}(j)])$
(4)

where ($m_{Q}$, $m_{Q^{\prime}}$) are the oxygen masses, and $\mathbf{q}$ the
wave vector. The equation of motion can be derived to be

$\sum_{\Gamma^{\prime}\gamma^{\prime}}D_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(\mathbf{q})e_{\Gamma^{\prime}}^{Q^{\prime}}(\mathbf{q}\nu)=[\omega(\mathbf{q}\nu)]^{2}e_{\Gamma}^{Q}(\mathbf{q}\nu)$
(5)

and via diagonalization of the equation of motion one obtains the phonon
spectra.

This methodology serves to evaluate the precision of the projective phonon
technique in comparison to the full phonon calculation of the 32 f.u.
supercell. It’s important to highlight that the inclusion of next-nearest
neighbors (NNN) EC is crucial for achieving a high degree of accuracy in the
phonon band structure at the $\Gamma$ and $X$ points.

While incorporating ECs beyond NNN can refine the phonon spectra comparison,
this does not necessarily yield a physically comprehensive solution, as the
octahedra beyond NNN are impacted by the supercell’s periodic boundary
conditions.

We note that the current has its main limitation in the DFPT implementation of
VASP, which only calculates the force constant at the $\Gamma$-point. An
accurate computation of the force constants at other $\mathbf{q}$ points
requires the use of significantly larger supercell structures. However,
scaling up the supercell size poses substantial computational challenges.
Fortunately, for BMRO extensive supercell calculations are not mandatory. As
demonstrated in Suppl. Fig. 3, we achieve satisfactory convergence at the
$\Gamma$ and $X$ high-symmetry points by incorporating on-site, nearest-
neighbor (NN), and NNN intersite elastic couplings.

Figure 3: Comparison between primitive cell phonons of Ba2MgReO6 computed with
phonopy and the reconstructed phonons reconstructed from the projected force
constant matrix $\Phi_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(ij)$. The red
lines mark the $E_{g}$ and $T_{2g}$ phononic bands, while the dotted light
blue are the reconstructed bands from onsite + NN EC and the solid blue line
mark the reconstructed bands from onsite + NN + NNN.

###  II.3 Phonon Band-structures Ba2WReO6 and Ba2MgReO6

Figure 4: Phonon bandstructure of Ba2WReO6. Frequencies are in THz. Figure 5:
Phonon bandstructure of Ba2MgReO6. Frequencies are in THz.

##  III Equation of State

To compute the volumes at different pressure we fit the pressure as function
of volume with the Birch-Murnaghan Equation of State (EOS)

$P=\frac{3}{2}B_{0}\Big{[}\Big{(}\frac{V}{V_{0}}\Big{)}^{-7/3}-\Big{(}\frac{V}{V_{0}}\Big{)}^{-5/3}\Big{]}\times\Big{\\{}1+\frac{3}{4}(B^{\prime}_{0}-4)\Big{[}\Big{(}\frac{V}{V_{0}}\Big{)}^{-2/3}-1\Big{]}\Big{\\}}$
(6)

The volumetric calculations were done in DFT by employing the VASP package
with the 8 f.u. cubic supercell, starting from the fully relaxed solution
described in Section I.4. The LDA potential is used, together with a k-mesh of
$12\times 12\times 12$ and energy cutoff of 600 eV with energy convergence of
the self consistent cycle of 10-5 eV. The results in Table 4 have been
compared with both experimental and first-principles calculation from Ref. and
Ref. Dar  _et al._ (2018) respectively.

Table 4:  Calculated bulk modulus ($B_{0}$ in GPa) and derivative of bulk modulus ($B^{\prime}_{0}$).  | DFT (this work) | DFT Dar  _et al._ (2018)
---|---|---
$B_{0}$ (GPa) | 197.47 | 158.27
$B^{\prime}_{0}$ | 4.51 | 4.80
Table 5:  Calculated lattice constant and bondlength distances.  | DFT (this work) | Exp. Ref Hirai  _et al._ (2020) | DFT Dar  _et al._ (2018)
---|---|---|---
$a$ (Å) | 7.99292 | 8.0802 | 8.1790
Re-O (Å) | 1.92 | 1.96 | 1.95
Mg-O (Å) | 2.08 | 2.08 | 2.14
Ba-O (Å) | 2.83 | 2.85 | 2.89
Figure 6: Fit of Pressure vs volumetric DFT solution with Birch-Murnaghan EOS.

From the EOS we derived for both the Experimental and DFT-relaxed structures,
the evolution of the lattice constant with pressure. In the following table we
report the values used to calculate $H_{IEI}$ and $H_{IEC}$.

Table 6:  Calculated lattice constant from Ba2MgReO6 EOS employed for the calculation of IEI and EC.  Pressure (GPa) | IEI $a$ (Å) | IEC $a$ (Å)
---|---|---
-10 | 8.3058 | 8.22258
-5 | 8.1810 | 8.09898
0 | 8.0802 | 7.99292
5 | 7.9954 | 7.91524
10 | 7.9221 | 7.84270

##  IV Intersite Exchange Interactions

In Suppl. Fig. 7 we plot the IEI in BMRO at zero pressure as a color map; the
IEI values are also listed in Suppl. Table 7 for several values of pressure.
One sees that the quadrupole IEI block includes, apart from antiferroic
$e_{g}$ coupling, also ferroic $t_{2g}$ ones. These $t_{2g}$ ferroic IEI have
not been previously included in model electronic effective Hamiltonians of
$d^{1}$ DP Chen  _et al._ (2010); Svoboda  _et al._ (2021); Iwahara and
Chibotaru (2023); Kubo  _et al._ (2023). The $t_{2g}$ IEI are the reason why
the F$xy$ quarupole order dominates over the AF$x^{2}-y^{2}$ one when the EL
coupling is neglected.

One may also notice significant octupolar IEI, which are the largest couplings
among time-odd IEI.

Figure 7: Color map of the IEI $V_{KK^{\prime}}^{QQ^{\prime}}$ in BMRO at pressure 0 GPa for the [1/2,1/2,0] Re-Re pair. All values are in meV. The numerical list of $V_{KK^{\prime}}^{QQ^{\prime}}$ is given in Suppl. Table 7. Table 7:  Calculated IEI $V_{KK^{\prime}}^{QQ^{\prime}}$ for the $J_{eff}$=3/2 multiplet. First two columns list $Q$ and $Q^{\prime}$, respectively. Third and fourth column displays the $KQ$ and $K^{\prime}Q^{\prime}$ tensors in the Cartesian representation, respectively. The last three columns displays the values of IEI (in meV) for the [1/2,1/2,0] nearest-neighbor Re-Re bond in Ba2MgReO6 for pressure 10 GPa, ambient condition and -16 GPa, respectively.  Pressure | 10 GPa | Ambient | -16 GPa
---|---|---|---
| (Ba2CaReO6 volume)
Dipole-Dipole
-1 | -1 | y | y | -0.79 | -0.59 | -0.31
0 | 0 | z | z | 1.00 | 0.73 | 0.38
1 | -1 | x | y | 0.99 | 0.91 | 0.73
1 | 1 | x | x | -0.79 | -0.59 | -0.31
Quadrupole-Quadrupole
-2 | -2 | xy | xy | -2.07 | -1.59 | -0.94
-1 | -1 | yz | yz | -2.67 | -2.14 | -1.38
0 | -2 | z2 | xy | 1.19 | 1.08 | 0.84
0 | 0 | z2 | z2 | 3.05 | 2.43 | 1.55
1 | -1 | xz | yz | -0.10 | -0.10 | -0.08
1 | 1 | xz | xz | -2.67 | -2.14 | -1.38
2 | 2 | x2-y2 | x2-y2 | -2.36 | -1.89 | -1.23
Octupole-Octupole
-3 | -3 | y(x2-3y2) | y(x2-3y2) | 1.16 | 0.88 | 0.50
-2 | -2 | xyz | xyz | -1.85 | -1.42 | -0.84
-1 | -3 | yz2 | y(x2-3y2) | -0.08 | -0.09 | -0.08
-1 | -1 | yz2 | yz2 | 1.23 | 1.03 | 0.75
0 | -2 | z3 | xyz | -0.64 | -0.60 | -0.49
0 | 0 | z3 | z3 | 3.13 | 2.56 | 1.75
1 | -3 | xz2 | y(x2-3y2) | -0.32 | -0.28 | -0.21
1 | -1 | xz2 | yz2 | -0.83 | -0.76 | -0.61
1 | 1 | xz2 | xz2 | 1.23 | 1.03 | 0.75
2 | 2 | z(x2-y2) | z(x2-y2) | -2.06 | -1.65 | -1.06
3 | -1 | x(3x2-y2) | yz2 | 0.32 | 0.28 | 0.21
3 | 1 | x(3x2-y2) | xz2 | 0.08 | 0.09 | 0.08
3 | 3 | x(3x2-y2) | x(3x2-y2) | 1.16 | 0.88 | 0.50
Dipole-Octupole
-1 | -3 | y | y(x2-3y2) | 0.07 | 0.08 | 0.07
-1 | -1 | y | yz2 | -2.25 | -1.79 | -1.15
-1 | 1 | y | xz2 | -0.23 | -0.21 | -0.17
-1 | 3 | y | x(3x2-y2) | 0.36 | 0.31 | 0.23
0 | -2 | z | xyz | 0.69 | 0.61 | 0.47
0 | 0 | z | z3 | -0.15 | -0.12 | -0.09
1 | -3 | x | y(x2-3y2) | -0.36 | -0.31 | -0.23
1 | -1 | x | yz2 | -0.23 | -0.21 | -0.17
1 | 1 | x | xz2 | -2.25 | -1.79 | -1.15
1 | 3 | x | x(3x2-y2) | -0.07 | -0.08 | -0.07

##  V Intersite Elastic Couplings

Table 8:  The Table displays the calculated values of EC
$\Phi_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(ij)$ for $(ij)$ being the on-site
$(ii)$ component, the irreducible ligand components $(ij)$ connecting NN Re-Re
octahedral centers with $\mathbf{R}=[1/2,1/2,0]$ and the irreducible ligand
components $(ij)$ connecting NNN Re-Re centers with $\mathbf{R}=[1,0,0]$.
First two columns list $Q$ and $Q^{\prime}$, respectively. Units are in eV/Å2.
On site elastic couplings
---
| xy | yz | z2 | xz | x2-y2
| 9.633 | 9.633 | 21.822 | 9.633 | 21.822
Nearest-Neighbors EC for $\mathbf{R}=[0.5,0.5,0]$
| xy | yz | z2 | xz | x2-y2
xy | 0.584 | 0 | -0.463 | 0 | 0
yz | 0 | -0.229 | 0 | -0.250 | 0
z2 | -0.463 | 0 | 0.532 | 0 | 0
xz | 0 | -0.250 | 0 | -0.229 | 0
x2-y2 | 0 | 0 | 0 | 0 | -1.745
Next-Nearest-Neighbors IEC for $\mathbf{R}=[1,0,0]$
| xy | yz | z2 | xz | x2-y2
xy | -0.041 | 0 | 0 | 0 | 0
yz | 0 | 0.027 | 0 | 0 | 0
z2 | 0 | 0 | 0.446 | 0 | -0.082
xz | 0 | 0 | 0 | -0.041 | 0
x2-y2 | 0 | 0 | -0.082 | 0 | 0.541

##  VI Mean-field approach to the IEI+EL effective Hamiltonian

We solve the full many-body effective Hamiltonian (eq. 1 of the main text)
using a generalized single-site mean-field theory. Both the IEI (eq. 2 of the
main text) and EC (eq. 4 ibid) intersite terms are mean-field decoupled; the
resulting single-site JT problem in electronic and elastic mean fields is
treated fully quantum-mechanically. This mean-field approach to the
cooperative JT effect is described in a number of references Englman and
Halperin (1970); Gehring and Gehring (1975); Polinger (2009); here we
implement it for the most general linear JT coupling and mean-field potential
within the $J_{\mathrm{eff}}$=3/2 space. For consistency, we outline the
derivation of this mean-field scheme and fully specify the solution of the
single-site JT problem below.

We label the sites in a given ordered phase by the translational vector
$\bm{l}$ labeling magnetic unit cells and the basis (sublattice) vector
$\bm{\tau}$ within the magnetic unit cell. The intersublattice IEI and EC are

$V_{KK^{\prime}}^{QQ^{\prime}}(\bm{\tau}\bm{\tau}^{\prime})=\sum_{\langle
ij\rangle}\sum_{\bm{l}}V_{KK^{\prime}}^{QQ^{\prime}}(ij)\delta\left({\bf
R}_{ij}-\bm{l}-\bm{\tau}^{\prime}+\bm{\tau}\right)$

and

$\Phi_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(\bm{\tau}\bm{\tau}^{\prime})=\sum_{\langle
ij\rangle}\sum_{\bm{l}}\Phi_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(ij)\delta\left({\bf
R}_{ij}-\bm{l}-\bm{\tau}^{\prime}+\bm{\tau}\right),$

respectively.

The mean-field decoupling of the IEI and EC terms then reads

$H_{\mathrm{IEI}}^{\mathrm{MF}}=\sum_{\bm{\tau}KQ}\left(O_{K}^{Q}(\bm{\tau})-\frac{1}{2}\langle
O_{K}^{Q}(\bm{\tau})\rangle\right)V_{KQ}^{\mathrm{MF}}(\bm{\tau}),$ (7)

$H_{\mathrm{EC}}^{\mathrm{MF}}=\sum_{\tau\Gamma}\sum_{Q\in\Gamma}\left(q_{\Gamma
Q}(\bm{\tau})-\frac{1}{2}\langle q_{\Gamma
Q}(\bm{\tau})\rangle\right)\Phi_{\Gamma Q}^{\mathrm{MF}}(\bm{\tau}),$ (8)

where the respective IEI and EC mean fields are

$V_{KQ}^{\mathrm{MF}}(\bm{\tau})=\sum_{K^{\prime}Q^{\prime}\tau^{\prime}}V_{KK^{\prime}}^{QQ^{\prime}}(\bm{\tau}\bm{\tau}^{\prime})\langle
O_{K^{\prime}}^{Q^{\prime}}(\bm{\tau}^{\prime})\rangle$

and

$\Phi_{\Gamma
Q}^{\mathrm{MF}}(\bm{\tau})=\sum_{\Gamma^{\prime}\tau^{\prime}}\sum_{Q^{\prime}\in\Gamma^{\prime}}\Phi_{\Gamma\Gamma^{\prime}}^{QQ^{\prime}}(\bm{\tau}\bm{\tau}^{\prime})\langle
q_{\Gamma^{\prime}Q^{\prime}}(\bm{\tau}^{\prime})\rangle.$

The single-site JT problem in mean-field at a given site $\bm{\tau}$ is then
given by the single-site Hamiltonian (eq. 3 of the main text) plus the mean-
field terms:

$H_{1s}=\sum_{\Gamma}\sum_{Q\in\Gamma}\left[\frac{p^{2}_{Q}}{2M_{\Gamma}}+\Phi_{\Gamma}\frac{q^{2}_{\Gamma
Q}}{2}-g_{\Gamma}O^{Q}_{2}q_{\Gamma
Q}\right]+H_{\mathrm{IEI}}^{\mathrm{MF}}+H_{\mathrm{EC}}^{\mathrm{MF}},$ (9)

where from now on we omit the $\bm{\tau}$ label for brievity.

We diagonalize the harmonic oscillator terms in (9) and quantize accordingly
the displacements $q_{\Gamma Q}$ in the JT and EC terms. Omitting the constant
$-\frac{1}{2}\langle...\rangle$ terms in (7) and (8), the resulting
Hamiltonian for the single-site JT problem reads

$H_{1s}=\sum_{\Gamma}\sum_{Q\in\Gamma}\left[\hbar\omega_{\Gamma}(\hat{n}_{\Gamma
Q}+1/2)+(\Phi_{\Gamma
Q}^{\mathrm{MF}}-g_{\Gamma}O^{Q}_{2})A_{\Gamma}(a^{\dagger}_{\Gamma
Q}+a_{\Gamma Q})\right]+\sum_{KQ}V_{KQ}^{\mathrm{MF}}O_{K}^{Q},$ (10)

where we introduced creation(annihilation) operators $a^{\dagger}_{\Gamma
Q}$($a_{\Gamma Q}$) for the local octahedral vibration mode $\Gamma Q$,
$\hat{n}_{\Gamma Q}=a^{\dagger}_{\Gamma Q}a_{\Gamma Q}$. The oscillator
strength $\hbar\omega_{\Gamma}$, where
$\omega_{\Gamma}=\sqrt{\Phi_{\Gamma}/M_{\Gamma}}$, calculated with ab initio
on-site elastic couplings $\Phi_{\Gamma}$ (Table 8) is equal to 75.3 (49.8)
meV for the $E_{g}$ ($T_{2g}$) modes, respectively. The quantized
displacements $A_{\Gamma}(a^{\dagger}_{\Gamma Q}+a_{\Gamma Q})=q_{\Gamma Q}$,
where $A_{\Gamma}=\sqrt{\hbar/(2M_{\Gamma}\omega_{\Gamma})}$ is equal to
4.17(5.12) pm for the $E_{g}$($T_{2g}$) modes.

In order to solve the Hamiltonian (10) we introduce a normalized mixed
(vibronic) basis $|\phi_{i}>=|M_{i}\rangle\times|\bm{n}_{i}\rangle$, where
$\langle\phi_{i}|\phi_{j}\rangle=\delta_{ij}$, the electronic state
$|M\rangle\equiv|J_{\mathrm{eff}};M\rangle$ is the eigenstate of
$\hat{J}_{\mathrm{eff}}$ with the projection $M$,
$|\bm{n}\rangle\equiv|n_{E_{g}z^{2}},n_{E_{g}x^{2}-y^{2}},n_{T_{2g}xy},n_{T_{2g}xz},n_{T_{2g}yz}\rangle\equiv\tilde{\bm{n}}_{\Gamma
Q},n_{\Gamma Q}\rangle$ is the local phonon modes occupancies. For
convenience, we introduce here another notation, $\tilde{\bm{n}}_{\Gamma Q}$,
to designate the occupancies of all modes except the $\Gamma Q$ one. In this
basis, the matrix elements of different operators in (10) read

$\displaystyle\langle\phi_{i}|O_{K}^{Q}|\phi_{j}\rangle=\delta_{\bm{n}_{i},\bm{n}_{j}}\left[O_{K}^{Q}\right]_{M_{i}M_{j}},$
$\displaystyle\langle\phi_{i}|\hbar\omega_{\Gamma}(\hat{n}_{\Gamma
Q}+1/2)|\phi_{j}\rangle=\delta_{i,j}\hbar\omega_{\Gamma}(n_{\Gamma Q}+1/2),$
(11) $\displaystyle\langle\phi_{i}|a^{\dagger}_{\Gamma Q}+a_{\Gamma
Q}|\phi_{j}\rangle=\delta_{M_{i},M_{j}}\delta_{\tilde{\bm{n}}_{\Gamma
Q}^{i},\tilde{\bm{n}}_{\Gamma Q}^{j}}\left(\delta_{n_{\Gamma Q}^{i},n_{\Gamma
Q}^{j}+1}\sqrt{n_{\Gamma Q}^{i}}+\delta_{n_{\Gamma Q}^{i}+1,n_{\Gamma
Q}^{j}}\sqrt{n_{\Gamma Q}^{j}}\right).$

We cut the basis $\\{|\phi\rangle\\}$ at a chosen maximum phonon number
$N_{\mathrm{ph}}^{\mathrm{MAX}}=\mathrm{MAX}\left[\sum_{\Gamma Q}n_{\Gamma
Q}\right]$. Once all matrix elements $\langle\phi_{i}|H_{1s}|\phi_{j}\rangle$
have been evaluated with the help of (11), we diagonalize $H_{1s}$ and
recalculate the electronic moments $\langle O_{K}^{Q}\rangle$, displacements
$\langle q_{\Gamma Q}\rangle$, total and free energy for a given site. Having
solved the single-site JT problem for all sites $\bm{\tau}$ in the structure,
we recalculate the exchange (7) and elastic (8) mean fields to close the mean-
field cycle.

This mean-field approach is implemented in the framework of ”McPhase” package
Rotter (2004) as an in-house module. As initial guesses of the mean-field
procedure we employ all 1k structures realizable within the BMRO single
conventional unit cell that contains 4 magnetic Re sites.

##  VII Effect of the vibronic basis cutoff and Jahn-Teller regime

In order to benchmark the convergence of the numerical solution of the JT
problem (10) with respect to the vibronic basis cutoff
$N_{\mathrm{ph}}^{\mathrm{MAX}}$, we solved this problem, without the mean-
field terms, varying $N_{\mathrm{ph}}^{\mathrm{MAX}}$. The JT problem in this
case reduces to the single-site term (eq. 3 of the main text), which in the
quantized form reads

$H_{\mathrm{JT}}=\sum_{\Gamma}\sum_{Q\in\Gamma}\hbar\omega_{\Gamma}(\hat{n}_{\Gamma
Q}+1/2)-g_{\Gamma}O^{Q}_{2}A_{\Gamma}(a^{\dagger}_{\Gamma Q}+a_{\Gamma Q}).$
(12)

We employ the values of parameters ($g_{\Gamma}$, $\omega_{\Gamma}$, and
$A_{\Gamma}$) that are calculated for BMRO at ambient pressure (see the
previous section and Suppl. Fig. 2b).

Figure 8: Total energy (a) and entropy (b) of the single-site JT problem
ground state as a function of the basis cutoff
$N_{\mathrm{ph}}^{\mathrm{MAX}}$. The size of the corresponding basis is
indicated by the upper tick labels of panel a.

In Suppl. Figs. 8a and 8b we display the calculated total energy (omitting the
zero-point contribution) and entropy vs. $N_{\mathrm{ph}}^{\mathrm{MAX}}$. One
sees that already for $N_{\mathrm{ph}}^{\mathrm{MAX}}=1$, corresponding to the
weak coupling treatment of the JT problem, the total energy is quite well
converged, while the entropy still corresponds to the ground-state quadruplet
and is not affected by including local phonons into the basis. In order to
explore the limits of the weak-coupling treatment, we also ran the same
calculation with the JT coupling strength artificially increased. With the JT
coupling scaled by the factors 1.5 and 2, the weak-coupling treatment clearly
becomes insufficient for the total energy. The ground-state degeneracy
markedly increases once the JT coupling becomes twice larger, indicating that
the JT ground-state solutions are now localized within the ”trough” of the JT
potential energy surface Bersuker and Polinger (1989).

##  VIII Relaxed BMRO structure

    # CIF file created by FINDSYM, version 7.1.3

    data_findsym-output
    _audit_creation_method FINDSYM

    _cell_length_a     7.9991900000
    _cell_length_b     7.9991900000
    _cell_length_c     7.9991900000
    _cell_angle_alpha  90.0000000000
    _cell_angle_beta   90.0000000000
    _cell_angle_gamma  90.0000000000
    _cell_volume       511.8444957459

    _symmetry_space_group_name_H-M "F 4/m -3 2/m"
    _symmetry_Int_Tables_number 225
    _space_group.reference_setting ’225:-F 4 2 3’
    _space_group.transform_Pp_abc a,b,c;0,0,0

    _atom_site_label
    _atom_site_type_symbol
    _atom_site_symmetry_multiplicity
    _atom_site_Wyckoff_symbol
    _atom_site_fract_x
    _atom_site_fract_y
    _atom_site_fract_z
    _atom_site_occupancy
    _atom_site_fract_symmform
    Mg1 Mg   4 a  0.00000  0.00000  0.00000  1.00000 0,0,0
    Re1 Re   4 b  0.50000  0.50000  0.50000  1.00000 0,0,0
    Ba1 Ba   8 c  0.25000  0.25000  0.25000  1.00000 0,0,0
    O1  O   24 e  0.25951  0.00000  0.00000  1.00000 Dx,0,0

    # end of cif

## References

  * Hubbard (1963) J. Hubbard, Proc. Roy. Soc. (London) A 276, 238 (1963).
  * Aichhorn  _et al._ (2009) M. Aichhorn, L. Pourovskii, V. Vildosola, M. Ferrero, O. Parcollet, T. Miyake, A. Georges, and S. Biermann, Phys. Rev. B 80, 085101 (2009).
  * Aichhorn  _et al._ (2011) M. Aichhorn, L. Pourovskii, and A. Georges, Phys. Rev. B 84, 054529 (2011).
  * Aichhorn  _et al._ (2016) M. Aichhorn, L. Pourovskii, P. Seth, V. Vildosola, M. Zingl, O. E. Peil, X. Deng, J. Mravlje, G. J. Kraberger, C. Martins, _et al._ , Computer Physics Communications 204, 200 (2016).
  * Blaha  _et al._ (2018) P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka, J. Luitz, R. Laskowski, F. Tran, and L. D. Marks,  _WIEN2k, An augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties_ (Karlheinz Schwarz, Techn. Universität Wien, Austria,ISBN 3-9501031-1-2, 2018).
  * Fiore Mosca  _et al._ (2021) D. Fiore Mosca, L. V. Pourovskii, B. H. Kim, P. Liu, S. Sanna, F. Boscherini, S. Khmelevskyi, and C. Franchini, Phys. Rev. B 103, 104401 (2021).
  * Pourovskii  _et al._ (2021) L. V. Pourovskii, D. F. Mosca, and C. Franchini, Phys. Rev. Lett. 127, 237201 (2021).
  * Pourovskii (2016) L. V. Pourovskii, Phys. Rev. B 94, 115117 (2016).
  * Liechtenstein  _et al._ (1987) A. Liechtenstein, M. Katsnelson, V. Antropov, and V. Gubanov, Journal of Magnetism and Magnetic Materials 67, 65 (1987).
  * Katsnelson and Lichtenstein (2000) M. I. Katsnelson and A. I. Lichtenstein, Phys. Rev. B 61, 8906 (2000).
  * Pourovskii (2023) L. V. Pourovskii, Phys. Rev. B 108, 054436 (2023).
  * Campbell  _et al._ (2006) B. J. Campbell, H. T. Stokes, D. E. Tanner, and D. M. Hatch, J. App. Crystal. 39, 607 (2006).
  * Stokes  _et al._ (2022) H. T. Stokes, D. M. Hatch, and B. J. Campbell, “ISODISTORT, ISOTROPY Software Suite,” https://iso.byu.edu (2022).
  * Momma and Izumi (2011) K. Momma and F. Izumi, Journal of Applied Crystallography 44, 1272 (2011).
  * Delange  _et al._ (2017) P. Delange, S. Biermann, T. Miyake, and L. Pourovskii, Phys. Rev. B 96, 155132 (2017).
  * Brooks  _et al._ (1997) M. Brooks, O. Eriksson, J. Wills, and B. Johansson, Physical review letters 79, 2546 (1997).
  * Iwahara and Chibotaru (2023) N. Iwahara and L. F. Chibotaru, Phys. Rev. B 107, L220404 (2023).
  * Soh  _et al._ (2023) J.-R. Soh, M. E. Merkel, L. Pourovskii, I. Živković, O. Malanyuk, J. Pásztorová, S. Francoual, D. Hirai, A. Urru, D. Tolj, D. Fiore-Mosca, O. Yazyev, N. A. Spaldin, C. Ederer, and H. M. Rønnow, “Spectroscopic signatures and origin of a hidden order in ba2mgreo6,” (2023), arXiv:2312.01767 [cond-mat.str-el] .
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Togo  _et al._ (2023) A. Togo, L. Chaput, T. Tadano, and I. Tanaka, J. Phys. Condens. Matter 35, 353001 (2023).
  * Togo (2023) A. Togo, J. Phys. Soc. Jpn. 92, 012001 (2023).
  * Bersuker (2006) I. Bersuker,  _The Jahn-Teller Effect_ (Cambridge University Press, Cambridge, 2006).
  * Pásztorová  _et al._ (2023) J. Pásztorová, W. H. Bi, R. Gaal, K. Krämer, I. Živković, and H. M. Rønnow, Journal of Solid State Chemistry 326, 124184 (2023).
  * Diao  _et al._ (2013) C.-L. Diao, C.-H. Wang, N.-N. Luo, Z.-M. Qi, T. Shao, Y.-Y. Wang, J. Lu, F. Shi, and X.-P. Jing, Journal of the American Ceramic Society 96, 2898 (2013), https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.12429 .
  * Rodrigues  _et al._ (2018) J. Rodrigues, D. Bezerra, and A. Hernandes, Journal of Raman Spectroscopy 49, 1822 (2018), https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/pdf/10.1002/jrs.5460 .
  * Dar  _et al._ (2018) S. A. Dar, V. Srivastava, U. K. Sakalle, and G. Pagare, Computational Condensed Matter 14, 137 (2018).
  * Hirai  _et al._ (2020) D. Hirai, H. Sagayama, S. Gao, H. Ohsumi, G. Chen, T.-h. Arima, and Z. Hiroi, Phys. Rev. Res. 2, 022063 (2020).
  * Chen  _et al._ (2010) G. Chen, R. Pereira, and L. Balents, Phys. Rev. B 82, 174440 (2010).
  * Svoboda  _et al._ (2021) C. Svoboda, W. Zhang, M. Randeria, and N. Trivedi, Phys. Rev. B 104, 024437 (2021).
  * Kubo  _et al._ (2023) H. Kubo, T. Ishitobi, and K. Hattori, Phys. Rev. B 107, 235134 (2023).
  * Englman and Halperin (1970) R. Englman and B. Halperin, Phys. Rev. B 2, 75 (1970).
  * Gehring and Gehring (1975) G. A. Gehring and K. A. Gehring, Reports on Progress in Physics 38, 1 (1975).
  * Polinger (2009) V. Polinger, “Orbital ordering versus the traditional approach in the cooperative jahn–teller effect: A comparative study,” in  _The Jahn-Teller Effect: Fundamentals and Implications for Physics and Chemistry_, edited by H. Köppel, D. R. Yarkony, and H. Barentzen (Springer Berlin Heidelberg, Berlin, Heidelberg, 2009) pp. 685–725.
  * Rotter (2004) M. Rotter, Journal of Magnetism and Magnetic Materials 272-276, Supplement, E481 (2004).
  * Bersuker and Polinger (1989) I. B. Bersuker and V. Z. Polinger,  _Vibronic Interactions in Molecules and Crystals_ (Springer Berlin Heidelberg, 1989).


# A Framework for Guided Motion Planning

Amnon Attali  Stav Ashur  Isaac Burton Love  Courtney McBeth  James Motes
Marco Morales  Nancy M. Amato

(September 2022)

###### Abstract

Randomized sampling based algorithms are widely used in robot motion planning
due to the problem’s intractability, and are experimentally effective on a
wide range of problem instances. Most variants bias their sampling using
various heuristics related to the known underlying structure of the search
space. In this work, we formalize the intuitive notion of _guided search_ by
defining the concept of a _guiding space_. This new language encapsulates many
seemingly distinct prior methods under the same framework, and allows us to
reason about guidance, a previously obscured core contribution of different
algorithms. We suggest an information theoretic method to evaluate guidance,
which experimentally matches intuition when tested on known algorithms in a
variety of environments. The language and evaluation of guidance suggests
improvements to existing methods, and allows for simple hybrid algorithms that
combine guidance from multiple sources.

##  I Introduction

A general paradigm in Artificial Intelligence is that one can solve provably
worst case intractable search problems using heuristics that are effective in
the average case, i.e., in most real world problems. This same principle holds
true in robotics, where much of the difficulty stems from the geometry and the
implicit nature of the search space, called the configuration space (c-space).
The family of sampling based motion planning (SBMP) algorithms searches the
continuous c-space using random sampling and local extensions in order to find
collision free paths.

Most such algorithms implicitly perform guided search. Early algorithms, such
as Rapidly-exploring Random Trees (RRT) [23], showed that exploration through
Voronoi bias, i.e., choosing which node to expand with probability
proportional to the volume of its Voronoi cell, is highly effective in the
absence of arbitrarily narrow passages. As such, much subsequent work has
centered around the narrow passage problem [3, 28]. Due to the intractable
nature of motion planning, tangible improvements often came from making
assumptions about the underlying space, such that the human engineer could
encode case specific heuristics into the search. Some work, rather than
improve runtime, focused attention on the types of solutions output by the
algorithm, for example searching for paths with high clearance [16].

In this work we make this previously implicit notion of guidance explicit by
formally defining the concepts of _guiding space_ and _guided search_. We
define a guiding space as an auxiliary space that helps estimate the value of
exploring each configuration. This value helps our proposed general guided
search algorithm determine its next exploration step. Our definition of
guidance is naturally hierarchical, as any search space (e.g., c-space) can
itself be used as a guiding space. Intuitively, guidance is the bias
introduced into the exploration process. Our algorithm is similar to A*, which
performs heuristic guided search in _discrete_ spaces, but does so in c-space.

In formally discussing guiding spaces we make the source of guidance explicit.
Doing so allows us to identify guidance-generating components in many
seemingly distinct prior methods, which fit under three main subcategories;
robot modification, environment modification, and experience based guidance.
The perspective of guiding spaces helps isolate conceptual contributions of
different algorithms, and consequently yields simpler generalized
implementations that are also easily composed. Framing SBMP algorithms in this
common language enables us to define a new method for evaluation that focuses
on the quality of guidance.

We approach the task of evaluating guidance from an information theoretic
perspective, whereby the procedure for computing guidance conveys information
to the search process. We test our metric experimentally on a range of
algorithms and environments. The results match intuition for the suitability
of each algorithm to features in certain c-spaces. Moreover, the metric
captures properties of algorithms otherwise obscured by traditional metrics
such as runtime or number of samples. Measuring quality of guidance locally
gives insight into the value of combining guiding spaces, which we demonstrate
with a hybrid guiding space that outperforms its individual parts.

##  II Preliminaries

A motion planning problem consists of the tuple $M=(r,E,s,t)$, for robot $r$,
environment $E$, and task $(s,t)$, where the objective is to find the shortest
collision free path $P_{s,t}$ from starting configuration $s$ to goal
configuration $t$. Note that $M$ is the single-query variant of the motion
planning problem, and we will touch on the multi-query version later in this
section. Together, $(r,E)$ describe the configuration space $\mathcal{C}$, a
continuous space consisting of all possible configurations of $r$ in $E$. This
space is implicitly partitioned into two parts, $C_{free}$ and $C_{obs}$.
$C_{free}$ is the set of all configurations in which $r$ does not intersect
with either itself or obstacles in $E$, and
$C_{obs}=\mathcal{C}\setminus\mathcal{C}_{free}$.

SBMP methods, such as RRT, explore $\mathcal{C}$ to solve $M$ by building a
_search tree_ $T$. They do so by iteratively picking an existing node $v\in
T$, and then expanding the tree by adding a new node as a child of $v$. We
refer to these two operations as _node selection_ , which is performed by a
sampler, and _node expansion_ , which is performed by a local planner. We
refer to the set of points reachable by the local planner as
$N_{\mathcal{C}}(v)$ for $v\in T$, the neighborhood of $v$. Usually the
initial tree contains only the start node $s$ and the algorithm terminates
when a configuration sufficiently close to $t$ is added to the tree.

Note that in this work we depart from the classical partition into tree and
roadmap algorithms (e.g., PRM [22]). Our perspective is that roadmap based
algorithms (which build a general graph structure instead of a search tree)
preprocess $\mathcal{C}$ to produce a data structure that solves the multi-
query version of the problem. In particular, most variants of roadmap
algorithms do not depend on the task $(s,t)$. In other words, roadmap
algorithms do not solve single task search problems such as $M$ but rather
provide guidance for doing so in a given c-space, where the roadmap is the
guiding space. This perspective motivates our definition of guided motion
planning: a tree (or forest) based search algorithm that uses a guiding space
to inform exploration.

##  III Defining Guided Planning

(a) A 2D workspace with a stationary 2-joint manipulator, and 2 obstacles.
Three intermediates from the path found in the guiding space (b) are shown in
dashed gray, where the last is the target configuration. We observe that one
intermediate is invalid.

(b) The lazy planning guiding space for the environment (a) without obstacle
$o_{2}$. A red curve shows a path between the current and target
configurations shown in (a), and a blue wedge shows the resulted sampling
bias.

Figure 1: An illustration of guiding space usage.

We define a _search space_ $\mathcal{X}$ as a state space (set) with a
distance function
$d_{\mathcal{X}}:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}_{+}$. In
contrast, we define a _guiding space_ $\mathcal{S}$ as a state space which
assigns a value to each state, namely there exists some function
$h_{\mathcal{S}}:\mathcal{S}\rightarrow\mathbb{R}_{+}$, which we call the
_heuristic_. Note that a c-space $\mathcal{C}$ is a search space and that
every search space can be viewed as a guiding space, i.e., $h_{t}(x)=d(x,t)$
for $x\in\mathcal{C}$. In other words, distance to some state (goal) $t$ can
be used as an heuristic.

Consider the motion planning problem $M=(r,E,s,t)$ depicted in Figure 1, where
a two joint planar manipulator $r$ must navigate from state $s$ to $t$ while
avoiding obstacles $E=\\{o_{1},o_{2}\\}$. Now consider the related c-space in
Figure 1(b), $\mathcal{C}^{\prime}=(r,E^{\prime}),E^{\prime}=\\{o_{1}\\}$. Let
$\mathcal{S}$ be the guiding space whose heuristic is the distance to $t$ in
$\mathcal{C}^{\prime}$, i.e.,
$h_{\mathcal{S}}(x)=d_{\mathcal{C}^{\prime}}(x,t)$. As we will show next, this
guiding space can be used to guide planning in $\mathcal{C}$.

A _guiding space projection_
$P_{T}:\mathcal{X}\times\mathcal{X}\rightarrow\mathcal{S}$ is a mapping from
pairs of states in some search space $\mathcal{X}$ to states of a guiding
space $\mathcal{S}$. This mapping is conditioned on a search tree $T$, which
intuitively represents the known information regarding $\mathcal{C}_{free}$
and $\mathcal{C}_{obs}$. By relating tasks in $\mathcal{X}$ to states in
$\mathcal{S}$ the projection allows a guided search algorithm to use the
heuristic during search. Often for convenience we will drop the second term to
the projection, namely $P_{T}(x)=P_{T}(x,t)$.

In the example of Figure 1 the projection is
$P_{T}((x,t)\in\mathcal{C}\times\mathcal{C})=x\in\mathcal{C}^{\prime}$. Note
that this projection does not change the actual configuration $x$, and only
modifies its validity since configurations which collide only with obstacle
$o_{2}$ in $\mathcal{C}$ are valid in $\mathcal{C}^{\prime}$. Moreover, notice
that despite colliding with obstacle $o_{2}$, the path highlighted in red
still contains useful information for computing paths in $\mathcal{C}$. First,
the length of this path is a good estimate for the true global shortest path
length, since a slight deviation is collision free and optimal. Second,
regardless of global properties of the path, at certain configurations it
contains locally near-optimal controls. The key insight is that not only does
this path contain useful information, but also it is easier to compute than
the true optimal solution, since it corresponds to a simpler search problem.
The existence of quick to compute useful information is what motivates our
definition of guided motion planning.

A _guided motion planning problem_ $M_{P}=(r,E,s,t,P)$ is a motion planning
problem augmented with a projection $P$. We note that the image of $P$
describes a (potentially implicit) guiding space $\mathcal{S}$ and
corresponding $h$. Next we formalize both how this information is used and
where it comes from.

input : $M_{P}=(r,E,s,t,P)$

1 Initialize tree $T=\\{s\\}$

2 while  _$\neg$ found goal_ do

3 select node: $x_{min}=\operatorname*{argmin}_{v\in T}h(P_{T}(v))$

4 expand node: $x_{new}=\operatorname*{argmin}_{v\in
N_{\mathcal{C}}(v)}h(P_{T}(v))$

5 $is\\_valid=(x_{min},x_{new})\subset\mathcal{C}_{free}$

6 if $is\\_valid$: add $x_{new}$ to $T$ as child of $x_{min}$

7 else: add $(x_{min},x_{new})$ to $T.failed\\_expansions$

8 end while

Algorithm 1 Guided Search

###  III-A Guided Search

Given a search tree $T$, a baseline exploration algorithm might explore the
space by selecting an existing state on the tree uniformly at random and then
expanding from it according to some local planner. Unfortunately this approach
has poor exploration properties, as it is heavily biased towards regions that
have already been explored. Discrete search avoids this problem as newly
explored states are checked against a cached set of visited states.

Algorithm 1 shows a guided search procedure for solving a problem $M_{P}$. In
each iteration the algorithm maps nodes on the tree through the projection
$P_{T}$ then uses $h$ to determine the next node to explore, i.e., by
selecting the node with minimum $h$. In general, we can think of the
projection as providing guidance by computing a distribution over the nodes of
the tree from which we can sample to select how to expand. The tree then grows
according to the neighborhood defined by the local planner. This parallels the
well known heuristic guided discrete search algorithm A*.

###  III-B Computing a Guiding Space Projection $P$

Consider an unknown but fixed distribution $\mathcal{D}$ over robots,
environments, and tasks. Guided search consists of an initial _offline_
training phase, during which we train a projection $P^{M}_{T}$ on samples
$\\{M~{}|~{}M\sim\mathcal{D}\\}$ from this distribution. Note we are generally
not concerned with time complexity during training, and that many algorithms
skip this phase entirely. For example, RRT does not “train”, and performs the
same procedure regardless of c-space. More commonly, after making assumptions
about the c-space distribution, the human engineer encodes their prior
knowledge into a suitable procedure for computing the projection, thus
replacing the training phase. For example knowing that most obstacles in the
environment are fixed might motivate projecting onto a guiding space in which
collision with those obstacles has been precomputed [21]. Alternatively,
knowing that some degree of freedom has negligible impact on $\mathcal{C}$
might motivate projecting onto the other dimensions [24].

Next begins the _online_ phase, in which a specific instantiation
$M\sim\mathcal{D}$ of a motion planning problem is sampled from the
distribution and we conduct guided search on $M_{P_{T}}$ following Algorithm
1.

###  III-C Simple Examples of Guiding Spaces

One guiding space for $\mathcal{C}$ is itself, $\mathcal{S}=\mathcal{C}$,
where $P_{T}(x)=x$ is the identity and $h=d_{\mathcal{C}}(x,t)$ is the true
distance to the goal configuration. This highlights the general setting in
which a search space is used as a guiding space, in which for guidance to be
useful we expect that solving tasks in $\mathcal{S}$ be easier than in
$\mathcal{C}$. Similarly, $\mathcal{S}=\\{0\\}$ is a guiding space with a
trivial projection $P$ and uninformed $h$. Such trivial projections correspond
to fixed sampling distributions over $T$. We observe a natural trade-off
between guidance quality and computation difficulty - computing solutions when
$\mathcal{S}=\mathcal{C}$ is difficult (the original problem) but they provide
perfect guidance, whereas when $\mathcal{S}=\\{0\\}$ we can compute $h$ in
constant time but get uninformed guidance.

The classic implementation of RRT samples uniformly in $\mathcal{C}$ to select
which node on the tree to expand. This is clearly not equivalent to sampling a
node uniformly at random from the tree, and is known to consequently guide
exploration towards unexplored regions via Voronoi bias. While in practice the
guiding space is implicit in the sampling and nearest neighbor procedure, we
can compute it explicitly through the Voronoi decomposition, where the value
$h$ of each node on the tree is proportional to the volume of its
corresponding cell. Notice how the guiding space depends on the search tree,
as the nodes dictate the Voronoi cells.

As noted previously, a PRM in configuration space is a guiding space. If the
distribution $\mathcal{D}$ contains a single robot and environment, but is
stochastic over tasks, then precomputing a discrete roadmap in $\mathcal{C}$
is a common approach for answering multiple queries. The projection $P$ is
often taken to be the function which returns the nearest neighbor on the
graph. For any task $(s,t)$ and vertex $v=P(s)$, the heuristic
$h_{\mathcal{S}}(v)$ is the length of the shortest path on the graph between
$v$ and $P(t)$. This length is an approximation to the optimal solution for
that task and can sometimes provide imperfect guidance; for example not every
point in space might be locally navigable to a node on the graph.

##  IV Related work: categories of guidance

In this section we discuss existing work which does guided search. Note that
many of the works cited below do not discuss guidance as a contribution, yet
here we identify that they do implicitly perform guided search.

We separate guiding space methods into three main categories: robot
modification, environment modification and experience based guidance. We find
these subcategories natural and comprehensive, as the entire motion planning
problem is encapsulated in the structure of the configuration space, which can
be learned from experience, or is induced by its robot and environment
constraints.

###  IV-A Robot Modification

Recalling that a configuration space is the product of a robot $r$ and
environment $E$, $\mathcal{C}=(r,E)$, we define robot modification as
producing a guiding space $\mathcal{S}=(r^{\prime},E)$ for some new robot
$r^{\prime}$. In this paradigm, the projection function $P$ maps
configurations in the first robot’s c-space to configurations in the second.

Given such a mapping $P$, $h$ generally computes paths between configurations
in $\mathcal{S}$. Notice that the translation of a path back into
$\mathcal{C}$ is often costly and may be one-to-many (e.g., inverse-
kinematics), but isn’t necessary because we are only interested in the _value_
of each state. It is then the job of the guided search algorithm to translate
this value into feasible paths in $\mathcal{C}$.

There are many natural ways one might modify a robot description to obtain a
closely related configuration space. One of the most natural projections is to
define $\mathcal{S}$ as workspace itself, ignoring the robot entirely and
treating it as a fully controllable point robot in the environment [19, 13].
Some prior work ignores kinodynamic constraints on the robot [26], thereby
yielding a guiding space where the robot is locally controllable. For chain-
link robots, such as manipulators, a common guiding space involves focusing on
a subset of the degrees of freedom [24]. Such methods often use heuristics in
deciding which subsets to consider, e.g. focus on the first joints of a
manipulator. Finally, some work attempts to learn a subspace of the robot’s
degrees of freedom, e.g. using statistical projections such as PCA [25].

In general, we may think of robot modification methods as performing transfer
learning in which examples of one robot performing tasks are used to inform
plans for a second robot [14]. In the absence of a ground truth robot model,
one can be learned from data [9].

All the above robot modification methods produce a search space as opposed to
a guiding space. For example, workspace guidance methods often build a sparse
graph in workspace, such as the medial axis [13]. Therefore there is often a
secondary guided search that occurs in the modified robot space, e.g., the
medial axis graph. This both highlights the hierarchical nature of guidance
and helps isolate contributions - the medial axis is a useful structure for
guiding search in workspace, and workspace is useful for guiding search in
c-space.

###  IV-B Environment Modification

Environment modification refers to guiding spaces which map
$\mathcal{C}=(r,E)$ to $\mathcal{S}=(r,E^{\prime})$. Unlike in robot
modification, the associated projection $P$ is often simpler to compute, and
translation of paths in $\mathcal{S}$ back to $\mathcal{C}$ requires
relatively little work. In fact, in all work mentioned below, $P$ is nothing
but the identity function, which only changes the validity of some points in
$\mathcal{C}$ to produce $\mathcal{S}$.

Narrow passages are a known bottleneck for randomized sampling algorithms. As
such, some works modify obstacles in an attempt to widen passages [6, 29].
Lazy planning methods [10, 18] also perform a type of environment modification
by ignoring constraints during initial planning, and incorporating them later
as needed. When removing constraints, paths in $\mathcal{S}$ are not
necessarily valid in $\mathcal{C}$, and thus often such methods iteratively or
hierarchically compute a guiding space (e.g., the initial lazy graph), use $h$
to produce samples near the unvalidated path, and then repeat with a new
guiding space informed by the recently acquired samples (e.g., remove edges
found to be invalid).

In contrast to the removal of environment constraints, adding constraints is
less common, but can be done to simplify geometries and therefore collision
detection [17]. Moreover, paths found when adding constraints are always valid
in the original search space, and are thus easy to reuse. In other words, if
we take $h$ to be the cost of paths in the modified space, constraint removal
provides optimistic lower bounds whereas constraint addition produces
conservative upper bounds.

###  IV-C Experience Based Guidance

The most common version of experience based guidance assumes the dataset
contains a variety of environments all for the same robot. Consequently,
database queries often consist of evaluating the degree to which paths in the
dataset violate discovered constraints in the current motion planning instance
[7]. While some methods only query using the task endpoints [12], more recent
approaches propose to query (or filter) the database using local task features
derived from workspace occupancy [11]. Another approach projects trajectories
onto a lower dimensional abstract space via auto-encoding [20] or metrics like
dynamic time warping [1], which improves database query.

Some methods perform multiple queries, where the guiding space effectively
becomes the composite space of multiple independent motion planning problems
from the dataset, e.g., through gaussian mixture models [11].

At the extreme end we have learning based methods which train a neural network
from a dataset of experience to provide guidance, effectively combining all
past solutions through the weights of the neurons of a neural network. In
encoder-decoder models [27] we may think of the encoder as $P$, which maps the
input configuration space onto an abstract guiding space, and the decoder as
$h$, which solves the task in the abstract space to provide guidance, say
through predicting cost-to-go [8]. Many learning based methods have poor
generalization due to dependence on problem specific features, e.g., top down
2D occupancy maps [30]. Those with good generalization often scale poorly by
learning an equivalent to cost-to-go [8], which even in a discrete space
corresponds to computing or storing a prohibitively expensive all pairs
shortest path.

In all learning based approaches a major challenge is gathering appropriate
data, leading to methods which simplify the problem by learning local guidance
through learned local planners [15]. Some work [5] proposes to gather
appropriate data for training by using probabilistic roadmaps.

(a)

(b)

(c)

(d)

Figure 2: The environments in which we evaluate the different guiding space
algorithms. We draw the robot in its start and goal configurations in red.

(a)

(b)

(c)

(d)

Figure 3: Example guiding spaces of each algorithm. We display the starts and
goals of the rectangular robot for RRT, LazyPRM, and PathDatabase to indicate
that these guiding spaces are in c-space, whereas MedialAxis is in workspace.
For PathDatabase we display the paths from the database that pass near the
goal.

##  V Evaluation metric

###  V-A Sampling Efficiency

Often we evaluate sampling based motion planning methods using algorithm level
holistic metrics such as running time or number of samples. While minimizing
such metrics is a reasonable end goal, they can be sensitive to implementation
details, and focusing on them can obscure which aspect of a proposed
contribution is responsible for the results. In the context of guided motion
planning, we now present a new method to evaluate quality of guidance.

In this work we focus on _node selection_ as opposed to node expansion, and
note that evaluating node expansion can be reduced to evaluating a node
selection procedure between the set of possible expansions at a given node.

We propose to evaluate guidance by comparing the distribution of samples
produced by the guiding space to some target distribution, namely one derived
from an optimal sampling distribution with oracle access.

######  Definition 1 (Sampling efficiency)

Given a tree $T$, a target sampling distribution over the tree nodes $Q$, and
an empirical sampling distribution $P$, the sampling efficiency of $P$ is
defined as the information gain of using the target $Q$ instead of $P$,
equivalently described as the Kullback-Liebler divergence between the two
distributions,

$SE_{Q}(P)=D_{KL}(P||Q)=\sum_{v\in T}P(v)\log\frac{P(v)}{Q(v)}$

When $P$ and $Q$ are very similar $SE_{Q}(P)$ is close to zero, and grows as
they grow far apart. In the case where our empirical distribution is
deterministic, $P(u)=1,P(v)=0,v\neq u$, we have,

$SE_{Q}(P)=\sum_{v}P(v)\log\frac{P(v)}{Q(v)}=\log\frac{1}{Q(u)}=-\log Q(u),$

reducing to the negative log-likelihood of the empirical sample under the
target distribution.

Intuitively, we aim to measure the difference between these two distributions,
and while the KL-Divergence provides such a measure, it can be dominated by
small regions in space where the two distributions have largest difference
(e.g., when one distribution assigns zero probability to a region in which the
other is non-zero). In sampling based motion planning we often aren’t worried
about a few bad samples, so long as good ones are made with high enough
probability so the bounded and symmetric Jensen-Shannon Divergence metric
provides a smoother comparison, i.e.,

$D_{JS}(P||Q)=\frac{1}{2}\left[D_{KL}(P||M)+D_{KL}(Q||M)\right],$

for mixture distribution $M=\frac{P+Q}{2}$. $D_{JS}$ is therefore useful as it
is less influenced by outliers in the distributions. Alternatively we can
assume the target $Q$ assigns a minimum nonzero probability to every node,
$Q(v)\geq\epsilon$ for all $v\in T$, and therefore
$SE_{Q}(P)\leq\log\frac{1}{\epsilon}$. In this way $\log\frac{1}{\epsilon}$
captures the maximum penalty we associate with poor sampling decisions made by
the guiding space.

In shortest path problems we are interested in (efficiently) computing a
boundedly optimal $(s,t)$ path in $\mathcal{C}$. Let $\delta_{v}$ be the
suboptimality of $v\in T$,
$\delta_{v}=d_{T}(s,v)+d_{\mathcal{C}}(v,t)-d_{\mathcal{C}}(s,t)>0$, where
$d_{T}$ represents the distance between nodes on the tree and
$d_{\mathcal{C}}$ is the true shortest path between configurations in
$\mathcal{C}$. Now let $\tau_{v}$ be the work remaining for a search starting
at $v$, $\tau_{v}=d_{\mathcal{C}}(v,t)$. We propose to define $Q$ with respect
to these two quantities:

$Q_{T}(v)\propto\exp(-\frac{\delta_{v}}{\delta})\exp(-\frac{\tau_{v}}{\tau})$
$=\frac{e^{-(\delta_{v}/\delta+\tau_{v}/\tau)}}{Z},Z=\sum_{x}e^{-(\delta_{x}/\delta+\tau_{x}/\tau)}$

i.e., $Q_{T}(v)=softmin(\delta_{v}/\delta+\tau_{v}/\tau)$ for temperature
parameters $\delta,\tau$ which naturally balance the relative importance of
optimality and efficiency. The smaller these parameters the steeper the
distribution, assigning higher probability to samples with low $\delta_{v}$
and $\tau_{v}$, and more heavily penalizing poor samples.

We can now define the smoothed $\bar{Q}_{T}\geq\epsilon$ as

$\bar{Q}_{T}(v)\propto\exp(-(\delta_{v}/\delta+\tau_{v}/\tau))+\gamma,$

where

$\min_{u}\frac{e^{-(\delta_{u}/\delta+\tau_{u}/\tau)}+\gamma}{Z+\gamma|T|}\geq\epsilon\implies\gamma\geq\max_{u}\frac{Z(\epsilon-
Q_{T}(u))}{1-\epsilon|T|}$

(a)

(b)

(c)

(d)

Figure 4: We show examples of each guidance method on the main environment it
was tested on. Note that the visualizations are the 2D projections of the 3D
c-space of a planar rectangular robot. In all cases, our implementation
queries the local planner for valid edges to the goal, thus terminating the
search once a configuration has been found that is visible from the goal. The
final paths are highlighted in red.

###  V-B Axes for Evaluating Guidance

The metric described above measures the quality of each individual node
expansion step of the guided search algorithm, but there are other options for
evaluating guidance.

####  V-B1 Environment level

One can measure quality of guidance globally or in terms of local (possibly
hierarchical) environment regions. In global evaluation we take an expectation
of quality over tasks sampled from the environment. In local evaluation we
subdivide the environment and compute our metric on each region independently,
reporting the distribution of results over the full environment. Note that our
method for global evaluation naturally applies to any subdivision of the
environment, and we leave a discussion of appropriate environment partitioning
to future work.

####  V-B2 Sample level

One can measure quality of guidance by considering both single and multiple
samples. In single-step evaluation we measure the expected quality of a single
node expansion of a single search tree. Such evaluation not only depends on
the task and environment, but on the current exploration progress, namely
information regarding previously sampled points in the configuration space
encoded in the search tree. In multiple-sample evaluation we instead measure
the holistic efficiency of all samples produced by the guiding space over the
course of planning.

In other words, single-step evaluation considers distributions over trees $T$,
whereas in multiple-sample evaluation considers distributions over
$\mathcal{C}$, implicitly described by the distribution of valid trees
produced by an algorithm. We refer the reader to [4] for a more in depth
discussion of multiple-sample evaluation, and here summarize that in that
setting our metric reduces to measuring the entropy of the observed sampling
distribution. In this case lower entropy means better guidance.

(a)

(b)

(c)

(d)

Figure 5: RRT, LazyPRM, MedialAxis, and PathDatabase (Lightning) evaluations.
All plots are averages over 128 random seeds and evaluation is done with a
smoothed target distribution that is at least $\epsilon=0.0001$ at each node
on the tree. When a seed finishes the search we consider its sampling
efficiency to be $0$ for all future samples (i.e., sampling the goal). Runtime
in PathDatabase only includes search time and not database creation, and the
guidance quality is averaged over the entire run instead of shown per
iteration as in the other algorithms.

##  VI Refactoring Existing Algorithms

Existing algorithms such as those cited in Section IV are not implemented in
the form of Algorithm 1, rather each uses their respective guiding space in a
different way. In the next section we conduct experiments on four guided
search algorithms, RRT [23], LazyPRM [10], DR-RRT [13], and Lightning [7].
Figure 4 shows examples of trees produced by these algorithms for solving the
tasks depicted in Figure 2, and Figure 3 visualizes their guiding spaces.

LazyPRM builds a probabilistic roadmap _without_ performing collision
detection with the environment (or only doing so for vertices but not edges).
Our implementation grows a tree using the length of the shortest path on this
lazy graph as guidance for selecting which node to expand on the tree. When an
invalid sample is found in c-space, we pass this information back to the
guiding space (the lazy graph) by deleting nearby vertices/edges.

DR-RRT uses the medial axis skeleton computed in workspace to guide the robot.
The original proposed implementation uses the edges of this skeleton to define
“dynamic regions” along which to sample, and does so by creating balls around
these regions and projecting back into c-space. For robots whose workspace
configuration is a subset of the dimensions of their c-space this is easy, but
for more complex c-space the general form of such a procedure involves an
inverse-kinematics projection. Our new implementation, which we call
MedialAxisGuidance is more straightforward - we project c-space nodes $v$ with
forward-kinematics into workspace, then project these workspace poses onto the
nearest medial axis edge $e_{v}$. The corresponding guidance value is then the
distance along the skeleton edges from the projected edge to the projection of
the goal, $d(e_{v},e_{t})$. Whenever a node $v$ fails to expand due to
obstacle collision we increase the weight of $e_{v}$, thereby reducing the
chance that the guidance encourages exploration along that edge.

Lightning uses a database of paths from related motion planning problems to
guide exploration. Namely, given a motion planning problem, a path that was
produced by the most similar task (start and goal) is selected from the
database, and then copied into the current c-space. Where the path is invalid
the algorithm fixes the solution by running a bi-directional RRT from the two
ends of the breakage. We highlight that this sub-problem can be just as
difficult as the original, and moreover that the full power of the database is
never used. Our new implementation, which we call PathDatabaseGuidance uses
the _length_ of the queried path as the guidance value. As with LazyPRM,
failed sample information is passed down to the guiding space by filtering out
paths that pass close to invalid regions. If no node can locally plan to any
path, then we default to RRT until a new node is found from which there is a
visible path. We note that this implementation is relatively similar to [1],
but we query the database using task similarity as in [7] rather than dynamic
time warping, and more importantly we actively filter the database using
information from the search tree.

##  VII Experiments

(a)

(b)

(c)

(d)

Figure 6: Results of the TrapCup environment experiments in which the
rectangular mobile robot must navigate from the center-left to the center-
right. MedialAxis struggles with the local minima of the first obstacle, which
resembles the Trap environment. LazyPRM struggles with the local minima of the
second obstacle, which resembles the Cup environment. A hybrid guiding space
which uses both projections can escape both minima faster and has overall
better guidance.

We design environments intended to demonstrate both pros and cons of each
guiding space mentioned in Section VI, using a rectangular mobile base
(3-dimensional c-space) in a 2-dimensional workspace.

We show results in Figure 5, plotting quality of guidance (lower is better)
over the course of search progress, i.e., number of samples. RRT initially
performs similarly on all environments, doing well and then poorly as it
struggles to get through narrow passages. LazyPRM performs well on
SimplePassage but poorly on Cup, where it spends time exploring the dead end.
MedialAxisGuidance does well on SimplePassage but poorly on Trap, where a
valid passage in workspace is not valid in c-space. These first three
experiments demonstrate that guidance matches our intuition for when each
algorithm performs well or struggles, which can be seen qualitatively in
Figure 4. In addition, we highlight that RRT is comparable to other algorithms
initially, and especially on environments where MedialAxisGuidance or LazyPRM
lead the search astray.

Next we test PathDatabaseGuidance on a dataset of RandomSimplePassages, an
environment with a randomly placed narrow passage between two rooms. We show
how guidance improves as dataset size increases, while runtime gets worse (as
querying the database takes longer), showing how such traditional metrics can
obscure properties of algorithms because they are implementation dependent.
This matches the intuition mentioned in Section III-C regarding the trade-off
between guidance quality and the difficulty of computing it.

###  VII-A Multi-guidance search

Finally, we create a hybrid environment TrapCup which combines features form
both the Trap and Cup environments. We now design a hybrid guiding space that
is a combination of LazyPRM and MedialAxis. In other words, it creates both
guiding spaces (i.e., the cross product) and uses a rule to determine how to
select or combine the heuristics from the two. In the experiments shown in
Figure 6 we simply swap the source of the guidance each time the search
collides with an obstacle. Intuitively, each form of guidance faces certain
local minima and by using multiple guiding spaces we can quickly escape those
minima by passing control to another guiding space. Notice that the individual
guidance plots for each guiding space indicate that these local minima are
distinct (e.g., the curves have different peaks), and thus we see the hybrid
guiding space out-performs its parts. Moreover, using number of samples alone
as a metric would obscure this fact since MedialAxis and LazyPRM take about
the same number of samples to find the goal. These results parallel what one
would expect from running a multi-heuristic A* search [2] as opposed to single
heuristic search, and highlights the value of framing sampling based motion
planning as guided search; we now have a structured way of combining different
algorithmic ideas from the literature as multiple guiding spaces, whereas the
original implementation are not so conducive to modularization.

##  VIII Conclusion

In this work we presented the guiding space framework, in which an explicit
guiding space provides heuristic guidance to a tree based search algorithm in
a configuration space. We showed how many existing methods fit this framework,
albeit implicitly. By viewing such methods as guiding spaces, we can more
easily compare and combine them into hybrid algorithms. As such, we defined
_sampling efficiency_ as a way to measure the quality of biased node
selections of a search tree, and showed how this metric matches intuition for
a variety of algorithms and environments.

## References

  * [1] Sandip Aine, Charupriya Sharma, and Maxim Likhachev.  Learning to search more efficiently from experience: A multi-heuristic approach.  In Proceedings of the International Symposium on Combinatorial Search, volume 6, pages 141–145, 2015.
  * [2] Sandip Aine, Siddharth Swaminathan, Venkatraman Narayanan, Victor Hwang, and Maxim Likhachev.  Multi-heuristic a.  The International Journal of Robotics Research, 35(1-3):224–243, 2016.
  * [3] Nancy M Amato, O Burchan Bayazit, Lucia K Dale, Christopher Jones, Daniel Vallejo, et al.  Obprm: An obstacle-based prm for 3d workspaces.  In Proc. Int. Workshop on Algorithmic Foundations of Robotics (WAFR), pages 155–168, 1998.
  * [4] Amnon Attali, Stav Ashur, Isaac Burton Love, Courtney McBeth, James Motes, Diane Uwacu, Marco Morales, and Nancy M Amato.  Evaluating guiding spaces for motion planning.  arXiv preprint arXiv:2210.08640, 2022.
  * [5] Amnon Attali, Marco Morales, and Nancy M Amato.  Discrete state-action abstraction via the successor representation.  arXiv preprint arXiv:2206.03467, 2022.
  * [6] O. Burçhan Bayazit, Dawen Xie, and Nancy M. Amato.  Iterative relaxation of constraints: a framework for improving automated motion planning.  In 2005 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2005, pages 3433–3440. IEEE, 2005.
  * [7] Dmitry Berenson, Pieter Abbeel, and Ken Goldberg.  A robot path planning framework that learns from experience.  In 2012 IEEE International Conference on Robotics and Automation, pages 3671–3678. IEEE, 2012.
  * [8] Mohak Bhardwaj, Sanjiban Choudhury, and Sebastian Scherer.  Learning heuristic search via imitation.  In Conference on Robot Learning, pages 271–280. PMLR, 2017.
  * [9] Botond Bocsi, Lehel Csató, and Jan Peters.  Alignment-based transfer learning for robot models.  In The 2013 international joint conference on neural networks (IJCNN), pages 1–7. IEEE, 2013.
  * [10] Robert Bohlin and Lydia E. Kavraki.  Path planning using lazy PRM.  In International Conference on Robotics and Automation, ICRA, pages 521–528. IEEE, 2000.
  * [11] Constantinos Chamzas, Zachary Kingston, Carlos Quintero-Peña, Anshumali Shrivastava, and Lydia E Kavraki.  Learning sampling distributions using local 3d workspace decompositions for motion planning in high dimensions.  In 2021 IEEE International Conference on Robotics and Automation (ICRA), pages 1283–1289. IEEE, 2021.
  * [12] David Coleman, Ioan A Şucan, Mark Moll, Kei Okada, and Nikolaus Correll.  Experience-based planning with sparse roadmap spanners.  In 2015 IEEE International Conference on Robotics and Automation (ICRA), pages 900–905. IEEE, 2015.
  * [13] Jory Denny, Read Sandström, Andrew Bregger, and Nancy M. Amato.  Dynamic region-biased rapidly-exploring random trees.  In Proceedings of the Twelfth Workshop on the Algorithmic Foundations of Robotics, WAFR 2016,, volume 13 of Springer Proceedings in Advanced Robotics, pages 640–655. Springer, 2016.
  * [14] Coline Devin, Abhishek Gupta, Trevor Darrell, Pieter Abbeel, and Sergey Levine.  Learning modular neural network policies for multi-task and multi-robot transfer.  In 2017 IEEE international conference on robotics and automation (ICRA), pages 2169–2176. IEEE, 2017.
  * [15] Aleksandra Faust, Kenneth Oslund, Oscar Ramirez, Anthony Francis, Lydia Tapia, Marek Fiser, and James Davidson.  Prm-rl: Long-range robotic navigation tasks by combining reinforcement learning and sampling-based planning.  In 2018 IEEE international conference on robotics and automation (ICRA), pages 5113–5120. IEEE, 2018.
  * [16] Roland Geraerts and Mark H Overmars.  Creating high-quality paths for motion planning.  The international journal of robotics research, 26(8):845–863, 2007\.
  * [17] Mukulika Ghosh, Shawna L. Thomas, Marco Morales, Samuel Rodríguez, and Nancy M. Amato.  Motion planning using hierarchical aggregation of workspace obstacles.  In 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 5716–5721. IEEE, 2016.
  * [18] Kris Hauser.  Lazy collision checking in asymptotically-optimal motion planning.  In IEEE International Conference on Robotics and Automation, ICRA 2015, Seattle, WA, USA, 26-30 May, 2015, pages 2951–2957. IEEE, 2015\.
  * [19] Christopher Holleman and Lydia E Kavraki.  A framework for using the workspace medial axis in prm planners.  In IEEE International Conference on Robotics and Automation, volume 2, pages 1408–1413. IEEE, 2000.
  * [20] Brian Ichter and Marco Pavone.  Robot motion planning in learned latent spaces.  IEEE Robotics and Automation Letters, 4(3):2407–2414, 2019.
  * [21] Léonard Jaillet and Thierry Siméon.  A prm-based motion planner for dynamically changing environments.  In 2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)(IEEE Cat. No. 04CH37566), volume 2, pages 1606–1611. IEEE, 2004.
  * [22] Lydia E. Kavraki, Petr Svestka, Jean-Claude Latombe, and Mark H. Overmars.  Probabilistic roadmaps for path planning in high-dimensional configuration spaces.  IEEE Trans. Robotics Autom., 12(4):566–580, 1996.
  * [23] Steven LaValle.  Rapidly-exploring random trees: A new tool for path planning.  Research Report 9811, 1998.
  * [24] Tomas Lozano-Perez.  A simple motion-planning algorithm for general robot manipulators.  IEEE Journal on Robotics and Automation, 3(3):224–238, 1987.
  * [25] Arthur Mahoney, Joshua Bross, and David Johnson.  Deformable robot motion planning in a reduced-dimension configuration space.  In IEEE International Conference on Robotics and Automation, ICRA 2010, pages 5133–5138. IEEE, 2010.
  * [26] Erion Plaku, Lydia E. Kavraki, and Moshe Y. Vardi.  Discrete search leading continuous exploration for kinodynamic motion planning.  In Robotics: Science and Systems III, 2007. The MIT Press, 2007\.
  * [27] Ahmed H Qureshi, Anthony Simeonov, Mayur J Bency, and Michael C Yip.  Motion planning networks.  In 2019 International Conference on Robotics and Automation (ICRA), pages 2118–2124. IEEE, 2019.
  * [28] Kensen Shi, Jory Denny, and Nancy M Amato.  Spark prm: Using rrts within prms to efficiently explore narrow passages.  In 2014 IEEE International Conference on Robotics and Automation (ICRA), pages 4659–4666. IEEE, 2014.
  * [29] Vojtech Vonásek, Robert Penicka, and Barbora Kozlíková.  Computing multiple guiding paths for sampling-based motion planning.  In 19th International Conference on Advanced Robotics, ICAR 2019, pages 374–381. IEEE, 2019.
  * [30] Jiankun Wang, Wenzheng Chi, Chenming Li, Chaoqun Wang, and Max Q-H Meng.  Neural rrt*: Learning-based optimal path planning.  IEEE Transactions on Automation Science and Engineering, 17(4):1748–1758, 2020.