# Measurement of enhanced spin-orbit coupling strength for donor-bound
electron spins in silicon

Radha Krishnan  Beng Yee Gan  School of Physical and Mathematical Sciences,
Nanyang Technological University, 21 Nanyang Link, Singapore 637371  Yu-Ling
Hsueh  School of Physics, University of New South Wales, Sydney, NSW 2052,
Australia  A.M. Saffat-Ee Huq  School of Physics, University of New South
Wales, Sydney, NSW 2052, Australia  Jonathan Kenny  School of Physical and
Mathematical Sciences,
Nanyang Technological University, 21 Nanyang Link, Singapore 637371  Rajib
Rahman  School of Physics, University of New South Wales, Sydney, NSW 2052,
Australia  Teck Seng Koh  School of Physical and Mathematical Sciences,
Nanyang Technological University, 21 Nanyang Link, Singapore 637371  Michelle
Y. Simmons  Centre for Quantum Computation and Communication Technology,
School of Physics, University of New South Wales, Sydney, NSW 2052, Australia
Bent Weber<EMAIL_ADDRESS>School of Physical and Mathematical Sciences,
Nanyang Technological University, 21 Nanyang Link, Singapore 637371

††preprint: APS/123-QED

##  1 Abstract

While traditionally considered a deleterious effect in quantum dot spin
qubits, the spin-orbit interaction is recently being revisited as it allows
for rapid coherent control by on-chip AC electric fields. For electrons in
bulk silicon, SOC is intrinsically weak, however, it can be enhanced at
surfaces and interfaces, or through atomic placement. Here we show that the
strength of the spin-orbit coupling can be locally enhanced by more than two
orders of magnitude in the manybody wave functions of multi-donor quantum dots
compared to a single donor, reaching strengths so far only reported for holes
or two-donor system with certain symmetry. Our findings may provide a pathway
towards all-electrical control of donor-bound spins in silicon using electric
dipole spin resonance (EDSR).

##  2 Introduction

Donor bound electron spins [1, 2] in silicon are promising qubit candidates
for the implementation of quantum information processing due to their
exceptionally long spin relaxation ($T_{1}$) [3] and coherence times ($T_{2}$)
[4]. This is owing to a low concentration of nuclear spins in the crystalline
bulk matrix [4, 5], especially in isotopically enriched 28Si, and generally
weak coupling to charge noise due to the weak coupling of the electron’s spin
to orbital degrees of freedom [6]. For donor spin qubits that are embedded
deep within the silicon bulk crystal, far from any interface, inversion
asymmetry of the silicon lattice precludes Dresselhaus spin-orbit coupling
(SOC), and an isotropic Coulombic confinement potential gives rise to
approximately spheroid-symmetric envelope wave functions in the lowest valley-
orbit states, allowing for only weak Rashba SOC and higher-order electro-
magnetic SOC [7].

Traditionally, SOC has been considered as an undesired effect in spin qubits,
as it can dominate spin relaxation [8, 9, 10], lead to state leakage through
mixing of singlet ($S$) and triplet ($T^{-}$) states [11], or an overall
increased sensitivity to charge noise [12]. In turn, however, SOC may be
harnessed for rapid electrical control of individual spin qubits [13]. Much
effort has therefore recently been directed towards driving electric-dipole
spin resonance (EDSR) transitions [14, 15] via the Stark shift [16] or by
utilizing the disparity in the $g$-factors [17] for individual addressability
of qubits via a microwave field [18, 19, 20, 21, 22]. Demonstrations have so
far largely been restricted to hole spin qubits bound to silicon quantum dots
[23, 24] and acceptor atoms [25, 26] in which SOC-mediated enhanced coupling
to external AC electric fields [27, 28] may be combined with long spin life-
[29] and coherence times [30]. Enhanced SOC, however, has more recently also
been reported for electron spin qubits in metal-oxide-semiconductor (MOS)
quantum dots [16, 19, 31, 32, 24], as well as in donors [7, 33]. However, the
strength of the SOC ($\sim 100$ neV) reported [16, 19, 31, 32, 24] is usually
2-3 orders of magnitude weaker for electrons compared to holes [24, 26] (see
Table 1).

Several factors have been reported to enhance the weak SOC for electrons in
silicon [34, 12, 7]. For instance, electric fields at surfaces or interfaces,
such as Si/SiO2 [35] or SiGe/Si [34], can enhance Rashba SOC [36], while
breaking of mirror symmetry and atomic-scale roughness of interfaces has been
shown to produce a significant Dresselhaus SOC contribution [12]. In donor-
based systems, strong electric fields can displace the electron wave function
from the donor nucleus and can give rise to electric field-induced SOC [7].
Also, engineering two donor atoms at specific locations can change local
symmetry and result in both Rashba and Dresselhaus SOC [33].For hole spins,
SOC is known to be strongly enhanced, owing to the reduced symmetry of the
non-spheroid spin carrier wave function [37].

Figure 1: Si:P DQD Device. (a) Overview STM-image after hydrogen resist
lithography, showing two quantum dots ($D_{1,2}$), connected to source (S),
drain (D), and gate ($G_{1,2}$) electrodes. (b) The charge stability diagram
of the DQD shows a number of closed honeycomb domains, where $(M,N)$ indicate
the effective electron occupation. (c,d) Atomic-resolution STM images of the
two quantum dots allowing to estimate the approximate number of incorporated P
donors (e,f). (g-i) Self-consistent confinement potential (i) and tight-
binding wave-functions, calculated for the $N=15$ (g) and $N=10$ (h) charge
state in NEMO-3D.

Similarly, as we show in this work, SOC can be significantly enhanced by more
than 2 orders of magnitude for electrons in multi-donor quantum dots due to an
interplay of the atomically-sharp confinement potential with dopant placement
disorder, a small valley-orbit splitting and non-spherical manybody wave
functions. We draw our conclusions from the transport spectroscopy of a multi-
donor double quantum dot in the Pauli spin blockade regime, in which a
pronounced leakage is observed at the degeneracy point of the effective
valence spin singlet $S(0,2)$ and triplet $T^{-}(1,1)$ states. A transport
model of spin-dependent second-order tunneling confirms that this resonance
can only be explained by an SOC mediated coupling of the singlet and triplet
states, while atomistic tight-binding calculations [38] confirm the magnitude
of the SOC observed.

Table 1: Spin-orbit interaction strength for electrons and holes in silicon. Reference | Carrier | SOC strength
---|---|---
Harvey-Collard et al. [31] | Electron | (113$\pm$22) neV
Tanttu et al. [32] | Electron |  $\lesssim$ 113.73 neV
Jock et al. [35] | Electron |  $\lesssim$ 82.7 neV
Weber et al. [7] | Electron |  $\lesssim$ 1.05 $\mu$eV
Hsueh et al. [33] | Electron | 33.9$\mu$eV
Li et al. [24] | Hole | 110 $\mu$eV
Heijden et al. [26] | Hole | (36$\pm$5) $\mu$eV
This Work | Electron | (40$\pm$10) $\mu$eV

