# Spin hierarchy in van der Waals molecule formation via ultracold three-body
recombination

Jing-Lun Li  Institut für Quantenmaterie and Center for Integrated Quantum
Science and Technology IQST, Universität Ulm, 89069 Ulm, Germany  Paul S.
Julienne  Joint Quantum Institute, University of Maryland, and the National
Institute of Standards and Technology (NIST), College Park, MD 20742, USA
Johannes Hecker Denschlag  Institut für Quantenmaterie and Center for
Integrated Quantum Science and Technology IQST, Universität Ulm, 89069 Ulm,
Germany  José P. D’Incao  JILA, NIST, and the Department of Physics,
University of Colorado, Boulder, CO 80309, USA

###### Abstract

We theoretically investigate the product-state distribution of weakly bound
diatomic van der Waals molecules via ultracold three-body recombination of
bosonic alkali atoms. We find a two-level hierarchy of spin propensity rules
at zero magnetic field. The primary propensity rule states that nearly all
molecular products conserve the total hyperfine spin of reactant atomic pairs,
while molecular products not conserving the total spin are highly suppressed.
For the dominant molecular products, there is a secondary propensity to
conserve certain spin components of the reactant pair such as the atomic
hyperfine spins, or the total electronic or nuclear spins. The second
propensity varies across species and depends fundamentally on the interplay
between effective electronic exchange and hyperfine interactions. The spin
sensitivity of product-state distribution can potentially open up new avenues
for controlling state-to-state reaction rates in ultracold three-body
recombination.

Despite the intrinsically complex microscopic properties of the interatomic
interactions and symmetries [1, 2, 3], chemical reactions are often governed
by unexpectedly general and simple fundamental principles, conservation laws
and propensity rules [4, 5, 6, 7, 8]. The ability to study in detail the
likelihood of obtaining a particular product state from well-defined chemical
reactants is crucial for advancing our knowledge and uncovering novel
principles and mechanisms controlling chemical reactions. Due to recent
advancement in ultracold atomic and molecular gases, experiments can now
prepare reactants in a well-defined quantum state and detect products
resolving all quantum degrees of freedom including vibration, rotation,
electronic spin, and nuclear spin. This state-to-state resolution of chemical
reactions allows for in-depth investigations of novel propensity rules in an
unprecedented level of detail [9, 10, 11, 12, 13, 14, 15, 16]. Understanding
such fundamental reaction principles is crucial for developing control over
chemical reactions and their products [17, 18, 19, 20, 21, 22].

In ultracold atomic gases a key exothermic chemical reaction is three-body
recombination, where three atoms collide to form a diatomic molecule and a
free atom. This chemical reaction, with rate constant $L_{3}$, drastically
limits the lifetime and stability of Bose-Einstein condensates [23, 24, 25,
26, 27, 28, 29, 30, 31, 32, 33] and has been extensively used as a probe to
explore fundamental few-body phenomena, such as the Efimov effect [34, 35, 36,
37]. In recent years, detailed insights on recombination have been gained by
examining its molecular product-state distribution in ultracold Rb gases [10,
11, 12, 13]. Such studies have found that the most populated molecular
products are weakly bound molecules, conventionally referred to as van der
Waals (vdW) molecules [38, 39, 40, 41]. Over decades, the formation and
reactions involving vdW molecules have been extensively studied across a wide
range of physical chemistry processes [42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
52]. Recent studies [10, 11, 12, 13] have revealed that, at zero magnetic
field, the Rb2 vdW molecules produced via recombination possess the same
hyperfine spins as the initial atoms and that the corresponding state-to-state
reaction rates are roughly proportional to the inverse of their binding
energy, $1/E_{b}$. Nevertheless, a deeper understanding of such propensity
rules is necessary to assess their validity among other atomic species with
vastly different physical properties. In particular, the role of the intricate
relationship between the long- and short-range physics in the spin dynamics of
three-body recombination remains unclear. At large distances, hyperfine and
vdW interactions govern the spin dynamics while at short distances the
dynamics is dominated by the electronic exchange interactions. The interplay
between hyperfine and electronic exchange interactions, for instance,
determines the spin mixing in vdW molecules, potentially influencing how such
molecular states are formed via recombination.

In this Letter we present a theoretical study on the product-state
distribution of vdW molecules formed by ultracold three-body recombination to
characterize the role of the molecular spin in the reaction dynamics. We
consider the bosonic alkali atoms (7Li, 23Na, 39K, 41K, 85Rb and 87Rb and
133Cs) at zero magnetic field. We find a two-level hierarchy of spin
propensity rules. The first rule states that nearly all molecular products
conserve the total hyperfine spin of the reactant atomic pair. The rates for
such reactions are generally consistent with the $1/E_{b}$ propensity rule
[13]. In contrast, the formation of molecules not conserving the total
hyperfine spin is highly suppressed and violates the $1/E_{b}$ propensity
rule. As we will show, this points to a different formation process than that
of the dominant molecular states. The second propensity rule states express
the likelihood to conserve certain spin components of the reactant pair such
as the atomic hyperfine spins, or the total electronic or nuclear spins. We
introduce a dimensionless parameter $\xi_{\rm ex}$ that indicates both the
propensity strength and the kinds of spin components to which the propensity
applies. Here, $\xi_{\mathrm{ex}}$ is the ratio of the effective exchange and
hyperfine interactions of the atomic pair and is an improved generalization of
the parameter $\rho$ previously introduced in [12]. It determines the spin
structure of near-threshold molecular states and varies strongly for different
atomic species. For details of how $\xi_{\rm ex}$ is defined see the
Supplementary Material [53] and Ref. [54].

Our three-body studies are performed in the adiabatic hyperspherical
representation [55, 56], where the hyperradius $R$ gives the overall size of
the system while the set of hyperangles $\Omega$ describe its internal motion.
The adiabatic separation between hyperradial and hyperangular motions leads to
the three-body hyperradial Schrödinger equation

$\displaystyle\left[-\frac{\hbar^{2}}{2\mu}\frac{d^{2}}{dR^{2}}+U_{\nu}(R)\right]F_{\nu}(R)$
$\displaystyle~{}~{}~{}~{}~{}~{}+\sum_{\nu^{\prime}}W_{\nu\nu^{\prime}}(R)F_{\nu^{\prime}}(R)=EF_{\nu}(R),$
(1)

governing the hyperradial motion via the three-body potentials $U_{\nu}$ with
inelastic transitions driven by the nonadiabatic couplings
$W_{\nu\nu^{\prime}}$. Here, $E$ is the total energy, $\mu=m/\sqrt{3}$ the
three-body reduced mass, $m$ the atomic mass, and $\nu$ the set of quantum
numbers necessary to characterize each channel. In the present study we
include the atomic hyperfine structure and interatomic interactions given by
the electronic singlet and triplet Born-Oppenheimer potentials [57, 58, 59,
60, 61]. Such potentials are modified to restrict their number of molecular
states while still preserving the proper singlet and triplet scattering
lengths, $a_{s}$ and $a_{t}$, respectively [62, 63, 54]. In Fig. 1 we display
the three-body potentials for 87Rb atoms and describe the corresponding
process leading to three-body recombination. We note that the vdW interaction
between atoms defines the typical scale for length $r_{\rm
vdW}=\frac{1}{2}(mC_{6}\hbar^{2})^{1/4}$ and energy $E_{\rm
vdW}=\hbar^{2}/mr_{\rm vdW}^{2}$ of the system, where $C_{6}$ is the vdW
dispersion coefficient. Compared to this scale, our considered collisional
energy of free atoms $E$ is extremely small $E\ll E_{\rm vdW}$, while the vdW
molecular binding energies can be up to thousands of $E_{\rm vdW}$. In what
follows we will use vdW units to allow a direct comparison between different
atomic species.

Figure 1: Three-body adiabatic potentials, $U(R)$, for 87Rb atoms in van der
Waals units. In three-body recombination, three free atoms, characterized by
the spin product state $|f_{*}$=1
$m_{f}^{*}$=-1$\rangle$$|f_{*}m^{*}_{f}\rangle|f_{*}m^{*}_{f}\rangle$,
approach each other via the three-body continuum channels [i.e., channels
where $U(R)$$>$$0$ for $R$$\gg$$r_{\rm vdW}$]. As the atoms reach the
“reaction region” ($R\lesssim r_{\rm vdW}$ shaded area) the complex structure
of avoid-crossings drive inelastic transitions to the final atom-molecule
channels [$U(R)$$<$0 for $R$$\gg$$r_{\rm vdW}$] characterized by the
$|\Psi_{FM_{F}}(v,L_{R})\rangle|fm_{f}\rangle$ state. The different colors of
outcoming arrows indicate the variation of the spin state of molecular
products.

In the hyperspherical representation the reaction rates are obtained from the
solutions of Eq. (1), from which we determine the scattering $S$-matrix [56].
The total recombination rate, $L_{3}$, and the corresponding state-to-state
rates, $L_{3\beta}$, are given by [37]

$\displaystyle L_{3}=\frac{1}{2}\sum_{\alpha\beta}\frac{192\pi^{2}\hbar}{\mu
k^{4}}|S_{\beta\alpha}|^{2}=\sum_{\beta}L_{3\beta},$ (2)

where $k^{2}=2\mu E/\hbar^{2}$ and with $\alpha$ and $\beta$ running over
initial (three-body continuum) and final (atom-molecule) channels,
respectively. In this study, we analyze the state-to-state recombination rates
$L_{3\beta}$ according to the binding energy and spins of the corresponding
molecular product states. For all atomic species we choose the initial atomic
state to be the ”spin-stretched” state of the lowest hyperfine manifold,
$|fm_{f}\rangle$$\equiv$$|f_{*}m^{*}_{f}$=-$f_{*}\rangle$, where $f$ is the
atomic hyperfine spin and $m_{f}$ its azimuthal projection.Thus, we have
$f_{*}=3$ for 133Cs, $f_{*}=2$ for 85Rb and $f_{*}=1$ for all others (7Li,
23Na, 39K, 41K, and 87Rb). With this choice fast two-body losses are
prevented, providing a convenient experimental condition for the study of
product-state distribution. We note that our results are also valid for the
state $|f_{*}m^{*}_{f}$=$+f^{*}\rangle$.

