A quasiclassical method for calculating the density of states of ultracold collision complexes

TL;DR

This paper develops a quasiclassical method to calculate the density of states for ultracold collision complexes, validating it against quantum results and applying it to predict collision times and loss mechanisms in ultracold molecular systems.

Contribution

It introduces a new quasiclassical approach for DOS calculation in ultracold collisions and applies it to bialkali molecules, providing insights into collision times and loss mechanisms.

Findings

The quasiclassical method agrees well with quantum calculations.

Calculated DOS for NaK-NaK is 0.124 μK^{-1}.

Sticking times are much shorter than previous estimates, suggesting different loss mechanisms.

Abstract

We derive a quasiclassical expression for the density of states (DOS) of an arbitrary, ultracold, $N$-atom collision complex, for a general potential energy surface (PES). We establish the accuracy of our quasiclassical method by comparing to exact quantum results for the K$_2$-Rb and NaK-NaK systems, with isotropic model PESs. Next, we calculate the DOS for an accurate NaK-NaK PES to be 0.124~$\mu$K$^{-1}$, with an associated Rice-Ramsperger-Kassel-Marcus (RRKM) sticking time of 6.0~$\mu$s. We extrapolate the DOS and sticking times to all other polar bialkali-bialkali collision complexes by scaling with atomic masses, equilibrium bond lengths, dissociation energies, and dispersion coefficients. The sticking times calculated here are two to three orders of magnitude shorter than those reported by Mayle et al. [Phys. Rev. A 85, 062712 (2012)]. We estimate dispersion coefficients and collision rates between molecules and complexes. We find that the sticking-amplified three-body loss mechanism is not likely the cause of the losses observed in the experiments.

Figures

Tables

Table 1. Table 1: The estimated DOS, in μ 𝜇 \mu K -1 (RRKM sticking time, in μ 𝜇 \mu s) for all singlet polar bialkali molecules ( J = 0 𝐽 0 J=0 ) in a single hyperfine state. The star indicates the NaK-NaK sticking time has been determined accurately, without extrapolation.

	23Na	39K	87Rb	133Cs
7Li	0.0051 (0.25)	0.014 (0.67)	0.024 (1.17)	0.068 (3.3)
23Na	.	0.124(6.0)*	0.27 (12.9)	0.83 (40)
39K	.	.	0.48 (23.0)	1.50 (72)
87Rb	.	.	.	5.3 (253)

Equations

$$\tau =\frac{2\pi \hslash \rho }{N^{(0)}},$$

$$N^{(\mathrm{cl})}(E,{\mathbf{J}}_0)=$$

$$\times \delta [\mathbf{P}(\mathbf{p})]\delta [\mathbf{X}(\mathbf{x})]\delta [{\mathbf{J}}_0-\mathbf{J}(\mathbf{x},\mathbf{p})],$$

$${\rho }^{(\mathrm{cl})}(E,{\mathbf{J}}_0)$$

$$\times$$

$${\rho }_{JMp}(E)=g_{\mathbf{N}Jp}{\int }_{JM}{\rho }^{(\mathrm{cl})}(E,\mathbf{J})d\mathbf{J},$$

$$1=g_{\mathbf{N},J,1}+g_{\mathbf{N},J,-1}.$$

$${\mathbf{x}}_i=\mathbf{X}+\mathcal{R}(\mathbf{\Omega }){\mathbf{x}}_i^{(\mathrm{bf})}(\mathbf{q}).$$

$$N_{JMp}(E)=g_{\mathbf{N}Jp}{C}_{Nm}{\int }_{JM}d\mathbf{\Omega }d\mathbf{q}d\dot{\mathbf{\Omega }}d\dot{\mathbf{q}}|\det \mathcal{J}(\mathbf{q},\beta ){|}^2\theta \left[E-H(\mathbf{\Omega },\mathbf{q},\dot{\mathbf{\Omega }},d\dot{\mathbf{q}})\right],$$

$$C_{\mathbf{Nm}}=\frac{1}{h^{3N-3}({\sum }_i{N}_i{m}_i{)}^3}\prod _i\frac{m_i^{3N_i}}{N_i!}.$$

$$\mathbf{\omega }=\left(\begin{array}{ccc}0& -\sin \alpha & \cos \alpha \sin \beta \\ 0& \cos \alpha & \sin \alpha \sin \beta \\ 1& 0& \cos \beta \end{array}\right)\dot{\mathbf{\Omega }}.$$

$$N_{JMp}(E)=g_{\mathbf{N}Jp}{C}_{\mathbf{Nm}}{\int }_{JM}d\mathbf{\Omega }d\mathbf{q}d\mathbf{L}d\dot{\mathbf{q}}\frac{|\det \mathcal{J}(\mathbf{q},\beta ){|}^2}{\sin \beta \det \mathcal{I}(\mathbf{q})}\theta \left[E-H(\mathbf{\Omega },\mathbf{q},\mathbf{L},d\dot{\mathbf{q}})\right].$$

$$N_{JMp}(E)=g_{\mathbf{N}Jp}{C}_{\mathbf{Nm}}\int d\mathbf{\Omega }d\mathbf{q}\frac{|\det \mathcal{J}(\mathbf{q},\beta ){|}^2}{\sin \beta \det \mathcal{I}(\mathbf{q})}{\int }_{JM}d\mathbf{L}d\dot{\mathbf{q}}\theta \left[E-V(\mathbf{q})-T_{\mathrm{kin}}(\mathbf{\Omega },\mathbf{q},\mathbf{L},\dot{\mathbf{q}})\right].$$

$${\dot{\mathbf{x}}}_i=\dot{\mathbf{X}}+\mathbf{\omega }\times \mathcal{R}(\mathbf{\Omega }){\mathbf{x}}_i^{(\mathrm{bf})}(\mathbf{q})+\mathcal{R}(\mathbf{\Omega }){\dot{\mathbf{x}}}_i^{(\mathrm{bf})}(\mathbf{q},\dot{\mathbf{q}}).$$

$$T_{\mathrm{kin}}=\sum_{i}\frac{m_{i}}{2}\biggl{\{}\left[\bm{\omega}\times\mathcal{R}(\bm{\Omega})\bm{x}_{i}^{(\mathrm{bf})}(\bm{q})\right]^{2}+\left[\bm{\dot{x}_{i}}^{(\mathrm{bf})}(\bm{q},\bm{\dot{q}})\right]^{2}+\\ 2\left[\bm{\omega}\times\mathcal{R}(\bm{\Omega})\bm{x}_{i}^{(\mathrm{bf})}(\bm{q},\bm{\dot{q}})\right]\cdot\mathcal{R}(\bm{\Omega})\bm{\dot{x}}^{(\mathrm{bf})}_{i}(\bm{q},\bm{\dot{q}})\biggl{\}}.$$

$$\mathbf{L}=\mathcal{R}(\mathbf{\Omega }){\mathcal{I}}^{(\mathrm{bf})}(\mathbf{q})\mathcal{R}(\mathbf{\Omega }{)}^{-1}\mathbf{\omega }=\mathbf{J}-\mathcal{R}(\mathbf{\Omega })\mathbf{j}.$$