##  3 Results and Discussion

The multi-donor double quantum dot device is shown in Fig.1a, fabricated by
STM hydrogen resist lithography [39]. Both quantum dots $D_{1}$ and $D_{2}$
are $\sim$4 nm in diameter, separated by $17.5\pm 0.5$ nm (centre to centre),
each containing approximately $N_{D}\simeq 15$ phosphorus (P) donors [39] as
estimated from the size of the lithographically defined dots (Fig. 1c-f),
assuming that six contiguous dimers of the Si(001)-2$\times$1 surface are
required to incorporate a phosphorus atom into the silicon matrix [40]. Both
dots are mutually tunnel-coupled and are coupled to source (S) and drain (D)
electrodes which have been staggered to allow independent electrostatic
control of the individual dots’ electrochemical potentials by gate electrodes
G1 and G2 [39]. A DC magnetic field was applied within the plane of the
device, along the Si$\left<110\right>$ direction, i.e. $\sim 30\degree$ from
the axis connecting the two quantum dots (arrow in Fig. 1a).

Figure 1b shows the measured charge stability diagram, indicating electron
numbers $N$ and $M$ per dot. Absolute electron numbers cannot be determined as
the dots could not be fully depleted within the accessible gate-voltage range.
However, we expect $N,M\approx N_{\rm D}\simeq 15$ near charge neutrality
($V_{\rm G1}=V_{\rm G2}=0$ V) [39, 41].

Figure 2: Low-bias magnetospectroscopy of the Si:P DQD Device. (a-d) Close-up
of the charge transition highlighted by the dashed circle in Figs. 1b,
measured at $V_{D}=\pm 500\mu V$, and at $B=0$ T (a,b) and $B=6.1$ T (c,d),
respectively. (e,f) Co-tunneling current, measured along the dashed arrow in
(b) and plotted as a function of magnetic field and detuning, compared to the
corresponding energy eigenspectrum (e) of effective two-electron singlet and
triplet states. $t_{\rm c}$ denotes the tunnel coupling, and $t_{\rm SOC}$ the
spin-orbit coupling. The bias window is indicated in (f) and the data inverted
in $B$ for ease of comparison with (e).

Despite the large electron numbers, we observe well-defined effective two-
electron singlet and triplet states as summarized in Fig. 2. A close-up of the
charge transition highlighted by the black dashed circle in Fig. 1b is shown
in Figure 2a-d, recorded at a low bias of $V_{SD}=\pm 500~{}\mu$V and at two
different values of the applied in-plane magnetic field ($B=0$ T and $B=6.1$
T). At larger bias, the charge degeneracy (triple) points widen into bias
triangles [39], exhibiting conventional Pauli spin blockade (see Supporting
Information). Current along the outlines of the charge stable domains at
$B=100$ mT (dashed white lines) arises from higher-order (co-)tunneling.

Current resonances connecting the triple points reflect co-tunneling due to a
strong mutual tunnel coupling between the two quantum dots. At finite magnetic
fields $B=6.1$ T (Fig. 2c,d), we observe that the single inter-dot co-
tunneling resonance splits into two. This is more clearly observed in the
evolution of the co-tunneling current resonances as a function of magnetic
field, recorded along the dashed arrow in Fig. 2b, and plotted in Fig. 2f. The
corresponding energy eigenspectrum of SOC-perturbed effective $(1,1)$ and
$(0,2)$ singlet- ($S$) and triplet- ($T$) valence spin states [42] is shown in
Fig. 2e. For ease of comparison, the data was rescaled with respect to
detuning energy $\varepsilon$, assuming a lever arm ($\alpha=0.08\pm 0.01$),
extracted from the observed bias window ($500~{}\mu$eV), and with the magnetic
field axis inverted. We define zero detuning ($\varepsilon=0$) where the
singlet states $S(1,1)$ and $S(0,2)$ are degenerate, and hybridize due to the
tunnel coupling $t_{c}$. This is visible from a high intensity resonance at
low magnetic fields ($B<1$ T), that is symmetric upon bias reversal.

At finite magnetic field, the degeneracy of the $T(1,1)$ and $T(0,2)$ triplet
manifold is lifted, resulting in three non-degenerate levels each, split by
the Zeeman energy, $\Delta E_{Z}=S_{z}g\mu_{B}B$. Positive detuning lowers the
energy of the $(0,2)$ states, allowing the triplet states to hybridize at
$\varepsilon=1.21$ meV. This resonance becomes only visible at $B\sim 5$ T,
when the hybridized $T^{-}$ states enter the bias window (horizontal arrows in
Fig. 2f). As the detuning at which level alignment of $S(1,1)-S(0,2)$ and
$T(1,1)-T(0,2)$ states occurs does not depend on magnetic field, the
corresponding transport resonances do not change position. The finite magnetic
field also leads to a change in the $(1,1)$ ground state from singlet to
triplet, giving rise to a level crossing of $T^{-}(1,1)$ and the $S(0,2)$
states when $\Delta E_{Z}>t_{c}$ at $B\simeq 1$ T. An enhanced occupation
probability of the $T^{-}(1,1)$ over the $S(1,1)$ result in a suppression of
co-tunneling current into the drain lead, and can be understood as a
manifestation of Pauli spin blockade in the co-tunneling regime [43]. Current
at $B\lesssim 1$ T, also reflected in our transport calculations (Fig. 3),
likely has its origin in spin-flip co-tunneling with the S lead [44, 45].

Figure 3: Effects of SOC in second order transport. (a) Second order transport
model of the data in Fig. 2, with and without SOC-mediated coherent coupling
included (see main text and Supporting Information for detail). No significant
co-tunneling current is observed when SOC-mediated coupling is absent while a
clear $T^{-}(1,1)-S(0,2)$ transport resonance is observed in the presence of
SOC. (b) Line cuts comparing experiment (Fig. 2f) and theory in terms of
resonance line widths and intensities of for low, intermediate and high
fields, as indicated by the red arrows in (a). Red and green lines show our
results at $C_{\text{SOC}}=40~{}\mu\mathrm{eV}$ and $0~{}\mu\mathrm{eV}$,
respectively, and with Gaussian (charge) noise of magnitude
$\sigma_{\varepsilon}=17~{}\mu$eV [46]) in the detuning parameter. We see that
the charge noise has the effect of broadening and reducing the intensity of
the $S(0,2)-T^{-}(1,1)$ resonance, as can be seen by comparing to results
without charge noise (black dashed line).

At fields $B\gtrsim 1$ T, spin selection rules should prohibit co-tunneling at
the $T^{-}(1,1)-S(0,2)$ level crossing. Yet, a fainter but sharp and clear
current resonance is observed that shifts linearly as a function of magnetic
field. Such current resonance can only be explained by the presence of a spin
non-conserving coupling, such as inelastic co-tunneling [43, 47], nuclear
hyperfine interaction [48, 49, 50, 51, 52], or, spin orbit coupling [8, 31,
53, 54, 9].