Figure 2:  Recombination fraction $L_{3\beta}/L_{3}$ characterizing the
product-state distribution in terms of the molecular binding energies,
$E_{\beta}$, and spins at vanishing magnetic field for different atomic
species. Here, we indicate the initial atomic spin state for each species,
$|f_{*}m_{f}^{*}\rangle$, while the final molecular states (see legends) are
classified according to their molecular spin $|FM_{F}(f_{a}f_{b})\rangle$ or,
in some cases, $|FM_{F}(\chi_{1})\rangle$ or $|FM_{F}(\chi_{2})\rangle$ for a
mixed molecular state (see [64] for a detail description of our assignment).
The lower panels share the same legend. The straight solid lines denote the
$1/E_{\beta}$ scaling of $L_{3\beta}$, while the dotted lines indicate a
constant rate. The shaded area indicates the region where the partial rates
typically scatter around the $1/E_{\beta}$ scaling. For each atomic species,
we indicate the corresponding values of $L_{3}$ and $\xi_{\rm ex}$ (see
discussion in the text).

While the initial spin state for recombination is always a product state
$|f_{*}m^{*}_{f}\rangle|f_{*}m^{*}_{f}\rangle|f_{*}m^{*}_{f}\rangle$, the
final atom-molecule state is characterized by
$|\Psi_{FM_{F}}(v,L_{R})\rangle|fm_{f}\rangle$ where $\Psi_{FM_{F}}(v,L_{R})$
is the total (spatial and spin) molecular wavefunction with ($v,L_{R}$)
representing the vibrational and rotational molecular quantum numbers [65].
Our studies are performed at zero $B$-field such that the two-atom total
hyperfine spin, $F$ ($|f_{a}-f_{b}|\leq F\leq f_{a}+f_{b}$), and its
projection, $M_{F}=m_{f_{a}}+m_{f_{b}}$, are good quantum numbers for
characterizing the molecular states. However, since in a three-body system the
two-body total angular momentum quantum numbers $F$ and $M_{F}$ are not in
general conserved quantities, the question is what rules control the product
distribution in a dynamical process like three-body recombination. Similar to
Rb [12, 13], our findings across different atomic species provide clear
evidence of a spin propensity rule favoring recombination into molecular
states, $|\Psi_{FM_{F}}(v,L_{R})\rangle$, whose values for $F$ and $M_{F}$ are
the same as that of the atomic pairs forming the initial state, thus implying
the quasi-conservation of $(F,M_{F})$. This is the first, and strongest, rule
of our propensity hierarchy. In our case, since the initial atomic state is
$|f_{*}m_{f}^{*}\rangle$, any pair of atoms will have $F_{*}=2f_{*}$ and
$M_{F}^{*}=-2f_{*}$ and, as a result, the favorable molecular states are those
conserving $(F_{,}M_{F})$, i.e., with $(F,M_{F})=(F_{*},M_{F}^{*})$.

Our numerical results are shown in Fig. 2 for the recombination fraction,
$L_{3\beta}/L_{3}$, for each molecular state according to their corresponding
binding energies, $E_{\beta}$. We classify the states according to their
molecular spin $|FM_{F}(f_{a}f_{b})\rangle$ or, in some cases,
$|FM_{F}(\chi)\rangle$ for mixed molecular states, i.e., molecular states
whose spins are in a superposition of $|FM_{F}(f_{a}f_{b})\rangle$ states of
same $F$ and $M_{F}$ values [64]. The molecular states included in our model
[54] have $(v,L_{R})$ values (not shown in Fig. 2) ranging from $v=-1$ to -7
(counting down from the most weakly bound state) and $L_{R}=0$ up to 26, in a
total of about $100$ molecular states. We note that the data in Fig. 2
corresponding to our most deeply bound molecular states are likely to be
model-dependent since they don’t fully reproduce the actual molecular states
within that range of energies [54]. Figure 2 shows that for molecular states
with $(F_{*},M_{F}^{*})$ the rates are generally consistent with the
$L_{3\beta}\propto 1/E_{\beta}$ energy scaling propensity rule of Ref. [13]
(see $1/E_{\beta}$ solid lines in Fig. 2). However, our data set encounters
insufficient binding energy range and large fluctuations, making it
challenging to achieve a reliable fitting. (To more accurately test the
$1/E_{\beta}$ propensity rule it would be necessary to include significantly
more deeply bound molecular states, as done in the spinless model of Ref.
[13]. Nevertheless, for our spin-dependent model this would lead to a
considerable increase in the numerical demand.) For the molecular states whose
($F,M_{F}$) are not conserved, however, $L_{3\beta}$ is generally much smaller
and depends only weakly on $E_{\beta}$ (see horizontal dashed lines in Fig.
2). We also have found that, similarly to our previous studies in Refs. [10,
12, 13], the partial rates have no apparent dependence on $L_{R}$. It should
also be noted that the results in Fig. 2 were carried out by fixing the spin
of the third atom at $|f_{*}m^{*}_{f}\rangle$. However, we have carefully
examined the contribution from other spin states of the third atom and found
it is typically below one percent. This points to a largely mechanical role of
the third atom in recombination (see Supplemental Material [53]). Our result
is in contrast to the strong three-body spin-exchange recombination of 7Li at
large magnetic fields found in Refs. [66, 67]. We attribute this difference to
the fact that at finite magnetic fields $F$ and even $f$ are not good quantum
numbers anymore. A systematic study on the effects of finite magnetic fields
in the product state distribution is beyond the scope of the present study.

We now analyse the propensity to conserve spin components of the reacting
atomic pair. This leads to the second rule in the propensity hierarchy,
involving the $\xi_{\rm ex}$ parameter. In fact, for very small and very large
$\xi_{\rm ex}$ parameters the production rate depends on the spin overlap of
the initial atom pair and the molecular product.

Figure 2 shows that for Rb the $F=F^{*},M_{F}=F^{*}$ molecular states can be
generally further characterized by the $f_{a},f_{b}$ quantum numbers. This is
a consequence of the small $\xi_{\rm ex}\ll 1$ parameter [see Figs. 2(a) and
2(b)], indicating that for these molecular states hyperfine interaction
effectively dominates over exchange interaction (see supplemental material
[53]). There is a propensity to conserve the $f_{a},f_{b}$ quantum numbers.
This is because the recombination primarily occurs mechanically and requires
no spin flips [12]. Indeed, both our numerical calculations [Figs. 2(a) and
2(b)] and experimental observations in Refs. [11, 12, 13] have already shown
this. If one or both atomic hyperfine spins flip the rate typically drops by a
factor $\approx 20$, indicate the difficulty of flipping atomic hyperfine
spins.

For cases of intermediate $\xi_{\rm ex}\sim 1$ [as those of 133Cs, 41K, 23Na
and 39K in Figs. 2(c)-2(f)] electronic exchange and hyperfine interactions are
comparable and the $F=F^{*},M_{F}=F^{*}$ molecular states are typically not
well-described by the $f_{a},f_{b}$ quantum numbers. Their spin states are
mixed, $|FM_{F}(\chi_{1})\rangle$ or $|FM_{F}(\chi_{2})\rangle$ (see [64] for
a detailed description). As a consequence, a clear separation into preferred
spin states does not exist.

For systems where $\xi_{\rm ex}\gg 1$, as is the case for 7Li, the electronic
exchange interaction dominates over the hyperfine interaction. Most molecular
states are then well described in the $|FM_{F}[SI]]\rangle$ basis where $S$
and $I$ are the total electronic and nuclear spin quantum numbers,
respectively. The data for lithium are shown in Figs. 2(g) and Fig. 3 using
the basis sets $|FM_{F}(\chi)\rangle$ and $|FM_{F}[SI]]\rangle$, respectively.

Figure 3:  Molecular product-state distribution for three-body recombination
of ultracold 7Li atoms with molecular levels characterized in the
$|FM_{F}[SI]\rangle$ basis. The straight solid lines denote the $1/E_{\beta}$
scaling of $L_{3\beta}$ with prefactors in the ratio of 0.7875 : 0.1875 :
0.0250.

We find that the molecular product rates depend upon the overlap between the
initial and final spins, i.e., $|\langle
f_{*}m^{*}_{f},f_{*}m^{*}_{f}|FM_{F}[SI]\rangle|^{2}$. The initial atomic
scattering state $|1\text{-}1\rangle|1\text{-}1\rangle$ has a non-zero overlap
only with the molecular states $|2\text{-}2[02]\rangle$,
$|2\text{-}2[11]\rangle$, and $|2\text{-}2[13]\rangle$, with $|\langle
f_{*}m^{*}_{f},f_{*}m^{*}_{f}|F_{*}M_{F}^{*}[SI]\rangle|^{2}=$ 0.1875, 0.0250
and 0.7875, respectively [54]. According to this analysis the dominant
molecular states should be $|2\text{-}2[13]\rangle$ and
$|2\text{-}2[02]\rangle$, with $|2\text{-}2[11]\rangle$ and other spin states
be only weakly populated. This is, in fact, consistent with our numerical
calculations in Fig. 3.