$$T_{\mathrm{kin}}=\sum _i\frac{m_i[{\dot{\mathbf{x}}}_i^{(\mathrm{bf})}(\mathbf{q},\dot{\mathbf{q}}){]}^2}{2}-\frac{{\mathbf{j}}^T[{\mathcal{I}}^{(\mathrm{bf})}(\mathbf{q}){]}^{-1}\mathbf{j}}{2}+\frac{{\mathbf{J}}^T\mathcal{R}(\mathbf{\Omega })[{\mathcal{I}}^{(\mathrm{bf})}(\mathbf{q}){]}^{-1}\mathcal{R}(\mathbf{\Omega }{)}^{-1}\mathbf{J}}{2}.$$

$$[{\mathcal{D}}_i{]}_{jk}=m_i\sum _{lm}{ϵ}_{jlm}[{\mathbf{x}}_i^{(\mathrm{bf})}(\mathbf{q}){]}_l[{\mathcal{K}}_i(\mathbf{q}){]}_{mk},$$

$$T_{\mathrm{kin}}(\mathbf{\Omega },\mathbf{q},\mathbf{J},\dot{\mathbf{q}})={\dot{\mathbf{q}}}^T\mathcal{A}(\mathbf{q})\dot{\mathbf{q}}+\frac{{\mathbf{J}}^T\mathcal{R}(\mathbf{\Omega })[{\mathcal{I}}^{(\mathrm{bf})}(\mathbf{q}){]}^{-1}\mathcal{R}(\mathbf{\Omega }{)}^{-1}\mathbf{J}}{2},$$

$$\mathcal{A}(\mathbf{q})=\frac{1}{2}\sum _i{\mathcal{K}}_i(\mathbf{q}{)}^T{m}_i{\mathcal{K}}_i(\mathbf{q})-{\mathcal{D}}_i(\mathbf{q}{)}^T[{\mathcal{I}}^{(\mathrm{bf})}(\mathbf{q}){]}^{-1}{\mathcal{D}}_i(\mathbf{q}).$$

$$N_{JMp}(E)=$$

$$\times \theta \left[E-V(\mathbf{q})-{\dot{\mathbf{q}}}^T\mathcal{A}(\mathbf{q})\dot{\mathbf{q}}-\frac{{\mathbf{J}}^T\mathcal{R}(\mathbf{\Omega })[{\mathcal{I}}^{(\mathrm{bf})}(\mathbf{q}){]}^{-1}\mathcal{R}(\mathbf{\Omega }{)}^{-1}\mathbf{J}}{2}\right].$$

$${\int }_{J}d\mathbf{J}=4\pi {\int }_{B_J}^{B_{J+1}}|\mathbf{J}{|}^{2}d|\mathbf{J}|=\frac{4}{3}\pi (B_{J+1}^3-B_J^3),$$

$$B_{J+1}^3-B_J^3={\left(\frac{2J+1}{2J-1}\right)}^2(B_J^3-B_{J-1}^3).$$

$$B_J=[\frac{1}{3}J(2J-1)(2J+1){]}^{\frac{1}{3}}{B}_1.$$

$$B_J=(\frac{4}{3}{)}^{\frac{1}{3}}J{B}_1.$$

$$N_{JMp}(E)=g_{\mathbf{N}Jp}8{\pi }^3(2J+1){\hslash }^3{C}_{\mathbf{Nm}}\int d\mathbf{q}\frac{|\det {\mathcal{J}}^{\prime }(\mathbf{q}){|}^2}{\det \mathcal{I}(\mathbf{q})}\int d\dot{\mathbf{q}}\theta \left[E-V(\mathbf{q})-{\dot{\mathbf{q}}}^T\mathcal{A}(\mathbf{q})\dot{\mathbf{q}}\right],$$

$$N_{JMp}(E)=\frac{g_{\mathbf{N}Jp}8{\pi }^{3+\frac{D}{2}}(2J+1){\hslash }^3{C}_{\mathbf{Nm}}}{\Gamma (\frac{D}{2}+1)}\int d\mathbf{q}G(\mathbf{q})[E-V(\mathbf{q}){]}^{\frac{D}{2}}.$$

$$G(\mathbf{q})=\frac{|\det {\mathcal{J}}^{\prime }(\mathbf{q}){|}^2}{\det \mathcal{I}(\mathbf{q})\sqrt{\det \mathcal{A}(\mathbf{q})}}$$

$${\rho }_{JMp}(E)=\frac{g_{\mathbf{N}Jp}8{\pi }^{3+\frac{D}{2}}{\hslash }^3{C}_{\mathbf{Nm}}(2J+1)}{\Gamma (\frac{D}{2})}\int d\mathbf{q}G(\mathbf{q})[E-V(\mathbf{q}){]}^{\frac{D}{2}-1}.$$

$$\rho (E)=\frac{16\sqrt{2}{\pi }^3{C}_{\mathbf{Nm}}}{\Gamma (\frac{D}{2}+\frac{3}{2})}\int d\mathbf{q}G(\mathbf{q})\sqrt{\det \mathcal{I}(\mathbf{q})}\{\pi [E-V(\mathbf{q})]{\}}^{\frac{D}{2}+\frac{1}{2}}.$$

$$\rho (E)=\frac{16{\pi }^3{C}_{\mathbf{Nm}}}{\Gamma (\frac{D}{2}+1)}\int d\mathbf{q}G(\mathbf{q})\sqrt{\frac{\det \mathcal{I}(\mathbf{q})}{{\mathcal{I}}_{\mathrm{rot}}(\mathbf{q})}}\{\pi [E-V(\mathbf{q})]{\}}^{\frac{D}{2}},$$

Full text

I Introduction

Ultracold dipolar gases have applications ranging from quantum computation DeMille (2002); Yelin et al. (2006); Ni et al. (2018) and simulation of condensed matter systems Micheli et al. (2006); Büchler et al. (2007); Cooper and Shlyapnikov (2009), to controlled chemistry Krems (2008); Ospelkaus et al. (2010), and high-precision measurements to challenge the standard model Andreev et al. (2018). Ultracold polar, bialkali gases in their absolute ground state have been realized experimentally for nonreactive species such as the bosonic 87Rb133Cs Takekoshi et al. (2014); Molony et al. (2014) and 23Na87Rb Guo et al. (2016) molecules, and the fermionic 23Na40K Park et al. (2015); Seeßelberg et al. (2018). The lifetime of these molecules in the trap is less than a second for the bosonic species Takekoshi et al. (2014); Guo et al. (2016) and a few seconds for 23Na40K Park et al. (2015). The coherence time between hyperfine states of 23Na40K molecules has been shown to approach a second Park et al. (2017). There is potential for improving this further Park et al. (2017) meaning that the trap lifetime of the molecules limits the coherence time. Increasing the lifetime of these molecules is therefore pivotal to realizing applications of these ultracold dipolar gases.