Such coupling of $T^{-}(1,1)-S(0,2)$ states as observed can be understood from
the formation of spin-orbit perturbed valley-orbit eigenstates [55, 19]. From
first-order perturbation theory, the valley-states in each QD can be expressed
as admixtures of both spin polarities by virtue of the spin-orbit coupling,

$\displaystyle\ket{\Downarrow}=\ket{v_{i},\downarrow}-\frac{C_{\rm
SOC}}{\Delta E+\tfrac{1}{2}g\mu_{B}B}\ket{v_{j},\uparrow}+...$ (1)

and

$\displaystyle\ket{\Uparrow}=\ket{v_{i},\uparrow}-\frac{C^{*}_{\rm
SOC}}{\Delta E-\tfrac{1}{2}g\mu_{B}B}\ket{v_{j},\downarrow}+...$ (2)

Here, $E_{Z}=g\mu_{B}B$ is the Zeeman energy with the electron $g$-factor
$g\approx 2$, and $\Delta E$ is the valley-orbit splitting. $C_{\rm
SOC}=\left\langle v_{j},\uparrow\middle|H_{\rm
SOC}\middle|v_{i},\downarrow\right\rangle$ is the SOC matrix element. The
degree of admixture of the spin-valley states thus sensitively depends on the
magnitude of the valley-obit splitting $\Delta E$. Naturally assumed to be
large in atomically confined systems, the valley multiplicity in silicon,
especially at large charge occupation [56, 39], provides for comparatively
small valley-orbit gaps in the range of a few hundred $\mu$eV [56] to a few
meV [39], and thus comparable in magnitude to the much larger few-electron
quantum dots realized in SiGe or Si MOS [55, 19].

Given that the spin-valley states of the individual quantum dots contain
admixtures of both the spin polarities, they further give rise to a coherent
coupling between singlet and triplet states across the two quantum dots of
magnitude $C_{\text{SOC}}/\sqrt{2}$. At intermediate magnetic fields, where
these are ground states, the magnitude of the effective coupling can be
estimated from a Schrieffer-Wolff transformation, as

$t_{\text{SOC}}=t_{c}C_{\text{SOC}}/\sqrt{2}\left|\tfrac{1}{\varepsilon-g\mu_{B}B}-\tfrac{1}{\varepsilon-\Delta
E}\right|$ (3)

This is hence a second order coupling in which $S(0,2)$ is coupled to
$T^{-}(0,2)$ via SOC in the occupied dot, and $T^{-}(0,2)$ to $T^{-}(1,1)$ via
tunneling.

We confirm this notion from a transport model in the co-tunneling regime as
summarized in Fig. 3 (see Supporting Information for details). Our model
reproduces all essential features of the data, but most notably the presence
of a single pronounced co-tunneling resonance at intermediate magnetic fields,
observed only when the SOC-mediated coupling between the dots is present. This
is clearly reflected in the colorplots of Fig. 3a as well as the individual
line traces (Fig. 3b) at three different strengths of the magnetic field,
comparing the measured drain current with that calculated in second-order
transport.

We note that the overall current intensity is dominated by the product of dot-
lead tunnel rates $\Gamma_{S}$ and $\Gamma_{D}$, while the relative
intensities of the respective resonances depend non-linearly on these rates.
For example, the current intensity at the singlet resonance observed at low
magnetic field ($<1$ T) is strongly dependent on dot-lead tunnel rates as it
is determined through steady-state DQD populations. Spin blockade occurs when
the DQD is in a triplet (1,1) state, lifted by co-tunneling of a spin between
a dot and the adjacent lead. Such spin non-preserving processes, in particular
spin-flip co-tunneling [44], are therefore implicitly taken into account in
the model, but can only explain strong resonant current at $B<1$ T. For $B>1$
T, the onset of Pauli spin blockade sets an upper bound on $t_{c}$.

Our model also confirms that the resonance lineshapes are Lorentzian when
lifetime broadening is dominant. In this case, linewidths and peak intensities
are determined by the magnitude of $t_{c}$ and $C_{\text{SOC}}$.
Interestingly, however, the $S(0,2)-T^{-}(1,1)$ resonance also possesses a
Gaussian character which suggests an inhomogeneous broadening mechanism, for
instance due to charge noise in detuning. We can model this by convolving the
calculated current (see Supporting Information) with a Gaussian of half-width
$\sigma_{\varepsilon}=17~{}\mu$eV, in reasonable agreement with recent
measurements of bandwidth-integrated charge noise in epitaxial Si:P devices
[46, 2]. We note that the current lineshape through the $T^{-}(1,1)-S(0,2)$
channel is under-determined if both $t_{c}$ and $C_{\text{SOC}}$ are
unbounded. However, the upper bound on $t_{c}$ constrains $C_{\text{SOC}}$
from below, which can then be estimated. In a similar way, the upper bound on
$t_{c}$ allows the energy splitting $\Delta E$ to be estimated from the
lineshape and position of the triplet channel.

Figure 4: Calculated SOC strength. Strength of the spin-orbit coupling
$C_{\textrm{SOC}}$ as calculated from atomistic tight-binding, and plotted as
function of the electron number $N$. The decrease of the SOC strength towards
higher electron numbers likely reflects screening of the donor-based Coulombic
confinement potential. Insert shows the envelope wavefunction for N =15.

From fits to the data, we extract $C_{\text{SOC}}\approx 40~{}\mu\mathrm{eV}$,
more than two orders of magnitude larger than values previously reported for
electrons in silicon, reaching values previously only observed for holes
(compare Table 1). Such considerable enhancement of the SOC strength likely
arises as a combination of several factors, including the atomically-abrupt
and disordered Coulombic confinement potential of the dots (Fig. 1h) lowering
wave function symmetry and hence increasing the symmetry dependent SOC
transition dipole moment [57]. Indeed, we can confirm the magnitude of the SOC
strength by atomistic tight-binding calculations of the multi-donor dots, as
shown in Fig. 4, plotting $C_{\textrm{soc}}$ against electron occupation $N$.
Sixteen P donors, represented by Coulomb charges each with a cutoff potential
capturing the central-cell correction, are placed in the Si crystal within a 4
by 4 nm region. The tight-binding Hamiltonian is set up using 20 atomic
orbital sp3d5s* basis including nearest-neighbor interactions [58]. A self-
consistent Hartree method is used to capture the multi-electron wavefunction,
which has successfully reproduced the two-electron binding energy in a single
donor [59]. The SOC estimates have been extracted from the self-consistent TB
wave functions, directly, by considering calculations with spin-orbit
($\ket{\Uparrow}$, $\ket{\Downarrow}$) and without ($\ket{v_{i},\uparrow}$,
$\ket{v_{j},\downarrow}$) by setting the silicon tight-binding spin-orbit
parameter to the bulk value [60, 61] or zero, respectively. A projection from
the spin-orbit wave function onto the pure orbital basis then gives us $C_{\rm
SOC}$ as plotted in Fig. 4, reaching values between 11 $\mu$eV ($N=16$) and 90
$\mu$eV ($N=6$), matching the value extracted from the experiments (40
$\mu$eV) in their order of magnitude. For a range of different disordered
donor configurations within the dots, we find values $\sim$10-20 $\mu$eV at
$N\simeq 15$, confirming that these results are relative robust against dopant
placement disorder within the dots. While the model slightly underestimates
the SOC strength at $N\simeq 15$, a strong enhancement towards lower electron
occupation likely reflects the less effective screening of the Coulombic
confinement potential, amplifying the effects of donor disorder and wave
function asymmetry.