To further test our interpretation of the spin dependence of the product-state
distribution in terms of $\xi_{\rm ex}$, we have performed additional
calculations in which we artificially change the hyperfine splitting for 87Rb
atoms so that the corresponding $\xi_{\rm ex}$ reaches a value comparable to
that of 23Na and 7Li. We find that the product-state distribution of three-
body recombination of 87Rb changes drastically and, in fact, recovers the
hierarchical structure of other species with similar values of $\xi_{\rm ex}$
[53].

As previously mentioned, the rate of formation for non-conserving
$(F,M_{F})\neq(F_{*},M_{F}^{*})$ molecular states, as shown in Fig. 2,
exhibits a weak dependence on $E_{\beta}$. This behavior contrasts sharply
with that of the conserving $(F,M_{F})$ molecular states discussed so far. In
fact, at low $B$ fields, the formation of non-conserving ($F,M_{F}$) molecular
states in our model can only occur via a process in which the spins of the
atoms in the molecule are changed via the spin-exchange interaction with the
third atom. This is so because the spin-exchange among the atoms within the
molecule cannot change ($F_{*},M_{F}^{*}$). A possible process to form non-
conserving ($F,M_{F}$) molecular states can be through the formation of a
strongly coupled intermediate that conserves ($F,M_{F}$) but flips the atomic
spins ($f_{a}$ and/or $f_{b}$). This molecular state is energetically close to
the final non-conserving $(F,M_{F})$ molecular state, formed subsequently by
interacting with the third atom. For 87Rb atoms [Fig. 2(a)], for instance, the
formation of non-conserving $|3\text{-}2(21)\rangle$ molecules could be
accomplished via the intermediate state $|2\text{-}2(21)\rangle$ while the
formation of $|4\text{-}2(22)\rangle$ molecules would be more favorable via
the $|2\text{-}2(22)\rangle$ state. (Similar processes can also be
hypothesized for other atomic species in Fig. 2.) A more detailed analysis of
the formation of non-conserving ($F,M_{F}$) molecular states is beyond the
scope of the present work. Regardless of the details of how such processes
occur they should become more prevalent as electronic exchange interactions
become more important. In fact, from Fig. 2 it is clear that the increase of
the formation rates of such molecular states is consistent with the increase
of $\xi_{\rm ex}$ across the different atomic species.

In summary, our in-depth analysis of the spin product-state distribution of
vdW molecules reveals a two-level hierarchy of spin propensity rules in
ultracold three-body recombination. The spin hierarchy, characterizing the
propensity of changes in various spins of reactant atomic pairs during
molecule formation, strongly depends on the interplay between effective
electronic exchange and hyperfine interactions across alkali species. The
demonstrated spin sensitivity of three-body recombination can potentially
allow for the coherent control of this reaction via spin mixing [14]. Such
ability can, for instance, allow for the enhancement or suppression the
population of a particular family of molecular products based on the spin
propensity rules we uncovered in our study.

This work was supported by the Baden-Württemberg Stiftung through the
Internationale Spitzenforschung program (BWST, contract No. ISF2017-061) and
by the German Research Foundation (DFG, Deutsche Forschungsgemeinschaft,
contract No. 399903135). We acknowledge support from bwForCluster JUSTUS 2 for
high performance computing. J.P.D. also acknowledges partial support from the
U.S. National Science Foundation (PHY-2012125 and PHY-2308791) and NASA/JPL
(1502690).

##  I Supplemental Material

###  I.1 Three-body hyperspherical representation

In the hyperspherical adiabatic representation [55, 56] one of the major tasks
is to solve the hyperangular adiabatic equation $\hat{H}_{\rm
ad}\Phi_{\nu}(R;\Omega)=U_{\nu}(R)\Phi_{\nu}(R;\Omega)$, at fixed values of
$R$, in order to determine the three-body potentials $U_{\nu}$ and channel
functions $\Phi_{\nu}$, both of which are required for the study of the
solutions of Eq. (1). The adiabatic Hamiltonian $\hat{H}_{\rm ad}$,

$\displaystyle\hat{H}_{\rm ad}=\frac{\hat{\Lambda}^{2}(\Omega)+15/4}{2\mu
R^{2}}\hbar^{2}+\hat{V}_{T}(R,\Omega)+\hat{H}_{\rm hf}(B),$ (3)

contains the hyperangular kinetic energy via the hyperangular momentum
operator [55, 56], $\hat{\Lambda}$, the $B$-field dependent hyperfine atomic
energies, $\hat{H}_{\rm hf}$, as well as all the interatomic interactions of
the system, $\hat{V}_{T}$. In the present study we assume that the interaction
is given by

$\displaystyle\hat{V}_{T}(R,\Omega)=\hat{V}(r_{12})+\hat{V}(r_{23})+\hat{V}(r_{31}),$
(4)

where $r_{ij}$ is the distance between atoms $i$ and $j$. The interatomic
interaction, $\hat{V}$, is represented by

$\displaystyle\hat{V}(r_{ij})=\sum_{SM_{S}}|SM_{S}\rangle V_{S}(r_{ij})\langle
SM_{S}|,$ (5)

where $S$ is the total electronic spin and $M_{S}$ is azimuthal projection. In
our present study, we used realistic Bohr-Oppenheimer potentials for the
singlet ($S=0$) and triplet ($S=1$) interactions with an additional term to
control the number of bound states in the problem [54].

The solutions of the hyperangular adiabatic equation are obtained by expanding
the channel functions $\Phi_{\nu}$ on the basis of separated atoms hyperfine
spins states,
$|\sigma\rangle\equiv|f_{i}m_{f_{i}}\rangle|f_{j}m_{f_{j}}\rangle|f_{k}m_{f_{k}}\rangle$,

$\displaystyle\Phi_{\nu}(R;\Omega)=\sum_{\sigma}\phi^{\sigma}_{\nu}(R;\Omega)|\sigma\rangle.$
(6)

In practice we restrict the spin basis $|\sigma\rangle$ to those where two of
the atoms are found in allowed spins states for a given
$M_{F}=m_{f_{i}}+m_{f_{j}}$ while the third atom remains in the income spin
state. Applying this expansion to the hyperangular adiabatic equation results
in a coupled system of equations for the components of $\phi_{\nu}^{\sigma}$:

$\displaystyle\left[\frac{\hat{\Lambda}^{2}(\Omega)+15/4}{2\mu
R^{2}}\hbar^{2}+E^{\sigma}_{\rm
hf}(B)-U_{\nu}(R)\right]\phi^{\sigma}_{\nu}(R;\Omega)$
$\displaystyle+\sum_{\sigma^{\prime}}V_{T}^{\sigma\sigma^{\prime}}(R,\Omega)\phi^{\sigma^{\prime}}_{\nu}(R;\Omega)=0,$
(7)

where $E_{\rm hf}^{\sigma}$ is the sum of the hyperfine energies of the three
separated atoms in the $|\sigma\rangle$ spin state at the magnetic field $B$.

###  I.2 The exchange parameter $\xi_{\rm ex}$

In this section we briefly discuss the physics encapsulated in the definition
of the dimensionless parameter $\xi_{\rm ex}$, presented in more details in
Ref. [54]. We define it as the ratio of the effective electronic exchange
energy, $\tilde{E}_{\rm ex}$, to the effective hyperfine, $\tilde{E}_{\rm
hf}$, interaction,

$\xi_{\rm ex}\equiv\frac{\tilde{E}_{\rm ex}}{\tilde{E}_{\rm hf}},$ (8)

with $\tilde{E}_{\rm ex}$ and $\tilde{E}_{\rm hf}$ discussed below. This
definition of $\xi_{\rm ex}$ is inspired by the parameter $\rho$ introduced in
Ref. [12] and provides a simple way to characterize the spin properties of vdW
molecules across different species [54].

In order to estimate the characteristic electronic exchange interaction we
make use of the universality of the dimer spectrum for vdW interactions [68,
40]. The dimer levels universally span a specific energy interval for each
vibrational quantum number $v$ [68, 40]. The $v=-1$ level must be in the
energy range [-39.5$E_{\rm vdW}$, 0], the $v=-2$ levels in [-272.5$E_{\rm
vdW}$, -39.5$E_{\rm vdW}$], and so on, regardless the short-range detail of
the potentials. These universal energy intervals are known as vdW energy bins
[68]. The short-range detail determining the actual location of dimer levels
in vdW energy bins depends on the scattering length $a$ and can be
parameterized by the quantity $u$ defined as [69]

$u(a)=\tan^{-1}[\bar{a}/(a-\bar{a})]/\pi,$ (9)

where $\bar{a}\approx 0.9560$ $r_{\rm vdW}$. For instance, $u=0$ (or
$a=\pm\infty$) indicates that the levels appear at the boundaries of the
energy bins. Note that as $a$ changes and the dimer energy level moves with
the bin the value of $u$ varies within the range [-0.5, 0.5]. For alkali
atoms, where the interaction is characterized by the singlet, $a_{s}$, and
triplet, $a_{t}$, scattering lengths, the quantity $|u(a_{s})-u(a_{t})|$ will
describe the relative position of a pair of singlet and triplet levels in each
energy bin, varying between 0 and 1. In order consider a more general case
where the nearest singlet and triplet level can belong to the same or adjacent
energy bins we define $u_{st}=\min[|u(a_{s})-u(a_{t})|,1-|u(a_{s})-u(a_{t})|]$
to characterize the relative position of such states. Here, $u_{st}$ varies
from 0 to 1/2, where whenever $u_{st}\approx 0.5$, the nearest singlet and
triplet energy levels are well separated indicating that the electronic
exchange interaction between these levels is strong. For $u_{st}\ll 1/2$,
since energy difference between the singlet and triplet energy levels is
small, the electronic exchange interaction is weak between these levels.
According to our discussion above, the quantity

$\tilde{E}_{\rm ex}=u_{st}E_{\rm bin}$ (10)

will characterize the energy shift between the nearby singlet and triplet
levels caused by the electronic exchange interaction.