The mechanism limiting the lifetime is currently unknown, but it likely involves ultracold collisions between the molecules Park et al. (2015, 2017); Ye et al. (2018), which have been studied extensively in the literature Quéméner and Julienne (2012); Mayle et al. (2013). In Refs. Ye et al. (2018); Gregory et al. (2019) it is shown that the loss is equally fast as in the case of reactive collisions and that the diatom-diatom collisions are the rate determining step. The current hypothesis is that the loss mechanism involves the formation of long-lived complexes of pairs of diatoms Mayle et al. (2013). These diatoms have long sticking times because of their strong chemical interactions, which gives rise to a high density of states (DOS) and chaotic dynamics.

Croft et al. study ultracold reactive collisions for the triatomic K2+Rb system with converged quantum scattering calculations and report that this required over $300\,000$ hours of CPU time F. E. Croft et al. (2017). For four-atom systems such as NaK+NaK, the computation time will even be orders of magnitude larger, making such calculations unfeasible at this time. Mayle et al. Mayle et al. (2012, 2013) suggested using the Rice-Ramsperger-Kassel-Marcus (RRKM) formalism Levine (2005) to calculate the sticking time $\tau$, i.e., the lifetime of those collision complexes, from the DOS $\rho$,

[TABLE]

where $N^{(0)}$ is the number of states at the transition state. The transition state separates the collision complexes, treated classically, from the pair of colliding molecules, treated quantum mechanically. A surface dividing the two regions can be chosen at some $R=R^{(0)}$, where $R$ is the Jacobi scattering coordinate in the asymptotic region. For ultracold collisions of ground-state nonreactive molecules, there is asymptotically only one open channel, $N^{(0)}=1$. To define the DOS, we may choose $R^{(0)}$ as the smallest intermolecular distance at which $N^{(0)}=1$. However, in practice the DOS already converges for smaller $R$ between 20 and 50 $a_{0}$.

The RRKM theory assumes ergodic dynamics. This assumption was found to be valid for the K+KRb sytem Croft et al. (2017) and should also apply to strongly interacting four-atom sytems that have an even higher DOS. Mayle et al. use simple model potential energy surfaces (PESs), for which the DOS can be calculated quantum mechanically. However, their state counting contained an error, explained in Sec. II.1, that caused an overestimation of the DOS. Furthermore, their method is not applicable to realistic PESs that do depend on the molecular orientation and vibrational coordinates. Nevertheless, their observation that the DOS of ultracold collision complexes is very large and that the RRKM model is a useful tool to calculate the sticking times is very valuable. Furthermore, the DOS is used as a parameter in multichannel quantum defect theory Mayle et al. (2012, 2013), which can be used to describe ultracold scattering. The only needed change in these calculations is to insert the corrected DOS.

Since the DOS is large, a quasiclassical calculation of the DOS is expected to be accurate. In Sec. II.2 we derive a quasiclassical expression for the DOS of an $N$-atom collision complex. Our expression can be applied to PESs that depend on the molecular orientation and vibrational coordinates. In Secs. III.1 and III.2 we validate this method for isotropic vibrational coordinate-independent PESs, such as used by Mayle et al. Mayle et al. (2013), which allows comparing to converged quantum mechanical state counting. In Sec. III.2 we apply our method to our recently calculated NaK-NaK PES Christianen et al. (2019) to accurately compute the DOS of the NaK-NaK system. In Sec. III.3 we extrapolate our results to also estimate the DOS for other polar bialkali collision complexes. Finally, we show in Sec. III.4 that the sticking times are not large enough for a three-body loss mechanism to explain the experimental losses.

II Theory

II.1 Counting angular momentum states

To calculate the DOS quantum mechanically we count all quantum states in a finite energy interval and divide by the size of the interval. Calculating the quantum states is as difficult as solving the scattering problem, so approximations are necessary. In Ref. Mayle et al. (2012) a method was developed to count quantum states for three particle systems described by isotropic, bond-length independent interaction potentials. For such potentials, the DOS calculation is simplified as the angular and vibrational coordinates are uncoupled from one another, as well as from the intermolecular distance. Hence, the DOS can be computed essentially by computing the DOS for the one-dimensional radial problem and subsequently multiplying by the number of contributing ro-vibrational states.

When counting states, it is important to take into account angular momentum conservation, which is done most conveniently in a coupled representation. For atom-diatom systems, coupled states are denoted $|(jl)JM\rangle$, where $j$ is the diatom rotational quantum number, $l$ corresponds to the end-over-end angular momentum, $J$ to the total angular momentum, and $M$ to the projection of the total angular momentum on a space-fixed axis. Thus, for a given $J$ and $M$, there is exactly one quantum state for each pair $(j,l)$ that satisfies the triangular conditions $|j-J|\leq l\leq j+J$. Mayle et al. Mayle et al. (2012), however, counted all uncoupled basis functions $|jm_{j}lm_{l}\rangle$ that have nonzero overlap with specific $J$ (and $M=0$). Hence, for each pair $j,l$ they count $2\min(j,l)+1$ states, rather than one. Because rotational states with $j$ in the low hundreds can contribute energetically, this led to an overestimation of the DOS by two to three orders of magnitude. This has also been noted by Croft et al. F. E. Croft et al. (2017) who corrected this mistake.

Furthermore, also the parity of the collision complex, $p=(-1)^{j+l}$, is conserved. This means that we should only count the states with parity $p=(-1)^{j_{0}+l_{0}}$, where $j_{0}$ and $l_{0}$ are the initial $j$ and $l$ of the collision. This constraint was not taken into account by either Mayle et al. Mayle et al. (2012, 2013) or Croft et al. F. E. Croft et al. (2017).

In the presence of an external field, $J$ is no longer rigorously conserved. In this case, the number of contributing states can be counted by summing over $J$, which can lead to an increase of the DOS by approximately four orders of magnitude. When nonparallel electric and magnetic external fields are present, also cylinder symmetry is broken, and $M$ is no longer conserved, which can increase the DOS by approximately two more orders of magnitude.

The quantum mechanical state-counting method is limited to isotropic and bond length-independent potentials. In the following section, we describe a quasiclassical approach that is also applicable to more general potential energy surfaces (PESs). The quasiclassical approach is accurate precisely because the DOS is so large, meaning we are close to the classical limit.

II.2 Quasiclassical DOS calculation

In this section, we derive an expression for the DOS for an $N$-atom system with a general PES. We compute the number of quantum states from the classical phase-space volume. For a system of $N_{i}$ particles of type $i$, the total number of quantum states below a certain energy, $E$, with given total angular momentum $J_{0}$ and center of mass (C.O.M.) $\bm{X}=(0,0,0)$ is given by

[TABLE]

where $\theta(E)$ is the Heaviside step function, with $\theta(x)=0$ for $x<0$ and $\theta(x)=1$ for $x\geq 0$. The factor $\ \prod_{i}{N_{i}!}$ corrects for indistinguishibility of the particles. The DOS is the derivative of $N^{\mathrm{(cl)}}$ with respect to energy

[TABLE]

The restrictions on the C.O.M. position, ${\bm{X}}({\bm{x}})$, and momentum, ${\bm{P}}(\bm{p})$, ensure their conservation as the C.O.M. motion is uncoupled from the collision dynamics. Finally, the delta function in ${\bm{J}}({\bm{x}},{\bm{p}})$ restricts the classical total angular momentum.