##  4 Conclusion

To conclude, we have reported a strong enhancement of the spin-orbit coupling
(SOC) strength for spins in the manybody wavefunctions of multi-donor silicon
quantum dots, with an SOC energy scale as much as $40~{}\mu$eV – exceeding
previous reports for electrons in silicon quantum dots by two orders of
magnitude. We explain this enhancement by the abrupt and disordered nature of
the Coulombic confinement potential inherent to donor-based quantum dots,
reducing wave-function symmetry. Such strongly enhanced SOC in silicon –
similar in magnitude to that for holes – may provide a pathway towards all-
electrical control of donor-bound spins in silicon by electric dipole spin
resonance (EDSR) as well as may enhance coupling of donor spins to
superconducting resonators.

## Acknowledgements

This research is supported under the Singapore Quantum Engineering Programme
(QEP2.0) “Atomic Engineering of Donor-based Spin Qubits in Silicon”
(NRF2021-QEP2-02-P07). BW acknowledges a Singapore National Research
Foundation (NRF) Fellowship (NRF-NRFF2017-11).

## References

  * [1] T. F. Watson, S. G. J. Philips, K. E., D. R. Ward, S. P., V. M., D. E. Savage, M. G. Lagally, M. Friesen, M. A. Coppersmith, S. N. Eriksson, and V. L. M. K., “A programmable two-qubit quantum processor in silicon,” Nature, vol. 555, pp. 633––637, 2018.
  * [2] Y. He, S. K. Gorman, D. Keith, L. Kranz, J. G. Keizer, and M. Y. Simmons, “A two-qubit gate between phosphorus donor electrons in silicon,” Nature, vol. 571, no. 7765, p. 371, 2019.
  * [3] T. F. Watson, B. Weber, Y.-L. Hsueh, L. C. Hollenberg, R. Rahman, and M. Y. Simmons, “Atomically engineered electron spin lifetimes of 30 s in silicon,” Science Advances, vol. 3, no. 3, p. e1602811, 2017.
  * [4] W. M. Witzel, M. S. Carroll, A. Morello, Ł. Cywiński, and S. D. Sarma, “Electron spin decoherence in isotope-enriched silicon,” Physical Review Letters, vol. 105, no. 18, p. 187602, 2010.
  * [5] M. Veldhorst, J. Hwang, C. Yang, A. Leenstra, B. de Ronde, J. Dehollain, J. Muhonen, F. Hudson, K. M. Itoh, A. Morello, et al., “An addressable quantum dot qubit with fault-tolerant control-fidelity,” Nature Nanotechnology, vol. 9, no. 12, p. 981, 2014.
  * [6] F. A. Zwanenburg, A. S. Dzurak, A. Morello, M. Y. Simmons, L. C. Hollenberg, G. Klimeck, S. Rogge, S. N. Coppersmith, and M. A. Eriksson, “Silicon quantum electronics,” Reviews of Modern Physics, vol. 85, no. 3, p. 961, 2013.
  * [7] B. Weber, Y.-L. Hsueh, T. F. Watson, R. Li, A. R. Hamilton, L. C. Hollenberg, R. Rahman, and M. Y. Simmons, “Spin–orbit coupling in silicon for electrons bound to donors,” npj Quantum Information, vol. 4, no. 1, pp. 1–5, 2018\.
  * [8] Y. Song and S. D. Sarma, “Spin-orbit coupling induced two-electron relaxation in silicon donor pairs,” Physical Review B, vol. 96, no. 11, p. 115444, 2017.
  * [9] P. Scarlino, E. Kawakami, P. Stano, M. Shafiei, C. Reichl, W. Wegscheider, and L. Vandersypen, “Spin-relaxation anisotropy in a gaas quantum dot,” Physical Review Letters, vol. 113, no. 25, p. 256802, 2014.
  * [10] J. I. Climente, C. Segarra, and J. Planelles, “Spin–orbit-induced hole spin relaxation in inas and gaas quantum dots,” New Journal of Physics, vol. 15, no. 9, p. 093009, 2013.
  * [11] D. Stepanenko, M. Rudner, B. I. Halperin, and D. Loss, “Singlet-triplet splitting in double quantum dots due to spin-orbit and hyperfine interactions,” Physical Review B, vol. 85, no. 7, p. 075416, 2012.
  * [12] R. Ferdous, K. W. Chan, M. Veldhorst, J. Hwang, C. Yang, H. Sahasrabudhe, G. Klimeck, A. Morello, A. S. Dzurak, and R. Rahman, “Interface-induced spin-orbit interaction in silicon quantum dots and prospects for scalability,” Physical Review B, vol. 97, no. 24, p. 241401, 2018.
  * [13] S. Nadj-Perge, S. Frolov, E. Bakkers, and L. P. Kouwenhoven, “Spin–orbit qubit in a semiconductor nanowire,” Nature, vol. 468, no. 7327, pp. 1084–1087, 2010.
  * [14] M. Shafiei, K. Nowack, C. Reichl, W. Wegscheider, and L. Vandersypen, “Resolving spin-orbit-and hyperfine-mediated electric dipole spin resonance in a quantum dot,” Physical Review Letters, vol. 110, no. 10, p. 107601, 2013.
  * [15] E. I. Rashba, “Theory of electric dipole spin resonance in quantum dots: Mean field theory with gaussian fluctuations and beyond,” Physical Review B, vol. 78, no. 19, p. 195302, 2008.
  * [16] R. Ferdous, E. Kawakami, P. Scarlino, M. P. Nowak, D. Ward, D. Savage, M. Lagally, S. Coppersmith, M. Friesen, M. A. Eriksson, et al., “Valley dependent anisotropic spin splitting in silicon quantum dots,” npj Quantum Information, vol. 4, no. 1, pp. 1–8, 2018.
  * [17] R. Ruskov, M. Veldhorst, A. S. Dzurak, and C. Tahan, “Electron g-factor of valley states in realistic silicon quantum dots,” Physical Review B, vol. 98, no. 24, p. 245424, 2018.
  * [18] M. Veldhorst, R. Ruskov, C. Yang, J. Hwang, F. Hudson, M. Flatté, C. Tahan, K. M. Itoh, A. Morello, and A. Dzurak, “Spin-orbit coupling and operation of multivalley spin qubits,” Physical Review B, vol. 92, no. 20, p. 201401, 2015.
  * [19] A. Corna, L. Bourdet, R. Maurand, A. Crippa, D. Kotekar-Patil, H. Bohuslavskyi, R. Laviéville, L. Hutin, S. Barraud, X. Jehl, et al., “Electrically driven electron spin resonance mediated by spin–valley–orbit coupling in a silicon quantum dot,” npj Quantum Information, vol. 4, no. 1, pp. 1–7, 2018.
  * [20] W. Huang, M. Veldhorst, N. M. Zimmerman, A. S. Dzurak, and D. Culcer, “Electrically driven spin qubit based on valley mixing,” Physical Review B, vol. 95, no. 7, p. 075403, 2017.
  * [21] W. Huang, C. Yang, K. Chan, T. Tanttu, B. Hensen, R. Leon, M. Fogarty, J. Hwang, F. Hudson, K. M. Itoh, et al., “Fidelity benchmarks for two-qubit gates in silicon,” Nature, vol. 569, no. 7757, pp. 532–536, 2019\.
  * [22] F. Krauth, S. Gorman, Y. He, M. Jones, P. Macha, S. Kocsis, C. Chua, B. Voisin, S. Rogge, R. Rahman, et al., “Flopping-mode electric dipole spin resonance in phosphorus donor qubits in silicon,” Physical Review Applied, vol. 17, no. 5, p. 054006, 2022.
  * [23] D. Q. Wang, O. Klochan, J.-T. Hung, D. Culcer, I. Farrer, D. A. Ritchie, and A. R. Hamilton, “Anisotropic pauli spin blockade of holes in a gaas double quantum dot,” Nano Letters, vol. 16, no. 12, pp. 7685–7689, 2016.
  * [24] R. Li, F. E. Hudson, A. S. Dzurak, and A. R. Hamilton, “Pauli spin blockade of heavy holes in a silicon double quantum dot,” Nano Letters, vol. 15, no. 11, pp. 7314–7318, 2015.
  * [25] R. Maurand, X. Jehl, D. Kotekar-Patil, A. Corna, H. Bohuslavskyi, R. Laviéville, L. Hutin, S. Barraud, M. Vinet, M. Sanquer, et al., “A cmos silicon spin qubit,” Nature Communications, vol. 7, no. 1, pp. 1–6, 2016.
  * [26] J. van der Heijden, T. Kobayashi, M. G. House, J. Salfi, S. Barraud, R. Laviéville, M. Y. Simmons, and S. Rogge, “Readout and control of the spin-orbit states of two coupled acceptor atoms in a silicon transistor,” Science Advances, vol. 4, no. 12, p. eaat9199, 2018.
  * [27] D. V. Bulaev and D. Loss, “Electric dipole spin resonance for heavy holes in quantum dots,” Physical Review Letters, vol. 98, no. 9, p. 097202, 2007\.
  * [28] J. Salfi, M. Tong, S. Rogge, and D. Culcer, “Quantum computing with acceptor spins in silicon,” Nanotechnology, vol. 27, no. 24, p. 244001, 2016.
  * [29] D. V. Bulaev and D. Loss, “Spin relaxation and decoherence of holes in quantum dots,” Physical Review Letters, vol. 95, no. 7, p. 076805, 2005.
  * [30] C. Testelin, F. Bernardot, B. Eble, and M. Chamarro, “Hole–spin dephasing time associated with hyperfine interaction in quantum dots,” Physical Review B, vol. 79, no. 19, p. 195440, 2009.
  * [31] P. Harvey-Collard, N. T. Jacobson, C. Bureau-Oxton, R. M. Jock, V. Srinivasa, A. M. Mounce, D. R. Ward, J. M. Anderson, R. P. Manginell, J. R. Wendt, et al., “Spin-orbit interactions for singlet-triplet qubits in silicon,” Physical Review Letters, vol. 122, no. 21, p. 217702, 2019.
  * [32] T. Tanttu, B. Hensen, K. W. Chan, C. H. Yang, W. W. Huang, M. Fogarty, F. Hudson, K. Itoh, D. Culcer, A. Laucht, et al., “Controlling spin-orbit interactions in silicon quantum dots using magnetic field direction,” Physical Review X, vol. 9, no. 2, p. 021028, 2019.
  * [33] Y.-L. Hsueh, D. Keith, Y. Chung, S. K. Gorman, L. Kranz, S. Monir, Z. Kembrey, J. G. Keizer, R. Rahman, and M. Y. Simmons, “Engineering spin-orbit interactions in silicon qubits at the atomic-scale,” Advanced Materials, p. 2312736, 2024.
  * [34] M. Nestoklon, E. Ivchenko, J.-M. Jancu, and P. Voisin, “Electric field effect on electron spin splitting in si ge/ si quantum wells,” Physical Review B, vol. 77, no. 15, p. 155328, 2008.
  * [35] R. M. Jock, N. T. Jacobson, P. Harvey-Collard, A. M. Mounce, V. Srinivasa, D. R. Ward, J. Anderson, R. Manginell, J. R. Wendt, M. Rudolph, et al., “A silicon metal-oxide-semiconductor electron spin-orbit qubit,” Nature Communications, vol. 9, no. 1, pp. 1–8, 2018.
  * [36] S. Lee, H. Koike, M. Goto, S. Miwa, Y. Suzuki, N. Yamashita, R. Ohshima, E. Shigematsu, Y. Ando, and M. Shiraishi, “Synthetic rashba spin–orbit system using a silicon metal-oxide semiconductor,” Nature Materials, vol. 20, no. 9, pp. 1228–1232, 2021.
  * [37] S. Liles, R. Li, C. Yang, F. Hudson, M. Veldhorst, A. S. Dzurak, and A. Hamilton, “Spin and orbital structure of the first six holes in a silicon metal-oxide-semiconductor quantum dot,” Nature Communications, vol. 9, no. 1, pp. 1–7, 2018.
  * [38] G. Klimeck, S. S. Ahmed, H. Bae, N. Kharche, S. Clark, B. Haley, S. Lee, M. Naumov, H. Ryu, F. Saied, et al., “Atomistic simulation of realistically sized nanodevices using nemo 3-d—part i: Models and benchmarks,” IEEE Transactions on Electron Devices, vol. 54, no. 9, pp. 2079–2089, 2007.
  * [39] B. Weber, S. Mahapatra, T. F. Watson, and M. Y. Simmons, “Engineering independent electrostatic control of atomic-scale ( 4 nm) silicon double quantum dots,” Nano Letters, vol. 12, no. 8, pp. 4001–4006, 2012.
  * [40] M. Fuechsle, J. A. Miwa, S. Mahapatra, H. Ryu, S. Lee, O. Warschkow, L. C. Hollenberg, G. Klimeck, and M. Y. Simmons, “A single-atom transistor,” Nature Nanotechnology, vol. 7, no. 4, p. 242, 2012.
  * [41] B. Weber, Y. M. Tan, S. Mahapatra, T. F. Watson, H. Ryu, R. Rahman, L. C. Hollenberg, G. Klimeck, and M. Y. Simmons, “Spin blockade and exchange in coulomb-confined silicon double quantum dots,” Nature Nanotechnology, vol. 9, no. 6, p. 430, 2014.
  * [42] S. Tarucha, D. Austing, T. Honda, R. Van der Hage, and L. P. Kouwenhoven, “Shell filling and spin effects in a few electron quantum dot,” Physical Review Letters, vol. 77, no. 17, p. 3613, 1996.
  * [43] H. Liu, T. Fujisawa, T. Hayashi, and Y. Hirayama, “Pauli spin blockade in cotunneling transport through a double quantum dot,” Physical Review B, vol. 72, no. 16, p. 161305, 2005.
  * [44] W. Coish and F. Qassemi, “Leakage-current line shapes from inelastic cotunneling in the pauli spin blockade regime,” Physical Review B, vol. 84, no. 24, p. 245407, 2011.
  * [45] N. Lai, W. Lim, C. Yang, F. Zwanenburg, W. Coish, F. Qassemi, A. Morello, and A. Dzurak, “Pauli spin blockade in a highly tunable silicon double quantum dot,” Scientific Reports, vol. 1, p. 110, 2011.
  * [46] L. Kranz, S. K. Gorman, B. Thorgrimsson, Y. He, D. Keith, J. G. Keizer, and M. Y. Simmons, “Exploiting a single-crystal environment to minimize the charge noise on qubits in silicon,” Advanced Materials, vol. 32, no. 40, p. 2003361, 2020.
  * [47] S. Moriyama, J. Martinek, G. Ilnicki, T. Fuse, and K. Ishibashi, “Inelastic cotunneling mediated singlet-triplet transition in carbon nanotubes,” Physical Review B, vol. 80, no. 3, p. 033408, 2009.
  * [48] M. Borhani and X. Hu, “Two-spin relaxation of p dimers in silicon,” Physical Review B, vol. 82, no. 24, p. 241302, 2010.
  * [49] A. Johnson, J. R. Petta, J. Taylor, A. Yacoby, M. Lukin, C. Marcus, M. Hanson, and A. Gossard, “Triplet–singlet spin relaxation via nuclei in a double quantum dot,” Nature, vol. 435, no. 7044, pp. 925–928, 2005.
  * [50] J. R. Petta, A. C. Johnson, J. M. Taylor, E. A. Laird, A. Yacoby, M. D. Lukin, C. M. Marcus, M. P. Hanson, and A. C. Gossard, “Coherent manipulation of coupled electron spins in semiconductor quantum dots,” Science, vol. 309, no. 5744, pp. 2180–2184, 2005.
  * [51] E. Laird, J. R. Petta, A. C. Johnson, C. Marcus, A. Yacoby, M. Hanson, and A. Gossard, “Effect of exchange interaction on spin dephasing in a double quantum dot,” Physical Review Letters, vol. 97, no. 5, p. 056801, 2006\.
  * [52] J. R. Petta, J. Taylor, A. Johnson, A. Yacoby, M. Lukin, C. Marcus, M. Hanson, and A. Gossard, “Dynamic nuclear polarization with single electron spins,” Physical Review Letters, vol. 100, no. 6, p. 067601, 2008.
  * [53] T. Meunier, I. Vink, L. W. van Beveren, K. Tielrooij, R. Hanson, F. Koppens, H. Tranitz, W. Wegscheider, L. Kouwenhoven, and L. Vandersypen, “Experimental signature of phonon-mediated spin relaxation in a two-electron quantum dot,” Physical Review Letters, vol. 98, no. 12, p. 126601, 2007\.
  * [54] S. Amasha, K. MacLean, I. Radu, D. Zumbuhl, M. Kastner, M. Hanson, and A. Gossard, “Measurements of the spin relaxation rate at low magnetic fields in a quantum dot,” arXiv preprint cond-mat/0607110, 2006.
  * [55] C. Yang, A. Rossi, R. Ruskov, N. Lai, F. Mohiyaddin, S. Lee, C. Tahan, G. Klimeck, A. Morello, and A. Dzurak, “Spin-valley lifetimes in a silicon quantum dot with tunable valley splitting,” Nature Communications, vol. 4, no. 1, pp. 1–8, 2013.
  * [56] M. Fuechsle, S. Mahapatra, F. Zwanenburg, M. Friesen, M. Eriksson, and M. Y. Simmons, “Spectroscopy of few-electron single-crystal silicon quantum dots,” Nature Nanotechnology, vol. 5, no. 7, p. 502, 2010.
  * [57] C. Tahan and R. Joynt, “Relaxation of excited spin, orbital, and valley qubit states in ideal silicon quantum dots,” Physical Review B, vol. 89, no. 7, p. 075302, 2014.
  * [58] G. Klimeck, S. Ahmed, H. Bae, N. Kharche, S. Clark, B. Haley, S. Lee, M. Naumov, H. Ryu, F. Saied, M. Prada, M. Korkusinski, T. Boykin, and R. Rahman, “Atomistic simulation of realistically sized nanodevices using NEMO 3-d—part i: Models and benchmarks,” IEEE Transactions on Electron Devices, vol. 54, pp. 2079–2089, Sept. 2007.
  * [59] R. Rahman, G. Lansbergen, J. Verduijn, G. Tettamanzi, S. Park, N. Collaert, S. Biesemans, G. Klimeck, L. Hollenberg, and S. Rogge, “Electric field reduced charging energies and two-electron bound excited states of single donors in silicon,” Physical Review B, vol. 84, no. 11, p. 115428, 2011\.
  * [60] D. Chadi, “Spin-orbit splitting in crystalline and compositionally disordered semiconductors,” Physical Review B, vol. 16, no. 2, p. 790, 1977.
  * [61] T. B. Boykin, G. Klimeck, and F. Oyafuso, “Valence band effective-mass expressions in the sp 3 d 5 s* empirical tight-binding model applied to a si and ge parametrization,” Physical Review B, vol. 69, no. 11, p. 115201, 2004.