Similarly, we can characterize the hyperfine interaction by the energy
difference between different hyperfine molecular states. In general, the
hyperfine interaction shifts the dimer levels by the atomic hyperfine
splitting constant $E_{\rm hf}$ or 2$E_{\rm hf}$. Nevertheless, $E_{\rm hf}$
will not properly describe the energy difference of nearby hyperfine molecular
levels when $E_{\rm hf}\gtrsim E_{\rm bin}$. Instead, this energy difference
should be constrained by the size of the corresponding energy bin, typically
with an upper limit around $E_{\rm bin}/2$ [54]. Therefore, we use

$\tilde{E}_{\rm hf}=\min[E_{\rm hf},E_{\rm{bin}}/2]$ (11)

to characterize the effective hyperfine interaction. In fact, the definition
above is crucial for some species, for instance 133Cs, where the hyperfine
splitting is very large ($E_{\rm hf}\approx 3455E_{\rm vdW}$).

For the discussion in the present work, we use the value of $\xi_{\rm ex}$
with $E_{\rm bin}=39.5$ $E_{\rm vdW}$ of $v=-1$. This choice emphasizes the
nature of the most weakly bound molecular states in the system. In fact,
taking the dependence of $E_{\rm bin}$ on $v$ into account for the molecular
states we considered in our study, $\xi_{\rm ex}$ remains rather constant for
85Rb, 87Rb and 133Cs but increase with the increasing $v$ for 7Li, 23Na, 39K
and 41K [54]. As a result, our classification of 85Rb and 87Rb as atomic
species with weak electronic exchange ($\xi_{\rm ex}\ll 1$), and 7Li as a
strong electronic exchange species ($\xi_{\rm ex}\gg 1$) and 133Cs as
intermediate electronic exchange species ($\xi_{\rm ex}\lesssim 1$) remain
unchanged. The electronic exchange interaction is, however, underestimated by
the presented $\xi_{\rm ex}$ of $v=-1$ for 23Na, 39K and 41K. In fact, most
more deeply bound (if not all) 23Na2, 39K2 and 41K2 molecules are
characterized by $\xi_{\rm ex}\gg 1$ [54].

Figure 4:  Product-state distribution for three-body recombination of 87Rb
with artificially suppressed atomic hyperfine splitting. The resulting
$\xi_{\rm ex}$ has the same value as that of (a) 23Na and (b) 7Li (b). The
straight solid lines denote the $1/E_{\beta}$ scaling of $L_{3\beta}$, while
the dotted lines indicate a constant rate. The shadowed area indicates the
typical range of scatter of the individual partial rates around the
$1/E_{\beta}$ scaling lines. (b) shares the legend with (a).

To further demonstrate the predictive power of $\xi_{\rm ex}$, we artificially
change the hyperfine splitting of 87Rb while keeping the singlet and triplet
scattering lengths fixed so that the new value for $\xi_{\rm ex}$ is similar
to that for 23Na or 7Li. For the relationship between $\xi_{\rm ex}$ and the
atomic hyperfine splitting and the singlet and triplet scattering lengths, we
refer to Ref. [54]. We find that such an artificial change in 87Rb drastically
changes the product state distribution of recombination originally displayed
in the main text [Fig. 2(a)]. Despite the difference in total recombination
rate $L_{3}$, the new 87Rb product state distribution resembles that of the
real 23Na [Fig. 2(e)] or real 7Li [Fig. 2(g)], as is shown in Figs. 4(a) and
4(b), respectively. In particular, the rates for the molecular products in the
$|2\text{-}2(\chi_{1})\rangle$ state are more prominent than those for
$|2\text{-}2(\chi_{2})\rangle$ in Fig. 4(b) [similar to Fig. 2(g)] while they
are comparable in Fig. 4(a) [similar to Fig. 2(e)]. Nevertheless, a
perceivable discrepancy in the result of non-conserving ($F,M_{F}$) molecular
states between Fig. 4(a) and Fig. 2(e) indicates that the formation of these
states is more sensitive to the system’s interactions.

###  I.3 Three-body calculations beyond $M_{F}=M^{*}_{F}$

Figure 5:  Product-state of three-body recombination for 87Rb (a) and 7Li (b)
from numerical simulations including multiple spin states of the third atom.
The straight solid lines denote the $1/E_{\beta}$ scaling of $L_{3\beta}$,
while the dotted lines indicate a constant rate.The shadowed area indicates
the typical range of scatter of the individual partial rates around the
$1/E_{\beta}$ scaling lines.

Here we justify the validity of the approximation of fixing the third atom’s
spin state, which has been employed in most of our three-body calculations.
For that, we perform additional three-body calculations by expanding the spin
space of the third atom to
$\\{|f_{*}m^{*}_{f}\rangle,|(f_{*}+1)(m^{*}_{f}-1)\rangle\\}$ so that
molecular products with $M_{F}=M^{*}_{F}+1$ will be taken into account. Our
results for 87Rb and 7Li using this expanded basis are shown in Fig. 5. The
calculated $L_{3}=1.88\times 10^{-29}$ cm6/s of 87Rb deviates from the
previous value $L_{3}=1.25\times 10^{-29}$ cm6/s. We attribute this
discrepancy to numerical fluctuations, which have negligible influence on
$L_{3\beta}/L_{3}$. For 87Rb, the second minor discrepancy compared to the
previous result [in Fig. 2(a)] is that the recombination rates of non-
conserving ($F,M_{F}$) molecular products are more concentrated in Fig. 5(a).
This indicates again the formation of non-conserving ($F,M_{F}$) molecular
states is sensitive to the system’s interactions. Nevertheless, we find in
both atomic species the main conclusions drawn from the previous single
$M_{F}$ basis simulations do not change and the molecules with
$M_{F}=M^{*}_{F}+1$ are weakly populated. The total reaction fraction to these
molecules is below one percent for 7Li and on the order of 10${}^{\text{-}4}$
for 87Rb. The low reaction rates into the molecular states with $M_{F}\neq
M^{*}_{F}$ indicates that the third atom in general interacts mechanically
when the other two atoms form a molecule. As discussed in three-body
recombination of 85Rb [12], the third atom should come into proximity (at a
typical internuclear separation the classical out-turning point of the
molecular state) with the two atoms forming a bond, such that it facilitates
the release of excess energy during molecular formation. Since we consider
only weekly bound vdW molecules in the present study, the spin-independent vdW
interaction is more prominent than the electronic exchange interaction at such
internuclear separations [54]. As a result, when the third atom interacts, the
mechanical contribution should be more dominant than the spin-flip
contribution.