The RRKM sticking time, Eq. (1), scales with the DOS for specific total angular momentum and projection quantum numbers, $J$ and $M$, rather than the sharply defined classical angular momentum ${\bm{J}}$. Therefore, we need to determine integration bounds for the classical total angular momenta that correspond to the specific quantum numbers. The relevant DOS can then be obtained by integrating over this subset of phase space, which we denote symbolically as

[TABLE]

where $p$ denotes the parity. Quantum mechanically the parity of the total wavefunction is conserved during a molecular collision and is therefore a good quantum number. Classically, this parity is not well defined, so a quantum mechanical factor $g_{\bm{N},J,p}$ needs to be introduced, which is defined as the fraction of classical phase space with quantum numbers ${\bm{N}}$ and $J$ that is assigned parity $p$. This factor always obeys the relation:

[TABLE]

In some cases (see, e.g., Sec. III.1), indistinguishability of the atoms and angular momentum conservation restrict the parity, meaning that $g_{\bm{N}Jp}=\delta_{p,1}$ or $g_{\bm{N}Jp}=\delta_{p,-1}$. In most situations, however, half the DOS comes from even parity states and the other half from odd parity states and $g_{\bm{N},J,p}$ approaches $\frac{1}{2}$ for both parities in the limit of large rotational excitations of the collision partners.

For an arbitrary PES, we cannot analytically carry out the integrals over the internal degrees of freedom on which the electronic energy depends. However, assuming the potential depends only on the coordinates, not the momenta, we can carry out all integrals over momenta analytically. This is complicated by the restrictions on momentum and angular momentum conservation. This means that we need to switch to a coordinate system with the total angular momentum, and C.O.M. position and momentum as coordinates. We choose a coordinate system with the minimal number of remaining integrals, which is the number of internal coordinates $D=3N-6$. Next, we need to determine the integration bounds on the classical angular momenta. Finally, we are in a position to integrate over the momenta analytically. The following subsections discuss these three parts of the problem.

II.2.1 Coordinate transformations

We first transform the $3N$ Cartesian coordinates of the $N$ atoms, $\{{\bm{x}}_{i},\,i=1,\ldots,N\}$, to C.O.M. coordinates, $\bm{X}$, $zyz-$Euler angles for the orientation of the complex, $\bm{\Omega}=(\alpha,\beta,\gamma)$, and a set of $3N-6$ internal coordinates, $\bm{q}$, for which we use Jacobi coordinates,

[TABLE]

The body-fixed coordinates $\bm{x}_{i}^{(\mathrm{bf})}$ are transformed to space-fixed coordinates by the $3\times 3$ rotation matrix $\mathcal{R}(\bm{\Omega})$.

The integrals over the delta functions in the C.O.M. position and momentum can now be carried out, which leads to

[TABLE]

where $\mathcal{J}(\bm{q},\beta)$ is the Jacobian matrix for the coordinate transformation of Eq. (6), which is independent of Euler angles $\alpha$ and $\gamma$, and

[TABLE]

We replace the derivatives of the Euler angles by the angular momentum $\bm{L}$ associated with the rotation of the frame of the system. We use $\bm{L}=\mathcal{I}\bm{\omega}$, where $\mathcal{I}$ is the inertial tensor of the system and the angular velocity is

[TABLE]

The Jacobian determinant for this transformation is given by $1/\sin{\beta}$. This gives:

[TABLE]

We assume that the electronic energy depends only on the coordinates, $\bm{q}$, and not their derivatives. This means the above integral can be separated as

[TABLE]

During a molecular collision, the total angular momentum $\bm{J}$ is conserved. To impose this restriction, we replace the integral over $\bm{L}$ by an integral over $\bm{J}=\bm{L}+\mathcal{R}(\bm{\Omega})\bm{j}$. Here, $\bm{j}=\sum_{i}m_{i}\bm{x}_{i}^{(\mathrm{bf})}(\bm{q})\times\bm{\dot{x}_{i}}^{(\mathrm{bf})}(\bm{q},\bm{\dot{q}})$ is the angular momentum in the body-fixed frame, often called the “vibrational angular momentum”, and $\bm{\dot{x}_{i}}^{(\mathrm{bf})}(\bm{q},\bm{\dot{q}})=\mathcal{K}_{i}(\bm{q})\bm{\dot{q}}$, with $[\mathcal{K}_{i}]_{jk}=\partial[\bm{x}_{i}]_{j}/\partial\bm{q}_{k}$. The Jacobian determinant for the transformation from $\bm{L}$ to $\bm{J}$ is unity.

Furthermore, we need an expression for the kinetic energy $T_{\mathrm{kin}}(\bm{\Omega},\bm{q},\bm{J},\bm{\dot{q}})$. The time derivative of coordinates in the space-fixed coordinate system $\bm{\dot{x}}_{i}$ can be written in terms of the coordinates in the body-fixed frame as

[TABLE]

In the C.O.M. frame, the total kinetic energy can now be written as

[TABLE]

We define the inertial tensor, $\mathcal{I}^{(\mathrm{bf})}(\bm{q})$, in the body-fixed frame, and we use

[TABLE]

This yields

[TABLE]

We can write $\bm{\dot{x}}_{i}^{(\mathrm{bf})}(\bm{q},\bm{\dot{q}})=\mathcal{K}_{i}(\bm{q})\bm{\dot{q}}$ and $\bm{j}=\sum_{i}\mathcal{D}_{i}\bm{\dot{q}}$, with

[TABLE]

where $\epsilon_{ijk}$ is the Levi-Civita tensor. Therefore, we can generally write the kinetic energy as a quadratic form in $\bm{\dot{q}}$,

[TABLE]

where $\mathcal{A}(\bm{q})$ is given by

[TABLE]

Substituting this into Eq. (11) gives

[TABLE]

II.2.2 Angular momentum integration bounds

The next step is to determine the integration range for the classical vector $\bm{J}$ that corresponds to a specific quantum number $J$. Quantum mechanically we count the states $|(jL)JM\rangle$. The values $j$ can reach are typically very large for the strongly interacting systems we consider here Mayle et al. (2012), and we are interested in ultracold collisions, meaning $J$ is small. Therefore we use the approximation $j>J$ so that for each value of $j$, there are $2J+1$ allowed values of $L$ for each pair of quantum numbers $J$ and $M$. Since space is isotropic, the DOS does not depend on $M$. Therefore, we integrate over all phase-space regions corresponding to the allowed $M$-values for the given $J$ and subsequently divide by $2J+1$. The total integral over the region corresponding to quantum number $J$ scales as $(2J+1)^{2}$. Classically this is associated with the three dimensional integral over the total angular momentum vector, $\bm{J}$,

[TABLE]

where $B_{J}$ is the lower integration boundary of the classical region that corresponds to the quantum number $J$. It is not directly evident what those boundaries should be. However, we know the integral should be proportional to $(2J+1)^{2}$. We can therefore derive a recurrence relation

[TABLE]

This recursion relation can be solved with $B_{0}=0$ to yield