*arg max *arg min Department of Computer Science, Stanford University, USA and https://cs.stanford.edu/people/paulliu/paul.liu@stanford.eduhttps://orcid.org/0000-0002-9386-6609 Department of Computer Science, University of Waterloo<EMAIL_ADDRESS>No<EMAIL_ADDRESS>Department of Computer Science, University of Illinois at Urbana-Champaign, USA and https://davidzheng.web.illinois.edu/dwzheng2@illinois.eduhttps://orcid.org/0000-0002-0844-9457 Paul Liu, Jack Spalding-Jamieson, Brandon Zhang, Da Wei Zheng ¡ccs2012¿ ¡concept¿ ¡concept_id¿10010147.10010178.10010213.10010215¡/concept_id¿ ¡concept_desc¿Computing methodologies Motion path planning¡/concept_desc¿ ¡concept_significance¿500¡/concept_significance¿ ¡/concept¿ ¡concept¿ ¡concept_id¿10003752.10010061.10010063¡/concept_id¿ ¡concept_desc¿Theory of computation Computational geometry¡/concept_desc¿ ¡concept_significance¿500¡/concept_significance¿ ¡/concept¿ ¡/ccs2012¿ [500]Computing methodologies Motion path planning [500]Theory of computation Computational geometry github.com/jacketsj/cgshop2021-gitastrophe Kevin Buchin and Éric Colin de Verdière 2 37th International Symposium on Computational Geometry (SoCG 2021) SoCG 2021 SoCG 2021 June 7–11, 2021 Buffalo, NY, USA socg-logo.pdf 189 64 83