## References

  * Collins [2002] M. A. Collins, Molecular potential-energy surfaces for chemical reaction dynamics, Theoretical Chemistry Accounts 108, 313 (2002).
  * Truhlar [2007] D. G. Truhlar, Valence bond theory for chemical dynamics, Journal of Computational Chemistry 28, 73 (2007).
  * Qu  _et al._ [2018] C. Qu, Q. Yu, and J. M. Bowman, Permutationally invariant potential energy surfaces, Annual Review of Physical Chemistry 69, 151 (2018).
  * Berry [1966] R. S. Berry, Ionization of Molecules at Low Energies, The Journal of Chemical Physics 45, 1228 (1966).
  * Moore [1973] J. H. Moore, Investigation of the Wigner Spin Rule in Collisions of ${\mathrm{N}}^{+}$ with He, Ne, Ar, ${\mathrm{N}}_{2}$, and ${\mathrm{O}}_{2}$, Phys. Rev. A 8, 2359 (1973).
  * Fano [1985] U. Fano, Propensity rules: An analytical approach, Phys. Rev. A 32, 617 (1985).
  * Busto  _et al._ [2019] D. Busto, J. Vinbladh, S. Zhong, M. Isinger, S. Nandi, S. Maclot, P. Johnsson, M. Gisselbrecht, A. L’Huillier, E. Lindroth, and J. M. Dahlström, Fano’s Propensity Rule in Angle-Resolved Attosecond Pump-Probe Photoionization, Phys. Rev. Lett. 123, 133201 (2019).
  * Lee  _et al._ [1991] A. Lee, C. Enos, and A. Brenton, Collisional excitation of CO: a study of the Wigner spin rule, Int. J. Mass Spectrom. Ion Processes 104, 49 (1991).
  * Hu  _et al._ [2021] M.-G. Hu, Y. Liu, M. A. Nichols, L. Zhu, G. Quéméner, O. Dulieu, and K.-K. Ni, Nuclear spin conservation enables state-to-state control of ultracold molecular reactions, Nature Chemistry 13, 435 (2021).
  * Wolf  _et al._ [2017] J. Wolf, M. Deiß, A. Krükow, E. Tiemann, B. P. Ruzic, Y. Wang, J. P. D’Incao, P. S. Julienne, and J. Hecker Denschlag, State-to-state chemistry for three-body recombination in an ultracold rubidium gas, Science 358, 921 (2017).
  * Wolf  _et al._ [2019] J. Wolf, M. Deiß, and J. Hecker Denschlag, Hyperfine Magnetic Substate Resolved State-to-State Chemistry, Phys. Rev. Lett. 123, 253401 (2019).
  * Haze  _et al._ [2022] S. Haze, J. P. D’Incao, D. Dorer, M. Deiß, E. Tiemann, P. S. Julienne, and J. Hecker Denschlag, Spin-Conservation Propensity Rule for Three-Body Recombination of Ultracold Rb Atoms, Phys. Rev. Lett. 128, 133401 (2022).
  * Haze  _et al._ [2023] S. Haze, J. P. D’Incao, D. Dorer, J.-L. Li, M. Deiß, E. Tiemann, P. S. Julienne, and J. Hecker Denschlag, Energy scaling of the product state distribution for three-body recombination of ultracold atoms, Phys. Rev. Res. 5, 013161 (2023).
  * [14] S. Haze, J.-L. Li, D. Dorer, J. P. D’Incao, P. S. Julienne, E. Tiemann, M. Deiß, J. Hecker Denschlag, Controlling the few-body reaction pathway using a Feshbach resonance, submitted.
  * Hermsmeier  _et al._ [2021] R. Hermsmeier, J. Kłos, S. Kotochigova, and T. V. Tscherbul, Quantum Spin State Selectivity and Magnetic Tuning of Ultracold Chemical Reactions of Triplet Alkali-Metal Dimers with Alkali-Metal Atoms, Phys. Rev. Lett. 127, 103402 (2021).
  * Tscherbul and Krems [2006] T. V. Tscherbul and R. V. Krems, Controlling Electronic Spin Relaxation of Cold Molecules with Electric Fields, Phys. Rev. Lett. 97, 083201 (2006).
  * Yang  _et al._ [2019] H. Yang, D.-C. Zhang, L. Liu, Y.-X. Liu, J. Nan, B. Zhao, and J.-W. Pan, Observation of magnetically tunable Feshbach resonances in ultracold 23Na40K + 40K collisions, Science 363, 261 (2019).
  * Bigagli  _et al._ [2023] N. Bigagli, C. Warner, W. Yuan, S. Zhang, I. Stevenson, T. Karman, and S. Will, Collisionally stable gas of bosonic dipolar ground-state molecules, Nature Physics 19, 1579 (2023).
  * Son  _et al._ [2022] H. Son, J. J. Park, Y.-K. Lu, A. O. Jamison, T. Karman, and W. Ketterle, Control of reactive collisions by quantum interference, Science 375, 1006 (2022).
  * Chen  _et al._ [2023] X.-Y. Chen, A. Schindewolf, S. Eppelt, R. Bause, M. Duda, S. Biswas, T. Karman, T. Hilker, I. Bloch, and X.-Y. Luo, Field-linked resonances of polar molecules, Nature 614, 59 (2023).
  * Liu  _et al._ [2024] Y.-X. Liu, L. Zhu, J. Luke, J. J. A. Houwman, M. C. Babin, M.-G. Hu, and K.-K. Ni, Quantum interference in atom-exchange reactions, Science 384, 1117 (2024).
  * Luke  _et al._ [2024] J. Luke, L. Zhu, Y.-X. Liu, and K.-K. Ni, Reaction interferometry with ultracold molecules, Faraday Discuss. , (2024).
  * Anderson  _et al._ [1995] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Observation of Bose-Einstein condensation in dilute atomic vapor, Science 269, 198 (1995).
  * Davis  _et al._ [1995] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Bose-Einstein condensation in a gas of sodium atoms, Phys. Rev. Lett. 75, 3969 (1995).
  * Bradley  _et al._ [1997] C. C. Bradley, C. A. Sackett, and R. G. Hulet, Bose-Einstein condensation of lithium: Observation of limited condensate number, Phys. Rev. Lett. 78, 985 (1997).
  * Burt  _et al._ [1997] E. A. Burt, R. W. Ghrist, C. J. Myatt, M. J. Holland, E. A. Cornell, and C. E. Wieman, Coherence, Correlations, and Collisions: What One Learns about Bose-Einstein Condensates from Their Decay, Phys. Rev. Lett. 79, 337 (1997).
  * Inouye  _et al._ [1998] S. Inouye, M. Andrews, J. Stenger, H. Miesner, D. Stamper-Kurn, and W. Ketterle, Observation of Feshbach resonances in a Bose–Einstein condensate, Nature 392, 151 (1998).
  * Courteille  _et al._ [1998] P. Courteille, R. S. Freeland, D. J. Heinzen, F. A. van Abeelen, and B. J. Verhaar, Observation of a Feshbach resonance in cold atom scattering, Phys. Rev. Lett. 81, 69 (1998).
  * Stenger  _et al._ [1999] J. Stenger, S. Inouye, M. R. Andrews, H.-J. Miesner, D. M. Stamper-Kurn, and W. Ketterle, Strongly enhanced inelastic collisions in a Bose-Einstein condensate near Feshbach resonances, Phys. Rev. Lett. 82, 2422 (1999).
  * Roberts  _et al._ [2000] J. L. Roberts, N. R. Claussen, S. L. Cornish, and C. E. Wieman, Magnetic Field Dependence of Ultracold Inelastic Collisions near a Feshbach Resonance, Phys. Rev. Lett. 85, 728 (2000).
  * Marte  _et al._ [2002] A. Marte, T. Volz, J. Schuster, S. Dürr, G. Rempe, E. G. M. van Kempen, and B. J. Verhaar, Feshbach resonances in rubidium 87: Precision measurement and analysis, Phys. Rev. Lett. 89, 283202 (2002).
  * Weber  _et al._ [2003a] T. Weber, J. Herbig, M. Mark, H.-C. Nägerl, and R. Grimm, Bose-Einstein condensation of cesium, Science 299, 232 (2003a).
  * Weber  _et al._ [2003b] T. Weber, J. Herbig, M. Mark, H.-C. Nägerl, and R. Grimm, Three-body recombination at large scattering lengths in an ultracold atomic gas, Phys. Rev. Lett. 91, 123201 (2003b).
  * Braaten and Hammer [2006] E. Braaten and H.-W. Hammer, Universality in few-body systems with large scattering length, Physics Reports 428, 259 (2006).
  * Naidon and Endo [2017] P. Naidon and S. Endo, Efimov physics: a review, Reports on Progress in Physics 80, 056001 (2017).
  * Greene  _et al._ [2017] C. H. Greene, P. Giannakeas, and J. Pérez-Ríos, Universal few-body physics and cluster formation, Rev. Mod. Phys. 89, 035006 (2017).
  * D’Incao [2018] J. P. D’Incao, Few-body physics in resonantly interacting ultracold quantum gases, Journal of Physics B 51, 043001 (2018).
  * Howard [1993] B. J. Howard, The structure and dynamics of van der Waals molecules, in  _Structures and Conformations of Non-Rigid Molecules_, edited by J. Laane, M. Dakkouri, B. van der Veken, and H. Oberhammer (Springer Netherlands, Dordrecht, 1993) pp. 137–161.
  * Gao [1998] B. Gao, Solutions of the Schrödinger equation for an attractive ${1/r}^{6}$ potential, Phys. Rev. A 58, 1728 (1998).
  * Gao [2000] B. Gao, Zero-energy bound or quasibound states and their implications for diatomic systems with an asymptotic van der Waals interaction, Phys. Rev. A 62, 050702 (2000).
  * Gao [2004] B. Gao, Binding energy and scattering length for diatomic systems, Journal of Physics B 37, 4273 (2004).
  * Blaney and Ewing [1976] B. L. Blaney and G. E. Ewing, Van der Waals molecules, Annual Review of Physical Chemistry 27, 553 (1976).
  * Jeziorski  _et al._ [1994] B. Jeziorski, R. Moszynski, and K. Szalewicz, Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes, Chemical Reviews 94, 1887 (1994).
  * Wormer and van der Avoird [2000] P. E. S. Wormer and A. van der Avoird, Intermolecular Potentials, Internal Motions, and Spectra of van der Waals and Hydrogen-Bonded Complexes, Chemical Reviews 100, 4109 (2000), pMID: 11749342.
  * Koperski [2002] J. Koperski, Study of diatomic van der Waals complexes in supersonic beams, Physics Reports 369, 177 (2002).
  * Reilly and Tkatchenko [2015] A. M. Reilly and A. Tkatchenko, van der Waals dispersion interactions in molecular materials: beyond pairwise additivity, Chem. Sci. 6, 3289 (2015).
  * Hermann  _et al._ [2017] J. Hermann, R. A. J. DiStasio, and A. Tkatchenko, First-principles models for van der Waals interactions in molecules and materials: Concepts, theory, and applications, Chemical Reviews 117, 4714 (2017), pMID: 28272886\.
  * Brahms  _et al._ [2010] N. Brahms, T. V. Tscherbul, P. Zhang, J. Kłos, H. R. Sadeghpour, A. Dalgarno, J. M. Doyle, and T. G. Walker, Formation of van der Waals Molecules in Buffer-Gas-Cooled Magnetic Traps, Phys. Rev. Lett. 105, 033001 (2010).
  * Brahms  _et al._ [2011] N. Brahms, T. V. Tscherbul, P. Zhang, J. Kłos, R. C. Forrey, Y. S. Au, H. R. Sadeghpour, A. Dalgarno, J. M. Doyle, and T. G. Walker, Formation and dynamics of van der Waals molecules in buffer-gas traps, Phys. Chem. Chem. Phys. 13, 19125 (2011).
  * Quiros  _et al._ [2017] N. Quiros, N. Tariq, T. V. Tscherbul, J. Kłos, and J. D. Weinstein, Cold Anisotropically Interacting van der Waals Molecule: TiHe, Phys. Rev. Lett. 118, 213401 (2017).
  * Mirahmadi and Pérez-Ríos [2021a] M. Mirahmadi and J. Pérez-Ríos, On the formation of van der Waals complexes through three-body recombination, The Journal of Chemical Physics 154, 034305 (2021a).
  * Mirahmadi and Pérez-Ríos [2021b] M. Mirahmadi and J. Pérez-Ríos, Classical threshold law for the formation of van der Waals molecules, The Journal of Chemical Physics 155, 094306 (2021b).
  * [53] See Supplemental Material for details of our model and calculations, which includes Refs. [68, 69].
  * Li [2] J.-L. Li, P. S. Julienne, J. Hecker Denschlag, and J. P. D’Incao, Spin structure of diatomic van der Waals molecules of alkali atoms, in preparation.
  * Suno  _et al._ [2002] H. Suno, B. D. Esry, C. H. Greene, and J. P. Burke, Three-body recombination of cold helium atoms, Phys. Rev. A 65, 042725 (2002).
  * Wang  _et al._ [2011] J. Wang, J. P. D’Incao, and C. H. Greene, Numerical study of three-body recombination for systems with many bound states, Phys. Rev. A 84, 052721 (2011).
  * Julienne and Hutson [2014] P. S. Julienne and J. M. Hutson, Contrasting the wide Feshbach resonances in ${}^{6}\mathrm{Li}$ and ${}^{7}\mathrm{Li}$, Phys. Rev. A 89, 052715 (2014).
  * Knoop  _et al._ [2011] S. Knoop, T. Schuster, R. Scelle, A. Trautmann, J. Appmeier, M. K. Oberthaler, E. Tiesinga, and E. Tiemann, Feshbach spectroscopy and analysis of the interaction potentials of ultracold sodium, Phys. Rev. A 83, 042704 (2011).
  * Tiemann  _et al._ [2020] E. Tiemann, P. Gersema, K. K. Voges, T. Hartmann, A. Zenesini, and S. Ospelkaus, Beyond Born-Oppenheimer approximation in ultracold atomic collisions, Phys. Rev. Res. 2, 013366 (2020).
  * Strauss  _et al._ [2010] C. Strauss, T. Takekoshi, F. Lang, K. Winkler, R. Grimm, J. Hecker Denschlag, and E. Tiemann, Hyperfine, rotational, and vibrational structure of the ${a}^{3}{\Sigma}_{u}^{+}$ state of ${}^{87}\mathrm{Rb}$2, Phys. Rev. A 82, 052514 (2010).
  * Berninger  _et al._ [2013] M. Berninger, A. Zenesini, B. Huang, W. Harm, H.-C. Nägerl, F. Ferlaino, R. Grimm, P. S. Julienne, and J. M. Hutson, Feshbach resonances, weakly bound molecular states, and coupled-channel potentials for cesium at high magnetic fields, Phys. Rev. A 87, 032517 (2013).
  * Chapurin  _et al._ [2019] R. Chapurin, X. Xie, M. J. Van de Graaff, J. S. Popowski, J. P. D’Incao, P. S. Julienne, J. Ye, and E. A. Cornell, Precision test of the limits to universality in few-body physics, Phys. Rev. Lett. 123, 233402 (2019).
  * Xie  _et al._ [2020] X. Xie, M. J. Van de Graaff, R. Chapurin, M. D. Frye, J. M. Hutson, J. P. D’Incao, P. S. Julienne, J. Ye, and E. A. Cornell, Observation of Efimov Universality across a Nonuniversal Feshbach Resonance in ${}^{39}\mathrm{K}$, Phys. Rev. Lett. 125, 243401 (2020).
  * [64] For molecules with spins states $|\Sigma_{ab}(v,L_{R})\rangle$ having a projection $|\langle FM_{F}(f_{a}f_{b})|\Sigma_{ab}(v,L_{R})\rangle|^{2}>0.9$ we assign the spin state to $|FM_{F}(f_{a}f_{b})\rangle$ state, otherwise, they are assigned to a mixed state $|FM_{F}(\chi)\rangle$. To elaborate, the $|6\text{-}6(\chi_{1})\rangle$ state of 133Cs2 is a mixture of $|6\text{-}6(33)\rangle$, $|6\text{-}6(43)\rangle$, and $|6\text{-}6(44)\rangle$. The $|4\text{-}4(\chi_{1})\rangle$ state of 85Rb2 is a mixture of $|4\text{-}4(22)\rangle$ and $|4\text{-}4(32)\rangle$, or of $|4\text{-}4(22)\rangle$ and $|4\text{-}4(33)\rangle$. The $|2\text{-}2(\chi_{1})\rangle$ state is a mixture of $|2\text{-}2(11)\rangle$ and $|2\text{-}2(22)\rangle$ for 87Rb2, while it is a mixture of $|2\text{-}2(11)\rangle$, $|2\text{-}2(21)\rangle$, and $|2\text{-}2(22)\rangle$ for molecules of other species. For 41K2, 23Na2, 39K2, and 7Li2, the mixed state $|2\text{-}2(\chi_{2})\rangle$ is a mixture of $|2\text{-}2(21)\rangle$ and $|2\text{-}2(22)\rangle$. A mixed state has less than 0.1 component in the unspecified $|FM_{F}(f_{a}f_{b})\rangle$ state.
  * [65] Note that for simplicity our notation of the atom-dimer state, $|\Psi_{FM_{F}}(v,L_{R})\rangle|fm_{f}\rangle$, omits the atom-dimer relative motion and any additional quantum numbers that characterize the molecular eigenstate.
  * Li  _et al._ [2022] J.-L. Li, T. Secker, P. M. A. Mestrom, and S. J. J. M. F. Kokkelmans, Strong spin-exchange recombination of three weakly interacting ${}^{7}\mathrm{Li}$ atoms, Phys. Rev. Res. 4, 023103 (2022).
  * van de Kraats  _et al._ [2024] J. van de Kraats, D. J. M. Ahmed-Braun, J.-L. Li, and S. J. J. M. F. Kokkelmans, Emergent Inflation of the Efimov Spectrum under Three-Body Spin-Exchange Interactions, Phys. Rev. Lett. 132, 133402 (2024).
  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010).
  * Friedrich and Trost [2004] H. Friedrich and J. Trost, Working with WKB waves far from the semiclassical limit, Physics Reports 397, 359 (2004).