[TABLE]

If $J\gg 1$ then this expression approaches

[TABLE]

Because the angular momentum at quantum number $J$ is given by $\sqrt{J(J+1)}\hbar\rightarrow(J+\frac{1}{2})\hbar$, the expression for $B_{J}$ should go to $\hbar J$. This means that $B_{1}=\frac{3}{4}^{\frac{1}{3}}\hbar$. This value of $B_{1}$ is also consistent with the quasiclassical quantization, since the integral over $\bm{J}$ and its conjugate variable, $\bm{\Omega}$ should give $h^{3}$ for $J=0$. The integral over $\bm{J}$ with the given value of $B_{1}$ yields $\pi\hbar^{3}$. If this is combined with the integral over $\bm{\Omega}$, which gives $8\pi^{2}$, we obtain $8\pi^{3}\hbar^{3}=h^{3}$, as expected.

II.2.3 Carrying out the integration

Given the integration range corresponding to the total angular momentum $\bm{J}$ we can carry out the integral of Eq. (II.2.1). In the ultracold regime, without an external field breaking angular momentum conservation, $J$ is very small and the energy term $\bm{J}^{T}\mathcal{R}(\bm{\Omega})[\mathcal{I}^{-1}]^{(\mathrm{bf})}\mathcal{R}(\bm{\Omega})^{-1}\bm{J}$ is negligible compared to the interaction energy $V(\bm{q})$. The integral over $\bm{J}$ will therefore yield a constant value of $\pi(2J+1)\hbar^{3}$ and the integrand no longer depends on $\bm{\Omega}$. If the integration over $\bm{J}$ and $\bm{\Omega}$ is carried out, the following expression remains

[TABLE]

where $\mathcal{J^{\prime}}(\bm{q})=\mathcal{J}(\bm{q},\bm{\beta})/\sin\beta$. Note that $\mathcal{J}$ contains one factor $\sin\beta$. The matrix $\mathcal{A}$ is positive definite such that the integral over $\bm{\dot{q}}$ is the volume of a hyperellipsoid. Therefore, with $D$ the dimension of $\bm{q}$, the resulting expression is

[TABLE]

We call the factor

[TABLE]

the geometry factor. The DOS of the system, $\rho=dN/dE$ is given by

[TABLE]

In general, this integral has to be evaluated numerically.

II.2.4 The DOS in presence of external fields

Above we considered the case where $J$, $M$, and $p$ are rigorously conserved, as is the case for any collisional complex in the absence of external fields. However, in the presence of a single external field, the Hamiltonian has cylindrical symmetry, such that $J$ and $p$ are no longer conserved, but $M$ still is. If multiple external fields—say electric and magnetic—occur at an angle to one another, the cylindrical symmetry is also broken, and neither $J$ nor $M$ is rigorously conserved. In the limit of strong fields, all values of $J$ (and $M$) can be populated, whereas in the limit of weak fields, $J$ and $M$ are conserved as discussed above. For intermediate field strengths, the coupling between the different $J$ states is small, meaning that the full parameter space may not be explored within the sticking time. The statistical theory assumes ergodicity but does not quantify the field strength at which the dynamics becomes ergodic in $J$. Purely statistically, we can only treat the strong (or zero) field limit. We assume that even in the strong field limit the interaction of the molecules with the field is small compared to the interaction between the molecules.

When both $J$ and $M$ are not conserved, the phase-space integral is easier than in the case without a field, because we can treat the integration over $\bm{J}$ the same as the integration over $\bm{\dot{q}}$. This leads to a factor $\sqrt{\det[\mathcal{I}(\bm{q})^{-1}/2]}$ in the denominator of Eq. (25) (note that the label (bf) is dropped, since the determinant is invariant under rotation) and an increase of the exponent of the energy by $3/2$, yielding

[TABLE]

In case only $J$ is not conserved, but $M$ still is, the integral is more difficult since then $M$ introduces directionality in space. The derivation for this case is given in the Appendix Sec. VI.1. The result is

[TABLE]

where $\mathcal{I}_{\mathrm{rot}}$ is defined by a series expansion in Appendix Sec. VI.1 and can be interpreted as a weighted average of the eigenvalues of $\mathcal{I}$.

III Results

First, we establish the validity of our quasiclassical approach by considering simple model potentials for K2-Rb and NaK-NaK for which quantum calculations of the DOS are possible. Then we calculate the DOS for a realistic PES and use this result to estimate the DOS for other alkali dimer complexes. We assume for both K2Rb and NaK-NaK that all identical atoms are in the same hyperfine state, meaning they are indistinguishable. If the sticking time is long enough for transitions between hyperfine states to occur during collisions, the DOS increases not just by a factor corresponding to the number of hyperfine states, but because the hyperfine angular momentum couples with the rotational angular momentum, also higher $J$ and $M$ states become accessible, leading to an increase of the DOS by orders of magnitude.

III.1 K2-Rb

For a three-atom system, the geometry factor $G(\bm{q})$ is a simple expression in terms of Jacobi coordinates, $\bm{q}=(R,r,\theta)$. Here, $R$ is the distance beween Rb and the C.O.M. of K2, $r$ is the bond length of the diatom and $\theta$ is the polar angle. For a general three-atom system (AB+C), the expression for the field-free DOS becomes

[TABLE]

Here, $g_{\mathrm{ABC}}=\prod_{i}{N_{i}!}$ is a degeneracy factor to account for indistinguishability, $\mu=(m_{\mathrm{A}}+m_{\mathrm{B}})m_{\mathrm{C}}/(m_{\mathrm{A}}+m_{\mathrm{B}}+m_{\mathrm{C}})$ is the reduced mass of the three-body system and $\mu_{\mathrm{AB}}=m_{\mathrm{A}}m_{\mathrm{B}}/(m_{\mathrm{A}}+m_{\mathrm{B}})$ is the reduced mass of the diatom. This agrees with expressions in the literature for three-body systems, for example Al3 Peslherbe and Hase (1994), except for the degeneracy factor $g_{\mathrm{ABC}}$ and the parity dependent factor $g_{{\bm{N}}Jp}$ which were not taken into account there.

We use K2-Rb as a model system, for which A and B are K and C is Rb. Expressions for the kinetic energy, inertial tensor and the body-fixed angular momentum $\bm{j}$ are given in appendix VI.2. To test the quality of the quasiclassical approximation we use an isotropic, $r$-independent Lennard-Jones interaction potential, as in Refs. Mayle et al. (2012, 2013); F. E. Croft et al. (2017), such that the potential energy is given by

[TABLE]

Here, $C_{12}=C_{6}^{2}/(4D_{e})$ and $C_{6}$ are the Lennard-Jones parameters, and $V_{\mathrm{K}_{2}}$ is the diatom potential of K2. We use $C_{6}=8599~{}E_{h}a_{0}^{6}$ and $D_{e}=1630$ cm*-1*, which are twice the values of the K-Rb potential. This is the same potential as used by Croft et al. F. E. Croft et al. (2017), except that we use for $V_{\mathrm{K}_{2}}$ the diatom potential constructed for our previous work in Ref. Christianen et al. (2019). Just as in Ref. F. E. Croft et al. (2017) we only take into account even $j$ to account for the indistinguishability of the K-atoms. If $j$ is even and $J=0$, then $l=n$ and $p=1$, therefore $g_{\bm{N}Jp}=\delta_{p,1}$.