# Coordinated Motion Planning Through Randomized $k$-Opt

Paul Liu  Jack Spalding-Jamieson  Brandon Zhang  Da Wei Zheng

###### Abstract

This paper examines the approach taken by team gitastrophe in the CG:SHOP 2021
challenge. The challenge was to find a sequence of simultaneous moves of
square robots between two given configurations that minimized either total
distance travelled or makespan (total time). Our winning approach has two main
components: an initialization phase that finds a good initial solution, and a
$k$-opt local search phase which optimizes this solution. This led to a first
place finish in the distance category and a third place finish in the makespan
category.

###### keywords:

motion planning, randomized local search, path finding

###### category:

CG challenge

##  1 Introduction

For a set of unit square robots ${\cal R}$ each with start and target
locations $\\{s_{r}\\}_{r\in{\cal R}},\\{d_{r}\\}_{r\in{\cal
R}}\subset{\mathbb{Z}}^{2}$, and a set of unit square obstacles
$O\subset{\mathbb{Z}}^{2}$, the _coordinated motion planning_ problem asks for
an optimal sequence of simultaneous “moves” for each robot that brings the
robots from their start location to their target location. At each timestep,
robots can move to an adjacent grid cell or stand still (a no-op). These moves
are subject to the following constraints at all timesteps:

  * •