11institutetext: NeoInstinct S.A.
Rue traversière 3, 1018 Lausanne (Switzerland)
11email<EMAIL_ADDRESS>
# A unified view on Self-Organizing Maps (SOMs)
and Stochastic Neighbor Embedding (SNE)

Thibaut Kulak  Anthony Fillion  François Blayo

###### Abstract

We propose a unified view on two widely used data visualization techniques:
Self-Organizing Maps (SOMs) and Stochastic Neighbor Embedding (SNE). We show
that they can both be derived from a common mathematical framework. Leveraging
this formulation, we propose to compare SOM and SNE quantitatively on two
datasets, and discuss possible avenues for future work to take advantage of
both approaches.

###### Keywords:

Data visualization Stochastic Neighbor Embedding Self-Organizing Map Auto
organization Representation learning.

##  1 Introduction

Self-organizing Maps (SOMs) [5] have been widely used over the last decades on
diverse applications [11], and have emerged as a very efficient biologically-
inspired method for learning organized representations of data. More recently,
the Stochastic Neighbor Embedding (SNE) algorithm [4] has been proposed to
visualize a set of datapoints in a two-dimensional space, hence creating a
map. This algorithm, and notably its extension t-distributed Stochastic
Neighbor Embedding (t-SNE) [6] has become very popular for visualizing large
datasets. The goal of this work is to propose a unified view over those two
algorithms.

We will show that both SOM and SNE algorithms can be seen in the context of a
general mathematical formulation that attempts to find points in a 2D space
associated to weights in an observation space, where the weights represent the
data, and where neighborhood relations in the 2D space and the observation
space are respected (points that are close in the map should be close in the
observation space, and vice versa). In the context of this optimization
problem, we will see that SOM’s algorithm fixes the points in the 2D space,
and then optimizes the weights of the neurons, while SNE fixes the weights of
the neurons, and optimizes the positions of the points in the 2D space.

Such a unified view of those algorithms brings several contributions: for
SOMs, this sheds light on what is actually minimized and on how we might be
able to assess the quality of a map both in terms of organization and
representation, which remains so far non-trivial [2, 1]. For SNE, we propose a
formulation that showcases the fact that SNE is not only a visualization
technique but is also learning a representation of the data. Moreover, the
number of neuron weights can be different from the number of datapoints in our
formulation, which goes beyond the standard SNE approach that assumes one
neuron per datapoint. Most importantly, as we will discuss later, we hope that
bridging those algorithms will open the way to future works leveraging their
respective advantages.

Our paper is organized as follows: first, we introduce the mathematical
framework, and show that the SOM and SNE algorithms can be derived from it.
Then, we exploit the fact that both algorithms pursue a very similar objective
to compare their results quantitatively on two different datasets. Finally, we
conclude and discuss interesting avenues for future work leveraging our
unified formulation.

###  1.1 Notations and definitions

We assume there are $N$ points $\\{z_{i}\\}_{i=1}^{N}$ in a 2D space
(visualization space), corresponding to $N$ weights (neuron weights)
$\\{w_{i}\\}_{i=1}^{N}$ in a $d$-dimensional space (observation space). We
introduce two conditional distributions $p(i|x)$ and $q(i|x)$ defined for a
datapoint (stimulus) $x$. Those distributions are defined over the discrete
set $\llbracket 1;N\rrbracket$:

$p(i|x)=\frac{\exp{\Big{(}-\frac{||x-w_{i}||^{2}}{\sigma^{2}}}\Big{)}}{\sum_{j=1}^{N}\exp{\Big{(}-\frac{||x-w_{j}||^{2}}{\sigma^{2}}}\Big{)}}\quad\text{and}\quad
q(i|x)=\frac{\exp{\Big{(}-||z_{i}-z_{i^{*}(x)}||^{2}}\Big{)}}{\sum_{j=1}^{N}\exp{\Big{(}-||z_{j}-z_{i^{*}(x)}||^{2}\Big{)}}},$
(1)

where $\sigma$ is an hyperparameter111Note that $p(i|x)$ has been obtained
from Bayes rule by choosing $p(x|i)=\mathcal{N}(w_{i},\sigma^{2})$ and a
uniform prior over the discrete space $p(i)=1/N$, so $\sigma$ is actually the
chosen standard deviation of $p(x|i)$. Such model also defines a probability
distribution on the observation space
$p(x)=\sum_{i=1}^{N}p(x|i)p(i)=1/N\sum_{i=1}^{N}p(x|i)$, and
$i^{*}(x)=\arg\max_{i=1}^{N}p(i|x)$ is the index of the winning neuron (the
neuron that responds the most to the stimulus $x$).