For such an isotropic $r$-independent PES it is possible to converge the DOS quantum mechanically and to compare this to our quasiclassical results. In the quasiclassical calculation, the remaining integrals in Eqs. (27), (28), and (29), were computed numerically. The numerical integration was done using an integration grid of 56 equidistant points in $r$ ranging from 3.5 to 9 $a_{0}$, 171 equidistant points in $R$ ranging from 3 to 20 $a_{0}$, and 4 points in $\theta$ placed on a Gauss-Legendre quadrature. The large grids in $r$ and $R$ are needed to converge the low energy results. To find the DOS quantum mechanically, we exploit the separation of radial and ro-vibrational degrees of freedom permitted by the isotropic $r$-independent PES. We compute the DOS for the one-dimensional radial problem, subsequently multiply by the number of contributing ro-vibrational states, and finally determine the DOS by binning the quantum states in an interval of 10 cm*-1* and divide their number by the interval length.

In Fig. 1(a) we show the DOS as a function of the energy $E-E_{\mathrm{min}}$, where $E_{\mathrm{min}}$ is the energy of the minimum of the potential. The vertical dashed line indicates the classical dissociation limit for formation of Rb+K2. Quantummechanically, the dissociation energy liest at slightly higher energy because of the zero-point energy of K2. To compute the classical DOS we place the dividing surface at $R=R^{(0)}=20~{}a_{0}$. Above the classical dissociation limit the DOS keeps increasing when we move the dividing surface outwards, but only slowly.

The classical and quantum results agree closely with each other, especially below the classical dissociation limit. In the quantum case there are more fluctuations in the DOS, as expected. In Fig. 1(b) the DOSs are plotted on a double logarithmic scale, both for the case without field and with field(s). Again, the classical and quantum mechanical results agree very well. In the quantum case, the fluctuations become smaller as the DOS becomes larger.

The DOSs in 1(b) show a clear power-law dependence on the energy, $E$. Straight dashed lines with slopes from top to bottom 3, 2.5, and 1.5 are plotted alongside the DOSs to guide the eye. The integrand of Eq. (30) has exponent $1/2$, and when $J$ or $J$ and $M$ are not conserved these exponents become $3/2$ and 2, respectively. Here, we use an isotropic potential, which therefore does not depend on $\theta$. Each “harmonic” degree of freedom contributes $1/2$ to the exponent, such that the exponent would increase by one if the potentials as a function of $R$ and $r$ were perfectly harmonic. The slopes of the graphs are slightly higher near dissociation.

III.2 NaK-NaK

Next, we apply our method to the four-atom NaK-NaK system. We use the Jacobi coordinates $\bm{q}=(R,r_{1},r_{2},\theta_{1},\theta_{2},\phi)$, where $R$ is the NaK-NaK distance, $r_{1}$ and $r_{2}$ are the bond lengths, $\theta_{1}$ and $\theta_{2}$ the polar angles, and $\phi$ is the dihedral angle. In these coordinates, Eq. (27) for the general AB+CD DOS can be written as

[TABLE]

Unlike for the three-atom system, there is no simple analytical expression for $\det{\mathcal{I}(\bm{q})}$ and $\det{\mathcal{A}(\bm{q})}$ for the four-atom system, so we calculated them numerically. For the NaK-NaK system, A and C are K, B and D are Na, and $g_{\bm{N}Jp}=1/2$. The expressions for $\mathcal{I}$ and $\bm{j}$ are given in appendix VI.2. In the quasiclassical calculations for the isotropic PES, we use an equidistant grid in $R$ from 5 to 20 $a_{0}$ with 151 points, a grid of $r_{1}$ from 4.5 to 10 $a_{0}$ with 56 points, a grid of $r_{2}$ ranging from $r_{1}$ to $10.5~{}a_{0}$ with a spacing of 0.1 $a_{0}$. We use a four-point Gauss-Legendre quadrature in $\theta_{1}$ and $\theta_{2}$ and a two-point Gauss-Chebyshev quadrature in $\phi$. We choose $r_{2}>r_{1}$ and multiply the result by a factor two because of the symmetry. An additional factor two is included to compensate for $\phi$ from running up to $\pi$ instead of $2\pi$.

The realistic potential energy surface of NaK-NaK consists of three parts Christianen et al. (2019): two symmetrically equivalent NaK-NaK parts and one Na2-K2 part. Although one set of Jacobi coordinates can in principle describe all arrangements, integrating over these Jacobi coordinates is very difficult, because an increasingly fine angular grid is needed when going further into an arrangement that does not match the chosen coordinates. We therefore construct a separate integration grid in Jacobi coordinates for all three arrangements and add the integrals. In the NaK-NaK arrangement for the realistic potential, we use an equidistant grid in $R$ with 31 points placed from 5 to 20 $a_{0}$. For $r_{1}$ we use a grid of 15 points from 4.5 to $9~{}a_{0}$, and for $r_{2}$ the grid ranges from $r_{1}$ to $10.5~{}a_{0}$, with a spacing of $0.3~{}a_{0}$. A 24-point Gauss-Legendre quadrature between 0 and $\pi$ is used for $\theta_{1}$ and $\theta_{2}$, and an 8-point Gauss-Chebyshev quadrature between 0 and $\pi$ is used for $\phi$. For the Na2-K2 we use a similar grid.

Because there are some overlapping parts of the grids in the center of the PES, we assign a geometry dependent weighting factor to the integrands for each arrangement. This weighting factor $W(\bm{q})$ is based on the symmetrization function in our previous work Christianen et al. (2019). In the NaK-NaK arrangements $W(\bm{q})=W_{1}W_{2}$ or $W(\bm{q})=W_{1}(1-W_{2})$, and in the Na2-K2 arrangement: $W(\bm{q})=1-W_{1}$, with

[TABLE]

We take $W_{1}\equiv W(u_{1},1,1/4)$ with

[TABLE]

where $r_{ij}$ indicates the distance between atom $i$ and $j$. Atoms 1 and 2 are the K-atoms and atoms 3 and 4 are the Na-atoms and $W_{2}\equiv W(u_{2},1/2,1/16)$ with

[TABLE]

Figure 2 shows the DOS for both the isotropic $r$-independent PES (quantum and quasiclassical) and for the realistic PES. For the isotropic PES, the NaK monomer potentials of Ref. Christianen et al. (2019) were used, together with a Lennard-Jones intermolecular potential with parameters $C_{6}=8500$ $E_{h}a_{0}^{6}$. and $D_{e}=4534$ cm*-1*. First, we note that the quasiclassical-quantum correspondence is even better than in the case of K2-Rb. This is not surprising given the DOS is larger by about five orders of magnitude at the dissociation energy. This results in fewer quantum fluctuations in the DOS. At the dissociation energy, the difference in DOS between the isotropic and realistic PESs is about one order of magnitude, both with and without angular momentum conservation. The slope of the DOS in Fig. 2(b) is much larger and less constant for the realistic PES than for the isotropic PES. This is due to anharmonicity and anisotropy of the PES. At the dissociation energy, the DOS for the realistic PES is found to be (for $J=0$): 0.124 $\mu$K*-1* in the field-free case, 2.14 nK*-1* if only $M$ is conserved, and 5.12 pK*-1* when neither $J$ nor $M$ is conserved. These DOS values correspond to RRKM sticking times of 5.96 $\mu$s, 103 ms, and 24.5 s, for these three cases, respectively.