No robot shares a location with another robot or an obstacle.

  * •

No robot moves into a location previously occupied by another robot, unless
the other robot is moving in the same direction at that timestep.

An _optimal_ sequence of moves can refer to two different criteria:

  * •

MAX: The minimum makespan, i.e. the total number of timesteps needed.

  * •

SUM: The minimum total distance, i.e. the total number of position-changing
moves each robot takes.

We, team gitastrophe, explored this problem as a part of the 2021
Computational Geometry Challenge (CG:SHOP 2021). The challenge ranked teams
according to each of the two different optimality criteria. Our team ranked
first among all junior teams, first according to the total distance criterion,
and third according to the makespan criterion. Team Shadoks [1] ranked first
according to makespan and third according to total distance, while Team UNIST
[3] ranked second in both categories.

One important property of the test instances of the challenge is that
feasibility is always guaranteed — no obstacles enclose the starting or target
positions of any robots. In the following sections, we assume this property.
We refer the reader to the survey of the challenge [2] for more details on the
instances.

##  2 Methods

As with other teams who took part in the challenge, our strategy consisted of
two phases: initialization (sec:init) and solution optimization (sec:k-opt).

The initialization phase consisted of moving the robots out of the initial
problem grid, into an intermediate state that can easily be routed to their
end goals. This strategy was remarkably similar to the feasibility approaches
used by both Team Shadoks [1] and Team UNIST [3]. The solution optimization
phase used a novel local search strategy, which consisted of locally
reordering the movements of $k$ sampled robots while keeping all others fixed.

###  2.1 Initialization strategy

Our initialization strategy is simple and relies on the notion of _depth_
values (see def:depth). The depth values partition the robots into subsets
which can be routed in parallel. Intuitively, once a robot is routed, it will
not interfere with robots of smaller depth value.

((a)) small_005 start and intermediate positions

((b)) small_005 target and intermediate positions

((c)) medium_005 start and intermediate positions

((d)) medium_005 target and intermediate positions

Figure 1: A visualization of our initializations for the instances small_005
and medium_005. The grey squares are obstacles. On the left (Figures 1 and 1),
the colours of the filled boxes match the instances start positions to the
intermediate positions. On the right, (Figures 1 and 1), the colours match the
target positions and intermediate positions.

((a)) Octagon filler

((b)) Diamond filler

((c)) Quad rectangle filler

((d)) Rectangle filler

Figure 2: Visualizations of the different methods we used to generate
intermediate positions for the instance small_012. This visualization only
shows the matched potential intermediate positions, which is why the displayed
intermediate positions in Figure 2 do not resemble a rectangle.

######  Definition 2.1.

Given a set of locations $H\subset{\mathbb{Z}}^{2}$, the depth value $D(v)$
for $v\in{\mathbb{Z}}^{2}$ is computed recursively:

  * •

If $v\in H$, then $D(v)=0$.

  * •

Otherwise, $D(v)=1+\min_{u}D(u)$ where the $\min$ is over all non-obstacle
squares adjacent to $v$.

Such values are uniquely defined and can easily be found via breadth-first
search.

The initialization operates in a few phases:

  1. 1.

Compute a set of intermediate positions $H$, derived from a superset of
candidate positions which we call _filler_ shapes (see Figure 2). These filler
shapes are located outside of the bounding box of the robots’ start positions,
and arranged in a way such that no intermediate positions are adjacent.

  2. 2.

Compute a min-cost matching of the robots to their intermediate positions. The
cost of a matching is the sum of distances from the robots’ start and target
positions to their matched intermediate position. To improve the
initialization, we generate many more intermediate positions than robots.
Figure 1 shows an example of the robots’ start and target positions and the
corresponding matched intermediate positions.

  3. 3.

For each coordinate, compute the depth of that coordinate using $H$ (see
def:depth). Route the robots from their start position to their corresponding
intermediate positions one at a time, in order according to decreasing depth
values of their start positions. The routes here are constructed to avoid the
paths of the robots routed before it. In this phase we route robots with
starting locations of larger depth before smaller ones.

  4. 4.

Re-route the robots _from their start positions_ to their target positions’ in
decreasing order of their depth values. We start with the paths computed in
the previous phase, and reroute each robot individually. We do this by
removing their path and inserting a new path respecting the current paths of
all robots.