We can interpret those conditional distributions qualitatively:

  * •

$p(i|x)$ increases with the proximity of the weights to the stimulus in the
observation space, i.e., which neurons respond to the stimulus.

  * •

$q(i|x)$ increases with the proximity of the points to the stimulus in the 2D
(visualization) space

To give further intuition on those conditional distributions, we show in Fig.
1 those distributions for two different maps, one organized map that was
obtained using SOM’s algorithm and one non-organized map (that was obtained by
randomly permuting the indices of the neurons weights of the organized map).
We can see that when the map is organized $p(i|x)$ and $q(i|x)$ are close to
each other, which is not the case for the non-organized map. More precisely,
if we look at $p(i|x)$ and $q(i|x)$ in the visualization space, $q(i|x)$ looks
like a Gaussian distribution no matter if the map is organized or not, but
$p(i|x)$ looks like a Gaussian only if the map is organized, which indeed
means that points around the winning neuron also respond to the stimulus.
Similarly, if we look at $p(i|x)$ and $q(i|x)$ in the observation space,
$p(i|x)$ looks like a Gaussian distribution no matter if the map is organized
or not, but $q(i|x)$ looks like a Gaussian only if the map is organized, which
indeed means that the points that respond to the stimulus are those that are
close in the visualization space. We will see in the next sections that both
SOMs and SNE methods attempt to make these distributions close to each other,
while representing the data.

Figure 1: Visualization of an unorganized map (left) and an organized map
(right), along with the conditional distributions $p(i|x)$ and $q(i|x)$, both
in visualization space (points $\\{z_{i}\\}$ colored by the probabilities,
middle row) and observation space (weights $\\{w_{i}\\}$ colored by the
probabilities, bottom row). Purple color corresponds to low probabilities and
yellow to high probabilities

###  1.2 SOMs in this framework

Let us assume that stimulus $x$’s are samples from an unknown density
$p^{\text{data}}(x)$. We propose the following objective function, that is a
trade-off between organization (in the sense of making $p(i|x)$ and $q(i|x)$
close) and representation (maximizing $p(x)$):

$J^{\text{SOM}}\Big{(}\\{z_{i}\\}_{i=1}^{N},\\{w_{i}\\}_{i=1}^{N}\Big{)}=\int_{x}\Big{(}\mathcal{D}_{KL}(q(i|x),p(i|x))-\ln
p(x)\Big{)}\;p^{\text{data}}(x)\;dx$ (2)

The standard SOM algorithm (non-batch version) is a stochastic algorithm that
updates the neuron weights as a stimulus $x$ is presented. For a stimulus $x$
coming from $p^{\text{data}}(x)$, we can get the objective function for this
stimulus $x$ by removing the integral over the unknown density function
$p^{\text{data}}(x)$ (stochastic gradient descent):

$J^{\text{SOM}}\Big{(}\\{z_{i}\\}_{i=1}^{N},\\{w_{i}\\}_{i=1}^{N}\;|\;x\Big{)}=\mathcal{D}_{KL}(q(i|x),p(i|x))-\ln
p(x)$ (3)

Differentiating with respect to the neuron weights and negating gives us the
following update rule of the neuron weights:

$-\frac{\partial
J^{\text{SOM}}\Big{(}\\{z_{i}\\}_{i=1}^{N},\\{w_{i}\\}_{i=1}^{N}\;|\;x\Big{)}}{\partial
w_{i}}=\frac{2}{\sigma^{2}}(x-w_{i})q(i|x)$ (4)

We can see that this is exactly the update rule of SOM’s algorithm, where
neuron weights get moved towards the stimulus $x$ if they are in the
neighborhood of the winning neuron.

In SOM’s algorithm, there is typically a neighborhood parameter that is
monotically decreasing. In our framework we did not put a neighborhood
parameter on $q(i|x)$, because we can alternatively see it222Trivially:
$\frac{||z_{i}-z_{i^{*}(x)}||^{2}}{\sigma_{z}^{2}}=||\frac{z_{i}}{\sigma_{z}}-\frac{z_{i^{*}(x)}}{\sigma_{z}}||^{2}$
as a scaling of the points $\\{z_{i}\\}_{i=1}^{N}$ typically chosen on a grid.
Choosing a monotically decreasing neighborhood parameter is equivalent to
scaling monotically the points $\\{z_{i}\\}_{i=1}^{N}$ from a small to a big
grid with fixed neighborhood parameter.

In our view, the SOM algorithm is performing an optimization over a set of
points in a 2D space and their corresponding weights, and is supposing a
particular time-dependent form for the points locations in the 2D space (from
a small to a big grid). Viewing it this way will make it easier to compare
with SNE in the next section, which does not use any neighborhood parameter in
the visualization space.

###  1.3 SNE in this framework

In this section, we will show that the SNE algorithm can be derived from our
framework, hence show that both SOMs and SNE algorithms are solving the same
mathematical problem, and exhibit their differences.

Let us define for a stimulus $x$ the following objective function:

$J^{\text{SNE}}\Big{(}\\{z_{i}\\}_{i=1}^{N},\\{w_{i}\\}_{i=1}^{N}\;|\;x\Big{)}=\mathcal{D}_{KL}(p(i|x),q(i|x))-\ln
p(x)$ (5)

This objective is very similar to Eq. 3, the only difference is the use of the
forward KL divergence instead of the reverse KL divergence. Given that the KL
divergence is not symmetric, Eq. 3 and Eq. 5 differ only by the way they
measure the difference between $p(i|x)$ and $q(i|x)$.

Differentiating with respect to the 2D points locations gives us the following
update rule of the points:

$-\frac{\partial
J^{\text{SNE}}\Big{(}\\{z_{i}\\}_{i=1}^{N},\\{w_{i}\\}_{i=1}^{N}\;|\;x\Big{)}}{\partial
z_{i}}=2(z_{i^{*}(x)}-z_{i})(q(i|x)-p(i|x))$ (6)

This update rule is very similar from the SNE algorithm and can be interpreted
as such: if a neuron weight $w_{i}$ is far from the stimulus $x$ (i.e.,
$p(i|x)$ small), but is close to the winning neuron in the 2D visualization
space (i.e., $q(i|x)$ big) then the point is moved towards the location of the
winning neuron in the 2D space. Reversely, if a neuron weight $w_{i}$ is close
from the stimulus $x$ (i.e., $p(i|x)$ big), but is far from the winning neuron
in the 2D visualization space (i.e., $q(i|x)$ small) then the point is moved
away from the location of the winning neuron in the 2D space.

We have presented here the update rule equations for any set of neuron
weights. In the standard SNE algorithm, usually the neuron weights are equal
to the datapoints, because the number of neurons chosen is equal to the number
of datapoints to visualize. In that case, a stimulus $x$ is therefore a weight
$w_{j}$, and trivially $i^{*}(x)=j$, which gives us the usual SNE
formulation.333In practice, the SNE formulation is slightly different because
it is done in batch mode (using all datapoints $x$ for each update), the
conditional distributions $p(i|w_{j})$ and $q(i|w_{j})$ are set to zero for
$i=j$, and a different $\sigma_{i}$ is used for each stimulus $x=w_{i}$ to
account for different local variations in the data.

It is important to note that Eq. 5 differs from the standard SNE objective
function by the addition of a representation term $-\ln p(x)$, which only
depends on the neuron weights $\\{w_{i}\\}_{i=1}^{N}$ and not on the 2D
points, which is why we recover the same update rule as SNE. If, similarly to
our interpretation of SOM’s algorithm, we interpret SNE as an optimization
over both the neuron weights and the 2D points, this term makes a lot of sense
because the first step in SNE consists in choosing as many neurons as the
number of datapoints and putting a Gaussian $p(x|i)$ around it. This can be
seen as a kernel density estimation method, and hence as a minimization of
$-\ln p(x)$

We have shown that SOMs and SNE are two dual algorithms sharing the same
global objective function over the neuron weights and 2D points, the main
difference being in the way they approach this optimization. SOMs fix the 2D
points (typically, on a grid) and optimize the neuron weights, while SNE fixes
the neuron weights (typically one for datapoint) and optimizes the 2D points
locations. We outline their differences in Table 1.

Table 1: Main differences between SOM and SNE | SOM | SNE
---|---|---
Objective | Reverse KL + Data fitting | Forward KL + Data fitting
| $\mathcal{D}_{KL}(q(i|x),p(i|x))-\ln p(x)$ | $\mathcal{D}_{KL}(p(i|x),q(i|x))-\ln p(x)$
Given | 2D points $\\{z_{i}\\}_{i=1}^{N}$ are given | Neuron weights $\\{w_{i}\\}_{i=1}^{N}$ are given
| (typically on a grid) | (typically one for each datapoint)
Optimized | Neuron weights $\\{w_{i}\\}_{i=1}^{N}$ are learned | 2D points $\\{z_{i}\\}_{i=1}^{N}$ are learned
| (with stochastic gradient descent) | (with gradient descent)

##  2 Quantitative comparisons

We have shown previously that the SOM and SNE algorithms pursue a very similar
objective: making the conditional distribution of the neurons response to a
stimulus $p(i|x)$ and the neighborhood of the winning neuron $q(i|x)$ close to
each other, and fitting the data. Such unified view permits to compare
quantitatively the quality of maps learned with SOM and SNE, and shed light on
their performances. This section is organized as follows: first, we detail in
more details our approach for comparing those techniques. Then we compare them
on a toy two-dimensional dataset. Finally, we compare them on a more realistic
real-world dataset, the Fashion-MNIST dataset.