III.3 Extrapolating the NaK-NaK results

In this section we estimate the DOS for other bialkali-bialkali systems by extrapolating the accurate DOS we obtained for NaK-NaK. To find approximate scaling laws, we use the values of $\det\mathcal{I}(\bm{y})$, $\det\mathcal{A}(\bm{y})$ in a planar, antiparallel configuration with the bond lengths $r_{1,2}$ at their equilibrium distances $r_{0}$ and the intermolecular distance $R$ at the minimum distance of the Lennard-Jones potential $R_{0}=(C_{6}/D_{e})^{1/6}$. We assume the potential is isotropic and harmonic, with the force constants $k_{R}=72D_{e}/R_{0}^{2}$ and $k_{r}=\omega^{2}m_{\mathrm{1}}m_{\mathrm{2}}/(m_{\mathrm{1}}+m_{\mathrm{2}})$, with $\omega$ the vibration frequency of the diatom. Furthermore, we consider only one arrangement, meaning we drop one factor $1/2$ for the symmetry. Substituting this into Eq. (32) yields an approximate DOS:

[TABLE]

We use the DOS calculated using our realistic PES to determine the factor $C_{\mathrm{corr}}$ in the above expression, which is meant to correct for the anisotropy and anharmonicity of the PES. We find $C_{\mathrm{corr}}=0.23$. We fix the value of this correction factor, and subsequently evaluate Eq. (36) for all polar bialkali-bialkali systems. We use diatom properties from Ref. Fedorov et al. (2014), and $C_{6}$ coefficients and $D_{e}$ values from Ref. Byrd (2013). The resulting DOSs are listed listed in table 1. We see that—as expected from the equations—the DOS strongly increases when moving from lighter to heavier alkali systems. Here, the reduced mass plays a bigger role than the total mass, e.g., compare NaK to LiCs. We see that the sticking times of the collision complexes, in the absence of chemical reactions, change over three orders of magnitude when moving from 0.25 $\mu$s for NaLi, to 253 $\mu$s for RbCs. Note that for fermionic molecules, $s$-wave scattering is forbidden and that therefore $p$-wave scattering is the dominant mechanism. Therefore $J=1$ and the sticking time is increased by a factor $2J+1=3$, see Eq. (27).

III.4 Sticking-amplified three-body loss

We use the calculated sticking times to study one particular loss mechanism that has been hypothesized to be responsible for the losses observed experimentally: sticking-amplified three-body loss. Here, a free diatom collides with a collision complex, leading to energy transfer from the complex to the diatom and the escape of both the complex and the diatom from the trap. To estimate the rate of this three-body loss process, we need to estimate the rate of complex-molecule collisions and compare the resulting lifetime to the sticking time of the complex.

The rate of complex-molecule collisions can be estimated with a quantum capture model Mayle et al. (2013). The only unknown parameter here is the dispersion ($C_{6}$) coefficient for complex-molecule collisions, which sets the mean scattering length and rate. This dispersion coefficient can be calculated from the dynamic dipole polarizabilities $\alpha(i\omega)$ at imaginary frequencies of both collision partners, A and B, using the Casimir-Polder relation,

[TABLE]

Quantum mechanically, the polarizability for a given state $i$ can be calculated from a sum over states, $f$, where $\omega_{fi}$ is the transition frequency

[TABLE]

From Eq. (38) it is clear that the static dipole polarizability ($\omega=0$) is an upper limit for the polarizability. For a ground state molecule $\omega_{fi}$ is always positive. In the case of diatom-diatom collisions, the dispersion coefficient is mainly due to rotational dispersion Żuchowski et al. (2013) and $\alpha_{\mathrm{diatom}}(0)$ is given approximately by

[TABLE]

where $B$ is the rotational constant and $d$ the dipole moment. The complex is clearly not in the ground state, meaning that terms of the sum in Eq. (38) in energy above and below the energy level of the complex could cancel to some extent, leading to a much smaller polarizability. Quasiclassically, this can be quantified for the static dipole polarizability. We derive in Appendix Sec. VI.3 that this static dipole polarizability can be expressed in terms of the following expectation value

[TABLE]

Note that here $d$ is the total dipole moment of the complex, which depends on the geometry. This expression is remarkably similar to the expression for the free diatom. For the NaK-NaK system, the interaction energy can rise up to $4534$ cm*-1* Christianen et al. (2019). This means that the expectation value of $T_{\mathrm{kin}}$ is in the order of $10^{3}$ cm*-1*, which is four orders of magnitude larger than the rotational constant of NaK, which is $0.095$ cm*-1*. This means that the rotational dispersion contribution to the integral in Eq. (37), will be much smaller than in the diatom-diatom case. Therefore the electronic dispersion term is the most important contribution, which can be estimated to be twice the electronic dispersion coefficient for the diatom-diatom collisions. For NaK, this means that the dispersion coefficient for the complex-diatom collisions, will be $17\,000$ $E_{h}a_{0}^{6}$, which is an order of magnitude smaller than the value of $500\,000$ $E_{h}a_{0}^{6}$ for diatomic collisions, which may be counter-intuitive.

Using the multichannel quantum defect theory from Ref. Mayle et al. (2013) and taking the limit of $T\rightarrow 0$, this dispersion coefficient gives an $s$-wave rate coefficient of $1.1\cdot 10^{-10}$ cm3s*-1*. Multiplying this by a typical density of the diatoms ($4\cdot 10^{10}$ cm*-3* Park et al. (2015)) and taking the inverse gives the lifetime of the complex due to three-body loss. This lifetime is given by $\tau_{3b}=0.23$ s. The sticking time of the NaK-NaK complex for $J=1$ is approximately $18~{}\mu s$, so the complex dissociates much faster than it collides with a third NaK diatom. Therefore, sticking-amplified three-body losses are not the cause of the losses in typical experimentsPark et al. (2015, 2017). Accounting solely for this loss mechanism, the lifetime of the NaK gas in the trap in the experimental conditions would be in the order of hoursPark et al. (2015, 2017). For the RbCs gas of Ref. Takekoshi et al. (2014) it would be tens of minutes. For the NaRb gas such as reported in Ref. Ye et al. (2018), the loss would be on the timescale of a minute, due to the relatively high densities.