Although a natural alternative to step 4 is to simply route the robots from
the intermediate locations to their target locations, it is instead much more
efficient to re-route from the robots’ starting locations. Given the paths
produced in the third step, a feasible set of paths is guaranteed to exist in
the fourth step, since it’s clear that there exists a set of paths from each
robots’ intermediate position to their target positions: All that is needed is
a path from the intermediate position to the target position that works after
all other robots have finished moving. Such a path exists since robots
residing at higher depth values do not block the paths of robots at lower
depth values.

#### Initialization variations

Local search strategies are typically very sensitive to the initialization. We
used a number of variations to generate a family of different initializations:

  * •

Using different filler shapes in step 1, or different costs for matching in
step 3.

  * •

Finding random shortest paths, or approximate shortest paths for steps 3 and
4.

  * •

Swapping the robots’ start and target locations.

The optimization procedure was run on the initializations that scored the
highest.

###  2.2 Solution optimization

As previously described, the basic paradigm of our local optimization is to
choose $k$ robots, and improve their paths while keeping the paths of all
others fixed. This type of approach is often referred to as $k$-opt. This
approach has two distinct components: (1) the choice of the $k$ robots and (2)
the improvements of their paths.

#### Choosing the $k$ robots

We experimented with a few types of weighted sampling to choose the $k$
robots:

  * •

_Sampling by completion time._ Sample $k$ robots without replacement, where
the sampling probability for each robot is proportional to their completion
time (the time at which they make their last move).

  * •

_Sampling by closeness._ Sample the initial robot proportional to completion
time, then sample $k-1$ other robots based on proximity to the initial robot’s
path. The proximity of two paths is the number of time steps for which the two
paths are adjacent to each other.

  * •

_Sampling by constraints._ Sample the initial robot proportional to completion
time, then find a minimum completion time path from the start to the target
for the sampled robot with the relaxation that this robot is allowed to move
through $k-1$ other robots. The robots that the path moves through forms the
$k-1$ other robots in the sample.

In our experiments, sampling by completion time and by closeness seemed to
produce the best results. However, we note that team Shadoks was able to
exploit a variant of sampling by constraints to win the MAX category.

#### Path optimization

Given a sample of $k$ robots, one would like to jointly optimize the $k$
robots simultaneously. However, this approach causes the state space to grow
exponentially with the number of robots. Furthermore, due to the size of the
grid, the state space for path finding is already quite large to begin with.
Like other teams, our path finding was done in the grid-time graph where
states are characterized by the positions of the robots at each time $t$. For
the largest instance in the data set, the size of this state space is already
on the order of $10^{5}$. For these reasons, we were only able to scale to
$k=2$ in our code with the approach of joint optimization.

Instead, our main insight was to _approximate_ the joint optimization by $k$
individual single path optimizations. Our inspiration was the analogy of
sorting: given the $k$ robots, we route the robots one-by-one to their
targets, where the $j$-th robot routed ($j\leq k$) respects the path of robots
$1\ldots j-1$. Crucially, this avoids the exponential blow-up of the state
space as we’re routing each robot from start to target in succession. The
problem then reduces to finding a good ordering of the $k$ robots to reroute.
We had the most success by simply ordering the robots by decreasing completion
time. If an ordering was infeasible (e.g. due to some subset of the $k$ robots
completely blocking off the path from the remaining robots) or not an
improvement, it was discarded.

During the competition we found that there were advantages to relatively
smaller values of $k$ for faster computation. For our approximate joint
optimization, larger $k$ helped to get out of bad local optima. We iteratively
used all values of $k$ between $1$ and $7$ during the challenge.

Since the main component of our algorithm is path-finding, we used a number of
practical techniques to reduce its cost. We use A* with a bucket priority
queue, where the heuristic function was taken to be either the Manhattan
distance or the shortest path distance in the graph with only obstacles
without robots. To further speed up the search when we have an initial
feasible solution, we limit the path-finding algorithm to search locally
around the robot’s original path, by enforcing that no robot may deviate more
than $R$ steps away from the set of positions forming their initial path for a
fixed parameter $R$. Since the search is centred around the initial solution,
the solution space is guaranteed to have a feasible solution for any value of
$R$.

The optimizations above are agnostic to the optimization objective and allowed
for large gains following the initialization phase for both SUM and MAX.

##  3 Results

###  3.1 Computational environment

For the most part, our experiments were done on a desktop with a Core i7-2600.
In the last three weeks of the competition we used Stanford’s Sherlock High-
Performance Computing Cluster for SUM optimization. The University of
Waterloo’s Multicore lab also generously donated some night-time compute
during the last week of the competition in the form of two EPYC 7662s.

###  3.2 Experimental results

((a)) SUM optimization varying $k$ where $R=20$.

((b)) SUM optimization varying $R$ where $k=7$.

Figure 3: Varying $k$ and $R$ for SUM optimization on microbes_004. Final SUM:
43437.

((a)) MAX optimization varying $k$ where $R=20$.

((b)) MAX optimization varying $R$ where $k=7$.

Figure 4: Varying $k$ and $R$ for MAX optimization on microbes_004. Final MAX:
126.

To empirically determine the effect of $k$ and $R$ had on the optimization, we
ran our algorithm on the instance microbes_00004. As seen in Figure 3,
different values of $k$ had an effect on the rate of convergence and the value
of the local minimum found by our strategy. Experimentally, values of $k$
greater than $10$ did not improve our scores and also ran slower. For MAX (in
Figure 4), improvements were much more discrete.

As we varied $R$, we had a similar trade-off between performance and runtime
for each optimization step of our algorithm (Figures 3 and 4). Empirically,
choosing $R$ to be around $20$ achieved the best balance for most of the
instances.

Figure 5 compares the solutions we computed during the competition to the
trivial lower bound, which is given by the maximum shortest-path distance for
SUM and the sum of all shortest-path distances for MAX. There is a clear trend
that score decreases as density increases (and as $n$ increases, to a lesser
extent). This could be due to two factors: first, that the trivial lower bound
becomes a worse approximation of the optimal score as these values increase;
second, that our algorithm performs poorer on these instances, because more
computational time is required and local changes to $k$ robots at a time are
insufficient to traverse the solution space.

Figure 5: Plots of $\frac{\text{score}}{\text{lowerbound}}$ versus $n$ for
both SUM and MAX. Each point corresponds to one instance. The density of the
instance, as defined by the contest organizers, is indicated by colour.

## References

  * [1] Loïc Crombez, Guilherme D. da Fonseca, Yan Gerard, Aldo Gonzalez-Lorenzo, Pascal Lafourcade, and Luc Libralesso.  Shadoks coordinated motion planning.  2021\.
  * [2] Sándor Fekete, Phillip Keldenich, Dominik Krupke, and Joseph S. B. Mitchell.  Computing coordinated motion plans for robot swarms: The cg:shop challenge 2021.  37th International Symposium on Computational Geometry, SoCG 2021, 2021.
  * [3] Hyeyun Yang and Antoine Vigneron.  A simulated annealing approach to coordinated motion planning.  2021\.