###  2.1 Comparison approach

The objectives $J^{\text{SOM}}$ and $J^{\text{SNE}}$ we have introduced in
previous section depend on the hyperparameter $\sigma$ governing what it means
to be close in the observation space, and on the scaling of the 2D points
governing what it means to be close in the visualization space. We propose to
take inspiration from [12] and [3] to choose those parameters in a meaningful
way. The entropies in bits of the conditional distributions $p(i|x)$ and
$q(i|x)$ are continuous approximations of the number of neighbors of the
stimulus $x$ (a.k.a. perplexities [6]), respectively in the observation and
the visualization spaces. The method proposed in [12] for comparing the
organization of maps was to look at $k$ neighbors in the observation space,
and look if they are indeed the $k$ neighbors in the visualization space, and
vice versa, for a varying number of neighbors $k$. We propose to take a
similar approach and, for a stimulus $x$ and a given perplexity Perp
(continuous number of neighbors), choose $\sigma$ so that
$2^{H(p(i|x)}=\text{Perp}$, and scale the points $\\{z_{i}\\}_{i=1}^{N}$ so
that $2^{H(q(i|x)}=\text{Perp}$. We will monitor $J^{\text{SOM}}$ and
$J^{\text{SNE}}$ for different perplexities averaged over the datasets , to
assess the quality of the map at different scales. Such approach is therefore
a continuous equivalent of the topographic function introduced in [12] for
assessing the quality of a map organization, but quantifies the data
representation at the same time because the objectives we proposed contain not
only an organization term but also a data fitting term.

In all experiments, we have used an exponential decay for the neighborhood
parameter (inverse of the scaling parameter) and the learning
rate444$\sigma(t)=\sigma_{i}\big{(}\frac{\sigma_{f}}{\sigma_{i}}\big{)}^{t/t_{\text{max}}}$
and
$\lambda(t)=\lambda_{i}\big{(}\frac{\lambda_{f}}{\lambda_{i}}\big{)}^{t/t_{\text{max}}}$.
for the SOM algorithm as suggested in [9]. For SNE, We used the open-source
implementation of t-SNE [8] with a number of degrees of freedom equal to 10000
(hence, considering a Gaussian distribution instead of a t-distribution, as in
SNE) and the other default hyperparameters.

###  2.2 Toy 2D datasets

We compare here quantitatively the SOM and SNE algorithms on a toy two-
dimensional dataset of 1000 points: the moons dataset. A visualization of the
maps learned (in visualization and observation space) is shown in Fig. 2. A
grid of $10\times 10$ was used in these experiments555We chose
$\sigma_{f}=10^{-10}$, $\sigma_{i}=100$, $\lambda_{f}=10^{-10}$,
$\lambda_{i}=1$ and $t_{\text{max}}=100000$, for a fair comparison between SOM
and SNE we did use also 100 neurons for SNE (and not one neuron per
datapoint), those neuron weights have been learned using kMeans [7] with 100
clusters.

We can see in Fig. 3 that SNE clearly outperforms SOM for both objectives and
all perplexities. We believe that this sheds light on the fact that assuming
the points to lie on a grid is a very strong assumption of SOM’s algorithm
that drastically limits its performance capabilities in terms of the
objectives defined above.

SOM

kMeans+SNE

Observation space

SOM

kMeans+SNE

Visualization space

Figure 2: Maps learned on the moons dataset for the SOM and kMeans+SNE
algorithms

Figure 3: Quantitative comparisons on the moons dataset of the maps learned
with SOM and SNE, for the two objectives we proposed (left: SOM objective,
right: SNE objective)

###  2.3 Fashion-MNIST dataset

We propose here to compare quantitatively the SOM and SNE algorithms on a more
realistic dataset: a subset of 10k images taken randomly from the Fashion-
MNIST dataset [13]. We believe that SOM’s assumption of the points lying on a
grid might be less restrictive when using a very large number of neurons,
which would be coherent with [10] arguing that large maps must be used for an
interesting representation to emerge. We learn a $100\times 100$ SOM666 We
chose $\sigma_{f}=10^{-3}$, $\sigma_{i}=10$, $\lambda_{f}=10^{-3}$,
$\lambda_{i}=1$ and $t_{\text{max}}=120000$ on the Fashion-MNIST dataset, and
as we chose 10k images from the original dataset we can use standard SNE with
10k neurons (one for each datapoint), and hence propose a clear and fair
comparison between standard SOM and standard SNE with the same number of
neurons.

We can see qualitatively that our maps have been learned properly for both
methods in Fig. 4, as they have grouped similar images into clusters (note
that the labels used for coloring the points have not been used for training,
but merely for this qualitative visualization). It is interesting to observe
that the maps learned with SOM and SNE share similarities. The yellow and
light green clusters are adjacent and well separated, and are also close to
the two darker green clusters, that are more interrelated. (see left part of
the maps in Fig. 4). Except for the purple cluster that is clearly separated
in both maps, the other clusters on the right hand sides of the maps learned
seem very interrelated in both SOM and SNE maps. Such a qualitative evaluation
shows that both maps are coherent with each other, which also supports our
view that both approaches optimize a similar objective function.

We show the quantitative comparison in Fig. 5, where we can see that the SNE
algorithm outperforms the SOM algorithm for both objectives for all
perplexities. This suggests that, even when using large maps, the assumption
of the points lying on a grid in the visualization space might be too strong
and unrealistic, which limits the quality of the SOM. Note that SNE
outperforms SOM on both objectives considered, including the SOM objective.
This is particularly impressive because SNE is not actually optimizing this
very objective but a different one (with the forward KL divergence instead of
the reverse KL divergence, see Table 1).

SOM

SNE

Figure 4: Maps learned on the Fashion-MNIST dataset (in the 2D visualization
space, colored by the labels)

Figure 5: Quantitative comparisons on the Fashion-MNIST dataset of the maps
learned with SOM and SNE, for the two objectives we proposed (left: SOM
objective, right: SNE objective)

##  3 Conclusion

In this work, we have proposed a unified view on the two widely used
algorithms: Self-Organizing Maps (SOMs) and Stochastic Neighbor Embedding
(SNE). We have shown that both approaches can be seen as the minimization of a
very similar objective function attempting to organize the map and fit the
data, the main difference being that the SOM algorithm fixes the 2D points
(typically on a grid) and learns the neuron weights, while SNE fixes the
neuron weights (typically, one for datapoint) and learns the 2D points. They
can therefore be seen as dual algorithms, and can be compared quantitatively.

We proposed to compare them quantitatively on two datasets: a toy 2D dataset
and the more realistic Fashion-MNIST dataset. In both experiments, we have
observed that SNE significantly outperforms SOM for both objectives
considered. We believe this is due to the SOM assumption that the 2D points
are on a grid, which is very strong assumption that limits the quality of the
map organization.

Most importantly, we hope that our work will inspire methods to combine the
advantages of both approaches. SOMs are bio-inspired, have a fixed number of
neurons, but are not biologically plausible due to the infinite neighborhood
in the beginning. SNE, in its standard form, has as many neurons as the number
of data points, which is not biologically plausible and might pose some
problems in terms of complexity if it were to be applied online on infinitely
growing datasets, but the perplexity is fixed throughout learning and no
neighborhood parameter is required, which alleviates one of the fundamental
drawbacks of SOM’s algorithm: the need to tune and decay the neighborhood
parameter carefully. Also, as we proposed in some of the experiments, our
formulation of the SNE algorithm can be applied on a fixed number of neurons
whose weights have been learned previously (e.g., using the kMeans algorithm).

Methods jointly learning both the positions in 2D and the neuron weights could
definitely be worth considering, and open an interesting avenue for future
work.

## References

  * [1] Bauer, H.U., Pawelzik, K.R.: Quantifying the neighborhood preservation of self-organizing feature maps. IEEE Transactions on neural networks 3(4), 570–579 (1992)
  * [2] De Bodt, E., Cottrell, M., Verleysen, M.: Statistical tools to assess the reliability of self-organizing maps. Neural Networks 15(8-9), 967–978 (2002)
  * [3] Demartines, P., Hérault, J.: Curvilinear component analysis: A self-organizing neural network for nonlinear mapping of data sets. IEEE Transactions on neural networks 8(1), 148–154 (1997)
  * [4] Hinton, G.E., Roweis, S.: Stochastic neighbor embedding. Advances in neural information processing systems 15 (2002)
  * [5] Kohonen, T.: The self-organizing map. Proceedings of the IEEE 78(9), 1464–1480 (1990)
  * [6] Van der Maaten, L., Hinton, G.: Visualizing data using t-sne. Journal of machine learning research 9(11) (2008)
  * [7] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E.: Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830 (2011)
  * [8] Poličar, P.G., Stražar, M., Zupan, B.: opentsne: a modular python library for t-sne dimensionality reduction and embedding. bioRxiv (2019)
  * [9] Rougier, N., Boniface, Y.: Dynamic self-organising map. Neurocomputing 74(11), 1840–1847 (2011)
  * [10] Ultsch, A., Mörchen, F.: Esom-maps: tools for clustering, visualization, and classification with emergent som. Technical Report No. 46, Dept. of Mathematics and Computer Science, University of Marburg, Germany (2005)
  * [11] Van Hulle, M.M.: Self-organizing maps. Handbook of Natural Computing pp. 585–622 (2012)
  * [12] Villmann, T., Der, R., Herrmann, M., Martinetz, T.M.: Topology preservation in self-organizing feature maps: exact definition and measurement. IEEE Transactions on neural networks 8(2), 256–266 (1997)
  * [13] Xiao, H., Rasul, K., Vollgraf, R.: Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms. arXiv (2017)