The conclusion that three-body collisions are not the cause the experimental losses is based on the sticking time without fields and without taking into account hyperfine transitions. The conclusion may change in the presence of strong electric or magnetic fields, which cause $J$ to no longer be conserved. However, it is not clear from our calculations how strong the external fields need to be to affect the DOS. The DOS, and therefore the sticking time, can also strongly increase in case of hyperfine transitions of the collision complex. However, it is not directly clear whether these occur on the timescale of the sticking time, especially since there are no unpaired electronic spins and the hyperfine transitions must therefore be caused by coupling to the rotational states. The strongest hyperfine coupling is due to the nuclear quadrupole moments interacting with the changing electric field gradients during the collisions. Both the inclusion of hyperfine states into the model or $J$ not being conserved may cause the sticking times to be orders of magnitude larger and may cause the three-body collision mechanism to be more important.

IV Conclusion and discussion

We have derived a quasiclassical equation for the DOS of an ultracold, $N$-atom collision complex, for an arbitrary PES [Eq. (27)]. We have established the accuracy of our quasiclassical method by comparing to exact quantum results for the K2-Rb and NaK-NaK system, with isotropic $r$-independent model PESs. We have calculated the DOS for an accurate NaK-NaK PES to be 0.124 $\mu$K*-1*, with an associated RRKM sticking time of 5.96 $\mu$s. We extrapolate our results to the other bialkali-bialkali systems. The resulting DOS increases rapidly with atomic mass, but only up to 5 $\mu$K*-1* for the heaviest system RbCs, two orders of magnitude below what was reported previously Mayle et al. (2013). Using the resulting sticking times, we conclude a sticking-amplified three-body loss mechanism is not the cause of losses in the experiments.

V Acknowledgements

We thank James Croft and John Bohn for valuable discussions and Dongzheng Yang for carefully checking all derivations in the manuscript. T.K. is supported by NWO Rubicon grant 019.172EN.007 and an NSF grant to ITAMP.

VI Appendix

VI.1 DOS in presence of a field

We are interested in the DOS in the presence of, e.g., an electric field, where $J$ is no longer conserved, but $M$ still is. For the DOS calculation this means that we can no longer treat space isotropically and neglect the $J$-dependence of the kinetic energy. We modify Eq. (25) accordingly, resulting in

[TABLE]

For conserved $M=0$, we integrate over $J_{z}$ from $-1/2$ to $1/2$. This is inaccurate only for the $J=0$ state, but the contribution of $J=0$ to the total DOS is very small if all $J$ are accessible. Because $M_{z}=0$ we can neglect the kinetic energy associated with $J_{z}$, so

[TABLE]

where $M_{zz}$ is the minor of the $z,z$ element of $\mathcal{R}(\bm{\Omega})^{-1}\mathcal{I}(\bm{q})^{-1}\mathcal{R}(\bm{\Omega})$. By Cramer’s rule the minor $M_{zz}$ of a matrix ${\mathcal{X}}$ is equal to $({\mathcal{X}}^{-1})_{zz}\det{\mathcal{X}}$, so

[TABLE]

The denominator in this expression depends only on the first two Euler angles $\alpha$ and $\beta$ and we find

[TABLE]

where we have chosen the lower integration bound, the zeroes of $\bm{\Omega}$, such that the inertial tensor is diagonal and the eigenvalues are ordered in magnitude. The variables $\mathcal{I}_{1}(\bm{q})$, $\mathcal{I}_{2}(\bm{q})$, and $\mathcal{I}_{3}(\bm{q})$ are the eigenvalues of $\mathcal{I}(\bm{q})$, where $\mathcal{I}_{1}(\bm{q})$ is the largest and $\mathcal{I}_{3}(\bm{q})$ the smallest. This choice is possible since we integrate over all angles, such that the integral is independent of the starting point.

The integral over $\beta$ results in

[TABLE]

where

[TABLE]

and $\mathcal{I}_{\mathrm{rot}}$ is the “rotationally averaged” value of $\mathcal{I}_{zz}$. An analytical expression for this integral can be obtained by expanding the logarithms as a power series. Only even powers of $f(\alpha)$ remain and all resulting integrals can be calculated analytically, yielding

[TABLE]

where ${}_{2}F_{1}{}$ is the hypergeometric function. The sum converges rapidly as long as ${\mathcal{I}_{3}}$ is of the same order as $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$. The values this sum can assume lie between $1$ and $\pi/2$. If we substitute this result into Eq. (44), we obtain

[TABLE]

and

[TABLE]

VI.2 Example calculations for K2-Rb and NaK-NaK

For K2-Rb the kinetic energy (for $J=0$) can be written as:

[TABLE]

If we define the $\bm{x}^{\mathrm{(bf)}}$ coordinates by choosing the K2 molecule to be in the $xy$-plane and the $R$ to be along the $x$-axis, then $\bm{j}$ is given by

[TABLE]

and $\mathcal{I}$ is given by

[TABLE]

For the NaK-NaK system, the kinetic energy (for $J=0$) can be written as

[TABLE]

If we choose $r_{1}$ and $R$ to lie in the $xy$-plane, with $R$ along the $x$-axis, then $\bm{j}$ is given by

[TABLE]

For the elements of $\mathcal{I}$ we find

[TABLE]

VI.3 Polarizability of a complex

In Sec. III.4 we argue qualitatively that the static dipole polarizability of the complex is much smaller than for the diatoms. Here we use our quasiclassical formalism to express the polarizability as an expectation value over phase space. For the complex, we can calculate the expectation value of the static dipole polarizability in our quasiclassical framework. If we re-introduce the integral over the Euler angles $\bm{\Omega}$ in Eq. (27) and introduce an external field $\bm{F}$, we can write the DOS as

[TABLE]

where $\delta\rho(\bm{q},\bm{\Omega})$ is given by

[TABLE]

Here $f(\bm{q},\bm{\Omega},\bm{F})$ is a perturbation on the energy caused by the external field $\bm{F}$, which we assume to be small enough for $\mathcal{J}$ to still be (approximately) conserved. Then the expectation of $\chi(\bm{q},\bm{\Omega})$ can be calculated as

[TABLE]

The polarizability tensor $\alpha$ is given by

[TABLE]

where $\bm{\mathcal{E}}$ is the electric field and $\bm{d}(\bm{q},\bm{\Omega})$ the electric dipole moment. The electric dipole moment is given by the vector sum of the dipoles of the two NaK molecules. These molecular dipoles lie along the molecular axes, the directions of which depend on $\bm{q}$ and $\bm{\Omega}$.

The interaction energy of the system with an electric field is given by $-\bm{d}\cdot\bm{\mathcal{E}}$. The expectation value of the electric dipole moment, for a given $J$, is given by

[TABLE]

Because of the integration over $\bm{\Omega}$, the expectation value of the dipole moment in the weak field limit vanishes. For the polarizability, the off-diagonal components integrate to zero, but the diagonal components do not. The expectation values of the diagonal polarizability components (if we take $\bm{\mathcal{E}}=0$) are given by

[TABLE]

If we introduce the isotropic polarizability as $\alpha_{0}=\frac{1}{3}(\alpha_{xx}+\alpha_{yy}+\alpha_{zz})$, we obtain

[TABLE]

Comparing to Eq. (57), this expression can be written as

[TABLE]

For a diatom-diatom complex $D=6$ and therefore the static dipole polarizability becomes

[TABLE]