Complex collisions of ultracold molecules: a toy model

TL;DR

This paper introduces a simplified one-dimensional model for ultracold molecular collisions, revealing insights into the transition from chaotic to non-chaotic potentials and establishing scaling laws for bound state densities relevant to real three-dimensional systems.

Contribution

The study presents a toy model that captures key features of ultracold molecular collisions and derives scaling laws for bound state densities, bridging simplified models and real systems.

Findings

Adiabatic potential distribution transitions from chaotic to non-chaotic with separation.

Four-body bound states show no visible chaos in energy level spacing.

Scaling laws accurately predict bound state density dependence on parameters.

Abstract

We introduce a model to study the collisions of two ultracold diatomic molecules in one dimension interacting via pairwise potentials. We present results for this system, and argue that it offers lessons for real molecular collisions in three dimensions. We analyze the distribution of the adiabatic potentials in the hyperspherical coordinate representation as well as the distribution of the four-body bound states in the adiabatic approximation (i.e. no coupling between adiabatic channels). It is found that while the adiabatic potential distribution transitions from chaotic to non-chaotic as the two molecules are separated, the four-body bound states show no visible chaos in the distribution of nearest-neighbor energy level spacing. We also study the effects of molecular properties, such as interaction strength, interaction range, and atomic mass, on the resonance density and degree of chaos in the adiabatic potentials. We numerically find that the dependence of the four-body bound state density on these parameters is captured by simple scaling laws, in agreement with previous analytic arguments, even though these arguments relied on uncontrolled approximations. This agreement suggests that similar scaling laws may also govern real molecular collisions in three dimensions.

Equations

$$U_{\mathrm{Morse}}(r)=D[1-e^{-a(r-r_0)}{]}^2-D$$

$$U_{\mathrm{P}\ddot{\mathrm{o}}\mathrm{schl}-\mathrm{Teller}}(r)=-\frac{D}{{\cosh (r/r_0)}^2}$$

$${\rho }_1=\sqrt{\frac{{\mu }_{12}}{\mu }}(r_1-r_2)$$

$${\rho }_2=\sqrt{\frac{{\mu }_{34}}{\mu }}(r_3-r_4)$$

$${\rho }_3=\sqrt{\frac{{\mu }_{12,34}}{\mu }}\left(\frac{m_1{r}_1+m_2{r}_2}{m_1+m_2}-\frac{m_3{r}_3+m_4{r}_4}{m_3+m_4}\right),$$

$${\mu }_{12}=\frac{m_1{m}_2}{m_1+m_2}$$

$${\mu }_{34}=\frac{m_3{m}_4}{m_3+m_4}$$

$${\mu }_{12,34}=\frac{(m_1+m_2)(m_3+m_4)}{m_1+m_2+m_3+m_4}.$$

$${\rho }_1$$

$${\rho }_2$$

$${\rho }_3$$

$$R=\sqrt{{\rho }_1^2+{\rho }_2^2+{\rho }_3^2}.$$

$$\Bigl{[}-\frac{1}{2\mu}\frac{\partial^{2}}{\partial R^{2}}+\hat{H}_{\text{ad}}(R,\theta,\phi)\Bigr{]}\psi(R,\theta,\phi)=E\psi(R,\theta,\phi)$$

$${\hat{H}}_{\text{ad}}(R,\theta ,\varphi )=\frac{{\hat{L}}^2(\theta ,\varphi )}{2\mu R^2}+\hat{V}(R,\theta ,\varphi ),$$

$$\psi (R,\theta ,\varphi )=\sum _{\nu }{F}_{\nu }(R){\Phi }_{\nu }(R;\theta ,\varphi ),$$

$${\hat{H}}_{\text{ad}}(R,\theta ,\varphi ){\Phi }_{\nu }(R;\theta ,\varphi )=U_{\nu }(R){\Phi }_{\nu }(R;\theta ,\varphi ).$$

$$\displaystyle\Bigl{[}-\frac{1}{2\mu}\frac{d^{2}}{dR^{2}}+U_{\nu}(R)\Bigr{]}F_{\nu}(R)$$

$$\displaystyle-\frac{1}{2\mu}\sum_{\nu^{\prime}}\Bigl{[}2P_{\nu\nu^{\prime}}(R)\frac{d}{dR}+Q_{\nu\nu^{\prime}}(R)\Bigr{]}F_{\nu^{\prime}}(R)=EF_{\nu}(R)$$

$$P_{\nu {\nu }^{\prime }}(R)=⟨{\Phi }_{\nu }(R,\theta ,\varphi )|\frac{\partial }{\partial R}|{\Phi }_{{\nu }^{\prime }}(R,\theta ,\varphi )⟩$$

$$=\int d\varphi d\theta \sin \theta {\Phi }_{\nu }(R,\theta ,\varphi )\frac{\partial }{\partial R}{\Phi }_{{\nu }^{\prime }}(R,\theta ,\varphi )$$

$$Q_{\nu {\nu }^{\prime }}(R)=⟨{\Phi }_{\nu }(R,\theta ,\varphi )|\frac{{\partial }^2}{\partial R^2}|{\Phi }_{{\nu }^{\prime }}(R,\theta ,\varphi )⟩$$

$$=\int d\varphi d\theta \sin \theta {\Phi }_{\nu }(R,\theta ,\varphi )\frac{{\partial }^2}{\partial R^2}{\Phi }_{{\nu }^{\prime }}(R,\theta ,\varphi ).$$

$$P_{\mathrm{Brody}}(s)={\Gamma\Big{(}\frac{2+q}{1+q}\Big{)}}^{q+1}(1+q)\;s^{q}\exp(-bs^{q+1}).$$

$${\rho }_4\propto N_4{\rho }_a(E).$$

$${\rho }_a(E)\Delta E=$$

$${\int }_{-D}^{E+\Delta E-{ϵ}_1}d{ϵ}_2\rho ({ϵ}_2)\Theta ({ϵ}_1+{ϵ}_2-E),$$

$$\rho_{a}(E)\Delta E\sim\Big{(}\frac{N_{2}}{D}\Big{)}^{2}(D+E)\Delta E$$

$${\rho }_4\propto \frac{N_2^3}{D}.$$

$$N_2=\frac{1}{\pi }[\delta (E=0)-\delta (E=\infty )]$$

$$\delta (E)={\int }_{r_t}^{\infty }dr\sqrt{E-U(r)}-{\int }_0^{\infty }dr\sqrt{E}$$

$$N_2=\frac{1}{\pi }{\int }_{r_t}^{\infty }dr\sqrt{-U(r)}=\sqrt{D}{r}_0,$$

$${\rho }_4\propto r_0^3\sqrt{D}.$$

$${\hat{P}}_{12}{\rho }_1$$

$${\hat{P}}_{12}{\rho }_2$$

$${\hat{P}}_{12}{\rho }_3$$

$${\hat{P}}_{12}R$$

$${\hat{P}}_{12}\theta$$

Full text

I introduction

Ultracold molecules have emerged as a new platform for quantum science and technology Koch et al. (2018); Gadway and Yan (2016). They combine the unique, tunable, coherent setting of ultracold matter with strong dipolar interactions and numerous stable internal rotational and vibrational states. Consequently, this platform has wide-ranging applications, including exploring new phases of matter and nonequilibrium behavior Lepers et al. (2013); Wall et al. (2015); Hazzard et al. (2014); Yan et al. (2013); Carr et al. (2009); Baranov et al. (2012); Bohn et al. (2017); Lemeshko et al. (2013); Blackmore et al. (2018), enabling quantum computation DeMille (2002); André et al. (2006); Yelin et al. (2006); Herrera et al. (2014); Karra et al. (2016); Ni et al. (2018), performing precision measurements, such as measuring the electron electric dipole moment Kozlov and Labzowsky (1995); Flambaum and Kozlov (2007); Hudson et al. (2011); Baron et al. (2014); Cairncross et al. (2017), and studying chemical reactions in the quantum regime Ospelkaus et al. (2010); Ni et al. (2010); de Miranda et al. (2011); Balakrishnan (2016); Krems (2005, 2008); Liu et al. (2018); Guo et al. (2018); Meyer and Bohn (2010).

To fully realize these applications, it is necessary to theoretically understand the collisional behavior of molecules. For example, the density of bi-molecular bound states at the collision energy is important in determining the lifetime of a cloud of molecules Quéméner et al. (2008); de Miranda et al. (2011); Ospelkaus et al. (2010); Morita et al. (2018), the many-body physics in an optical lattice Doçaj et al. (2016); Wall et al. (2017a, b); Ewart et al. (2018), and chemical reaction rates. In addition to being necessary for applications, understanding molecular collisions is of fundamental scientific interest. An important question is whether and how quantum chaos occurs in molecular collisions: molecules are intermediate between atoms, which are mostly simple and non-chaotic, and thermodynamically large systems, which are frequently chaotic. It is believed that quantum systems which become chaotic in the classical limit display spectral fluctuations identical to those of random matrices generated from certain classes of statistical ensembles in random matrix theory Berry et al. (1977); Bohigas et al. (1984); Mehta (2004). A chaotic time-reversal symmetric system is expected to be described by a Gaussian orthogonal ensemble (GOE), with a characteristic level repulsion described by the Wigner-Dyson distribution. On the other hand, levels of non-chaotic, integrable systems follow a Poisson distribution. Manifestations of random matrix theory are well-studied in nuclear spectra and collisions Weidenmüller and Mitchell (2009); Mitchell et al. (2010); Brody et al. (1981).

Molecular collisions are much more complex than atomic collisions, because molecules have a much denser collection of electronic, rotational, vibrational, and hyperfine states, resulting in a several orders of magnitude higher density of resonances. This was suggested by Mayle et al. Mayle et al. (2012, 2013), and recently estimated by Christianen et al. Christianen et al. (2019) to be $0.124\mathrm{\mu K^{-1}}$ for NaK + NaK collision complexes, without including the hyperfine states that will further increase the resonance density. Thus, there will exist many accessible closed-channel bound states around the collision energy – even at the coldest available experimental temperatures. Consequently, while scattering resonances are routinely measured in ultracold atoms Chin et al. (2010), even for atoms with the densest resonance spectrum, such as that explored in lanthanide atoms (Dy and Er) Frisch et al. (2014); Maier et al. (2015a, b); Jachymski and Julienne (2013, 2015); Jachymski (2016); Makrides et al. (2018); Augustovičová and Bohn (2018); Yang et al. (2017) and predicted in alkaline-earth-like atoms (Yb) Green et al. (2016), it is much harder to resolve the resonances in molecular collisions due to the significantly higher density of states. In addition to being difficult to experimentally measure, molecule-molecule collisions are extremely challenging to model, and quantitative calculations remain impossible for diatomic molecules. Recent work has emphasized the high resonance density inherent to bialkali (and heavier) atom-molecule Yang et al. (2018); Croft et al. (2017a, b); Frye et al. (2016); Doçaj et al. (2016); Wall et al. (2017a, b); Mehta et al. (2018) and molecule-molecule Dawid et al. (2018); Gregory et al. (2019); Mehta and Morehead (2015); Mehta (2014) collisions, building on earlier research on lighter molecules formed from lighter atoms Forrey et al. (1998); Avdeenkov and Bohn (2001); Bohn et al. (2002); Tscherbul and Krems (2006); Simoni and Launay (2006); Quéméner et al. (2008); Tscherbul et al. (2009); Simoni et al. (2009); Żuchowski and Hutson (2010); von Stecher et al. (2009).

As a result of the experimental and theoretical challenges, open questions remain, such as: Is chaos universal in molecule collisions? What is the effect of the molecular size, the mass of the constituent atoms, and of the type of interactions?

In this paper, we present a simplified model of the collision of two identical molecules, which can shed light on the questions above while at the same time being amenable to straightforward numerical calculation. The principle simplification is to restrict the atoms to move in one dimension, interacting via a simple model potential (either Morse or Pöschl-Teller). This potential is chosen to be deep enough to harbor many two-atom bound states, mimicking the situation in real molecules. For concreteness, we choose each molecule to be composed of two distinguishable fermionic isotopes.

We find several results that may have implications for scattering of real molecules in three dimensions, which can be divided into three categories. First, we find that the statistics of the adiabatic potentials transition from chaotic to non-chaotic as the intermolecular separation increases, similar to results observed for atom-molecule scattering in calculations using realistic potential energy surfaces Frye et al. (2016). We characterize how the detailed behavior of this crossover depends on the parameters of the model potential. We show that there are clear trends in the average spacing of the adiabatic potentials and the chaoticity of their level spacing statistics as a function of hyperradius, collision energy, molecular mass, and range of the potential. We expect that these qualitative trends persist to three-dimensional systems.

The two remaining findings are more surprising. The second finding results from calculating the four-body bound states in the approximation that the adiabatic potentials are uncoupled. We find that the four-body level statistics show no apparent chaos even when the adiabatic potentials are strongly chaotic. This suggests caution when interpreting the level statistics of adiabatic potentials, which may not lead to corresponding statistics in the spectrum of the four-body complex, particularly if nonadiabatic couplings are weak or negligible. Lastly, and perhaps most interestingly, we derive a simple analytical relation between the four-body density of states and the two-atom interaction parameters. We show that, although the derivation of these analytic relations relies on uncontrolled approximations, they accurately capture the scaling with interaction parameters. This suggests that the analogous expressions for the three-dimensional system will allow one to simply estimate the effects of changing molecular species or other experimental parameters on complex molecule-molecule collisions.

The structure of this paper is as follows. Sec. II describes the problem we solve and the numerical method, which uses the tools associated with hyperspherical coordinates. Sec. III presents the results for the adiabatic potentials. Sec. IV shows the density of four-body bound states in the adiabatic (zero-coupling) approximation from numerically solving the Schrödinger equation, as well as their nearest-neighbor level statistics. It derives approximate analytic expressions for the scaling of the four-body density of states and shows that despite the uncontrolled approximations involved, they agree well with the numerical results. Sec. V summarizes the results, and gives an outlook, especially focusing on these results’ relevance to real three-dimensional molecule-molecule scattering.

II Molecule-molecule scattering and hyperspherical coordinates

We study a system of two identical diatomic molecules in one dimension. Each molecule consists of two fermionic atoms, giving two sets of identical fermions in total, labeled as (1, 2) and (3, 4). We take 1 and 3 to be distinguishable from 2 and 4. The masses of the distinguishable atoms are either set to be equal or different by a small amount. Throughout the paper, we fix the mass of atom 1 (and 2), i.e., $m_{1}=m_{2}=1$, and vary the mass of atom 3 (and 4) to atom 1 by setting $m_{3}/m_{1}=1$ or $1.3$.

We present results for pairwise atom-atom interactions with two potentials: Morse and Pöschl-Teller potentials. In either case, the model potential is assumed to be identical for all pairs of atoms.

[TABLE]

where $r\geq 0$ is the interatomic distance, $D$ is the depth of the potential, $r_{0}$ is the width of the potential, and $a$ is set to $1/r_{0}$, so that $U(0)/D$ is independent of $r_{0}$.

To numerically solve the Schrödinger equation and interpret its eigenstates, we transform it into hyperspherical coordinates. Hyperspherical coordinates have been used widely Delves (1958, 1960); Macek (1968); Bondi et al. (1983); Chuluunbaatar et al. (2008); Lin (1995); Suno et al. (2002); Nielsen et al. (2001); Esry et al. (1996); Mehta et al. (2009); Rittenhouse et al. (2011); Croft et al. (2017a) to study few-body physics and have proven to provide an efficient and powerful method for studying the full scattering and bound state problem. In this representation, the particle coordinates are rewritten into one hyperradial coordinate and a collection of angular coordinates (defined below).

In one dimension, the position and the mass of particle $i$ in the lab frame are labeled by $r_{i}$ and $m_{i}$, respectively. By setting the center-of-mass coordinate to zero, the remaining three degrees of freedom can be expressed in terms of Jacobi coordinates $\rho_{1}$, $\rho_{2}$, and $\rho_{3}$, defined as

[TABLE]

where

[TABLE]

Figure 1 illustrates these Jacobi coordinates for the system studied in this paper. The hyperspherical coordinates $(R,\phi,\theta)$ are related to the Jacobi coordinates by

[TABLE]

where the hyperradius is defined as

[TABLE]

In this formalism, the hyperradius $R$ effectively characterizes the size of the system, and the angular coordinates $\theta$ and $\phi$ describe the relative distances between pairs of atoms. For example, as $(\phi,\theta)\rightarrow(\pi/4,\pi/2)$, the system is separated into two constituent parts, with (1,3) and (2,4) at coincidental locations. When the system is in a state of two well-separated molecules, $R$ measures the intermolecular separation. In this coordinate system, the exchange and parity symmetries of the system can also be readily imposed via boundary conditions for $\theta$ and $\phi$, as discussed in Appendix A. For further information on the hyperspherical representation, refer to Ref. Rittenhouse et al. (2011).

With the hyperspherical coordinates, the Schrödinger equation of the four-body system in 1D can be decomposed into a hyperangular component and a hyperradial component, giving the hyperspherical Schrödinger equation for the reduced wavefunction $\psi(R,\theta,\phi)=R\;\Psi(R,\theta,\phi)$ (where $\Psi(R,\theta,\phi)$ is the wavefunction, and we set $\hbar=1$ throughout this paper),

[TABLE]

where $\hat{H}_{\text{ad}}$ is the adiabatic Hamiltonian,

[TABLE]

$\hat{L}$ is the (hyper)angular momentum operator, and $\hat{V}$ is the interaction between the particles. $\hat{V}$ is obtained by summing up all the pairwise interactions $U(r)$ for each pair of molecules with separation $r$, i.e., $V(R,\theta,\phi)=\sum_{i>j}{U(r_{ij})}$.

It is useful to rewrite $\psi$ as

[TABLE]

where $\Phi_{\nu}(R;\theta,\phi)$ are eigenstates of the adiabatic Hamiltonian at $R$ labeled by $\nu$. The hyperspherical adiabatic potentials are defined as the corresponding eigenvalues $U_{\nu}(R)$ from

[TABLE]

Finally, $F_{\nu}(R)$ are the solutions to

[TABLE]

where

[TABLE]

and

[TABLE]

The above reformulation of the Schrödinger equation is exact, but in this paper we will follow previous literature on complex molecular collisions and concentrate much of our analysis on the adiabatic potential $U_{\nu}(R)$. When we do solve the Schrödinger equation, we will restrict to the fully adiabatic (zero-coupling) approximation, setting $P_{\nu\nu^{\prime}}=Q_{\nu\nu^{\prime}}=0$.

To solve the hyperangular equation Eq. (9), B-splines are used as basis functions in both $\theta$ and $\phi$ coordinates. A detailed description is included in Appendix B. Similarly, the hyperradial equation (10) is solved with B-splines in $R$.

III Numerical results for adiabatic potentials

Figure 2(a) displays the two-atom bound states of Morse potential with $D=50$ and $r_{0}=1$ for mass ratios $m_{3}/m_{1}=1$ and $m_{3}/m_{1}=1.3$. For each case, there are 14 bound states.

Figure 2(c) shows the adiabatic potential curves $U_{\nu}(R)$ obtained by numerically solving the adiabatic equation at each hyperradius $R$ and imposing even parity boundary conditions (see Appendix A for more explanations on the parity symmetry of the system). The left and right plot in Fig. 2(c) correspond to Morse potential with $D=10$ and $D=100$, respectively. The density of adiabatic potential curves increases as the number of bound states for the two-atom interaction increases. Since our interest is in collisions of ground-state molecules, the vertical axis is scaled by $E_{b}$, which is the absolute value of the average lowest dimer-dimer energy at large $R$. $E_{b}$ is the average of $|E_{12}+E_{34}|$ and $|E_{13}+E_{24}|$, the energies of two possible dimer-dimer configurations. The density of adiabatic potentials increases dramatically as the depth of the model potential increases, predicting a high density of four-body bound states. This will be discussed in detail later.

The statistical distribution of adiabatic potentials $U_{v}(R)$ as a function of $R$ is particularly important for understanding quantum chaos Chang (1993); Capecchi et al. (2001); Daily et al. (2015), and one of the most informative and basic measures of the adiabatic potential statistics is shown in Fig. 2(b). Fig. 2(b) shows the probability distribution of nearest-neighbor spacing of the adiabatic potential curves at $R=3.8$, for the $D=100$ case shown in Fig. 2(c). The horizontal axis $s$ is scaled by the average level spacing. At the hyperradius shown, the level distribution is close to a chaotic distribution, with significant level repulsion readily apparent. 144 curves in an energy range $-E_{b}\pm E_{\delta}$ are included in the statistics, where $E_{b}=181.7$ and $E_{\delta}=30$, so that the density of states is approximately uniform over the energy window. The solid curve is the kernel density estimation of the probability distribution given by $p(s)=1/(\sqrt{2\pi}nh)\sum_{i=1}^{n}e^{-[(s-s_{i})/h]^{2}/2}$ where $s_{i}$ are data values and $h=0.2$, an alternative useful to the histogram to estimate $p(s)$. The Wigner-Dyson distribution corresponding to a Gaussian Orthogonal Ensemble and the Poisson distribution, which are predicted to describe chaotic and integrable systems, respectively, are drawn in dashed curves for comparison.

Figure 3 compares the level densities $\rho_{a}$ of the adiabatic potential curves as a function of hyperradius for different two-atom interactions, and demonstrates a universal trend of $\rho_{a}$ first increasing as $R$ increases, attaining a maximum, and then decreasing to a constant value at long range. The Morse potential was used for panels (a)-(e), where it is revealed that $\rho_{a}$ is independent of $D$ and $r_{0}$ for small $R$. This occurs because the kinetic energy dominates in this region, and therefore $\rho_{a}$ is independent of the two-atom interactions. As $R$ increases, $\rho_{a}$ first increases as a result of the decreasing splitting between eigenvalues of the kinetic energy ($\propto 1/R^{2}$), and then decreases and converges to a constant value as the spectrum converges to that of two independent molecules. The peak appears at around $R=6$ for Morse potential with $r_{0}=1$. Fig. 3(a) and Fig. 3(b) shows that $\rho_{a}$ decreases (increases) as the depth(width) of the potential increases. For reference, the number of the two-atom bound states for $D=100$ and $D=200$ are 20 and 28, respectively. Fig. 3(c) demonstrates that $\rho_{a}$ is independent of the width of energy range $E_{\delta}$ over which the density is calculated for a range of widths from $E_{\delta}=5$ to $30$. In Fig. 3(d), no significant shift of $\rho_{a}$ is observed with a change in mass difference. Fig. 3(e) shows that $\rho_{a}$ increases as the center of the energy window used to calculate $\rho$ shifts from low to high ($-1.2E_{b}$ to $-0.8E_{b}$). This suggests the increase of the complexity of the adiabatic potential curves as the system shifts to higher energies. Fig. 3(f) shows the similarity of these features for the Morse and Pöschl-Teller potential.

Figure 4 reveals that the nearest-neighbor spacing distribution of adiabatic potentials evolves from the Wigner-Dyson distribution at $R=3.8$ towards a Poisson distribution as $R$ increases. For each panel in Fig. 4, two curves are highlighted by thick lines to represent cases of small and large $R$, and the other lighter curves show other values of $R$ from $1.6$ to $10$. Each panel shows a different choice of model potential parameters or masses corresponding to the panels in Fig. 3. The chaotic to non-chaotic crossover is clear in panels (a-c) and (e), though with some deviations to be examined later. In panel (d), a similar trend is visible, with level repulsion and a somewhat Wigner-Dyson-like distribution emerge at small $R$, but at large $R$ the distribution deviates dramatically from Poisson, showing a bunching of levels. Panel (f) is the only case where the level repulsion at small $R$ is questionable, showing only a modestly stronger repulsion than for large $R$. We will examine the origins of the features in panels (d) and (f) shortly. (Notice that in these plots, the small dips near $s=0$ that persist to large hyperradii are artifacts of the Gaussian smooth kernel distribution fitting. A more quantitative measure of the chaoticity of the statistical distribution will be presented later.)

A few other features emerge from close examination of Fig. 4(a-f). In Fig. 4(d), where $m_{3}/m_{1}=1$, a spike near $s=0$ arises for large $R$. This is the result of near-degeneracies that occur for large $R$ when $m_{3}/m_{1}=1$. In Fig. 4(f), which shows a Pöschl-Teller potential with $r_{0}=2$, chosen to have the same number of two-atom bound states as panel (a), there is still a visible level repulsion at small $R$ relative to large $R$, but the shape of the curves is always less well described by both the Wigner-Dyson and Poisson distributions. The difference in the short-range behavior for the Pöschl-Teller case stems from the formation of bands of the adiabatic potential curves at short range. Since at small $R$, the value of Pöschl-Teller potential does not vary rapidly, the potential energy surface is nearly constant in the hyperangular coordinates $\theta$ and $\phi$. Hence, the adiabatic potentials separate into bands, where each band is composed of potentials with the same kinetic energy eigenvalue, and the bandwidth – the separation among these potential curves – results from the small differences in their potential energy.

Figure 5 shows a more detailed evaluation of the evolution of chaoticity observed in the level spacing statistics as a function of $R$, by plotting the Brody parameter $q$ Brody (1973); Brody et al. (1981) as a function of the hyperradius 111In Figure 5, because the potential curves are selected from a specific energy range, $q$ has rapid oscillations as $R$ changes, as levels move in and out of the selected energy window. A moving average is applied to smooth the curves to determine their general feature. The average $(R_{i},q_{i})$ is given by $R_{i}=\sum_{n=i}^{n=i+9}f_{n}R_{n}$ and $q_{i}=\sum_{n=i}^{n=i+9}f_{n}q_{n}$, where $R_{n}$ and $q_{n}$ are the hyperradii and Brody parameter evaluated at grid points indexed by $n$, with adjacent $R_{n}$ separated by $\Delta R=0.1$. $f_{n}$ are weights given by a normalized Gaussian with width 1.58.. The Brody distribution is a statistical distribution characterized by $q\in[0,1]$ that can be expressed as

[TABLE]

It reduces to the standard Poisson distribution at $q=0$ and to Wigner-Dyson at $q=1$. The $q$ plotted in Fig 5 is determined by a goodness-of-fit hypothesis test to Eq. (13).

Again, all cases shown in Fig 5 are broadly similar. Take Fig. 5(a) as an example: $q$ smoothly decreases as $R$ increases, at some point falling to zero. The maximum $q$ occurs for small $R$, peaking around $q\sim 0.5$ for $D=100$. Fig. 5(a) also shows no significant change in $q$ by changing the two-atom interaction strength over a range from $D=50$ to $200$. Fig. 5(b) conveys a similar message that $q$ is unaffected by a change in the width of two-atom potential from $r_{0}=0.5$ to $1.5$, after $R$ is rescaled by width $r_{0}$. Fig. 5(c) shows that $q$ is roughly independent of parity, and that when adiabatic curves of both even and odd parities are included, they become doubly degenerate for $R\gtrsim 4$, leading to a zero effective $q$. At small R, the value of q is approximately cut in half when both even and odd parity are included because those sets of spectral lines are uncorrelated with each other, and thus the characteristic level repulsion present in each set separately is not present in the combined set. Fig. 5(d) exhibits an increase in $q$ as the difference in the mass of distinguishable fermions increases. Fig. 5(e) shows that there is a small increase of $q$ as the system is examined at a higher energy. Fig. 5(f) shows that the trends in $q$ for the Pöschl-Teller potential are roughly similar, but $q$ is larger, and with very large oscillations. However the Pöschl-Teller results are poorly behaved and caution must be exercised in interpreting the results.

IV Four-body density of bound states: numerical results and analytic estimates

As mentioned in the introduction, an accurate, easy-to-use estimate of $\rho_{4}$ at the collision energy would be a powerful tool, for example, to determine the lifetime of cold molecule gases. Reference Wall et al. (2017b) derived an analytic formula, involving only simple power law dependences on the potential parameters, relating $\rho_{4}$ to the parameters of the atom-atom interaction for three-dimensional diatomic molecule collisions. A very recent paper by Groenenboom et al. Christianen et al. (2019) derived a similar formula, properly incorporating the angular momentum conservation laws appropriate for zero-field scattering and incorporating more features of the potential energy surface. However, these estimates rely on several uncontrolled approximations whose reliabilities are difficult to gauge.

Here we derive the analogous relations for the one-dimensional problem that we consider, and show by comparison with numerical results that these relations capture the power law scaling of $\rho_{4}$ on two-atom interaction parameters, such as $D$ and $r_{0}$, very accurately. This provides evidence that the analogous analytic formulas for realistic three-dimensional molecular scattering may accurately capture the scaling of $\rho_{4}$ with molecular parameters.

We derive an analytic estimate of $\rho_{4}$ at collision energy $E$ starting with the following relationship:

[TABLE]

where $\rho_{a}(E)$ is the density of adiabatic channels $U_{\nu}(R)$ with threshold energy $E$, and $N_{4}$ is the number of four-body bound states supported by a single channel. In other words, $N_{4}$ is the number of bound states of an effective interaction $U_{\nu}(R)$.

For this approximation to hold, it is important to assume that the adiabatic potentials are identical except for a shift in energy. Although this approximation is substantial, it is roughly true – the width and depth of the adiabatic potentials do not vary rapidly from channel to channel within the energy window of interest. Since at large $R$, the adiabatic potentials must asymptote to pairs of two-atom molecular bound states, their density can be related to the number of two-atom molecular bound states for the two-atom interaction, which we denote as $N_{2}$.

To estimate $\rho_{a}(E)$, we note the number of adiabatic channels with threshold between $E$ and $E+\Delta E$. This is the number of ways that a pair of diatomic molecules can have energy from $E$ to $E+\Delta E$ (assuming $E>-D$),

[TABLE]

where $\rho(E)$ is the two-atom density of states, and $\Theta$ is the Heaviside step function. We make another substantial approximation that $\rho(E)$ is roughly a constant independent of $E$, so $\rho\sim N_{2}/D$. Performing the integration in Eq. (IV) yields

[TABLE]

With the assumption that $U_{\nu}(R)$ has roughly the same shape as the two-atom interaction $U(r)$ (i.e. the shape may be different in details, but has a similar width and depth), we can approximate $N_{4}=c_{2}N_{2}$, where $c_{2}$ is an $O(1)$ constant roughly independent of molecule parameters (still another approximation). Combining with Eq. (14) and Eq. (16), and taking the $E\rightarrow 0$ limit for low energy collisions, we have

[TABLE]

The final step is to obtain $N_{2}$’s dependence on the two-atom potential. By Levinson’s theorem, $N_{2}$ of a given parity for angular momentum $l=0$ is the difference in the phase of scattering $\delta$ at zero and infinity energy

[TABLE]

with a reduced mass $\mu_{2}=1/2$. Using the WKB approximation,

[TABLE]

with $r_{t}$ the classical turning point. With $\delta(E=\infty)=0$,

[TABLE]

in the last step evaluating the integral explicitly for the Morse potential.

Finally, substituting Eq. (20) into Eq. (17), we obtain a scaling relation of $\rho_{4}$ with the depth and width of the two-atom interaction

[TABLE]

Despite these crude approximations, Figure 6 shows that the scaling predicted by Eq. (21) captures the numerically determined $\rho_{4}$ remarkably well. Fig. 6 plots the predicted power laws of $\rho_{4}$ as a function of two-atom interaction depth $D$ and width $r_{0}$ with an overall multiplicative fitting parameter, showing good agreement with the numerical results.

Finally, although we might expect that the chaotic (Wigner-Dyson) distribution of adiabatic channels would lead to a chaotic distribution of four-body bound states, Figure 7 shows that this is not the case – at least within the approximation of zero nonadiabatic couplings. Rather, the nearest-neighbor spacing distribution of the approximated four-body bound states shows a nearly Poissonian distribution of levels, with no apparent signs of chaos. This is surprising given the clear signatures of chaos in the adiabatic potentials over a significant range of $R$. One possible explanation is that the relative spacing of each adiabatic potential curve is large compared to the level spacing of the four-body bound states, so the chaotic signature carried by the hyperangular component is hidden in the bound state distribution. It is likely that including channel couplings $P_{\nu\nu^{\prime}}$ and $Q_{\nu\nu^{\prime}}$ will lead to a chaotic distribution of the four-body eigenstates, which would indicate that chaos in the four-body eigenstates is not directly inherited from the adiabatic potentials, but through a mechanism that relies on the nonadiabatic couplings.

V Conclusions

We have introduced a toy model to understand molecular collisions by considering a 1D system of two identical molecules with identical interactions between each pair of molecules. This allows us to explore a molecular system which has proliferating scattering resonances in a numerically exact fashion. We demonstrate that the hyperspherical adiabatic potentials transition from chaotic (Wigner-Dyson) to non-chaotic (Poisson) with increasing intermolecular separation for certain model potentials, despite the simplicity of the potentials, the reduced dimensionality, and even though the number of two-atom bound states is only $\sim 20$.

We further studied the dependence of the four-atom system’s properties on the two-atom interaction by exploring the parameters space of the model potential. We show that the system at short range becomes more Wigner-Dyson-like as the mass ratio of distinguishable particles $m_{3}/m_{1}$ increases. For both the case of $m_{3}/m_{1}=1$ and the case including both inversion parities, the short-range chaos is suppressed by the extra degeneracy in the system. It is found that changing the depth $D$ or width $r_{0}$ of the two-atom interaction within the order of magnitude does not affect the chaoticity of collisions measured by the Brody parameter. The Brody parameter is found to be greater at higher collision energy.

We also calculated and analyzed the four-body density of states $\rho_{4}$ in the uncoupled approximation. Two remarkable results emerge from this analysis of the four-body bound states. First, we derived a scaling relation $\rho_{4}\propto\sqrt{D}r_{0}^{3}$ under various (uncontrolled) approximations, analogous to the arguments of Ref. Wall et al. (2017b) that relates the density of four-body bound states to the parmeters characterizing the two-atom interaction. Despite the uncontrolled nature of the approximations leading to the simple scaling relations for $\rho_{4}$, the scaling quantitatively agrees with the numerical calculations. Second, the nearest-neighbor level spacing distribution of the four-body bound states shows no signs of chaos despite the clear chaos in the hyperangular coordinates that manifests in the adiabatic potential curves. We anticipate that including the nonadiabatic couplings in future work will lead to level repulsion between states supported by different adiabatic potential curves, revealing a signature of chaos in the full four-body spectrum.

It will be illuminating to explore realistic collisions in three dimensions. Nevertheless, there are several implications of this study for real molecule scattering. Simple (power law) relations between the four-body bound state density and two-atom interaction are found and can be explained by a simple theory. We can speculate based on this evidence that the analogous power laws derived in Ref. Wall et al. (2017b) may accurately predict $\rho_{4}$ for realistic molecule scattering in 3D. If so, this can provide a potent guide for understanding and guiding molecule experiments. Another sage for 3D systems derived from this study is a caution that contrary to expectation, chaos in the distribution in the adiabatic channels does not necessarily translate into the four-body bound state distribution.

Acknowledgements.

This material is based upon work supported with funds from the Welch Foundation, grant no. C-1872. K. R. A. H thanks the Aspen Center for Physics, which is supported by the National Science Foundation grant PHY-1066293, for its hospitality while part of this work was performed.

Appendix A: Symmetry and Boundary Conditions

The exchange symmetry of fermions and spatial inversion symmetry can be used to simplify the problem by reducing the range of angular coordinates of the Schrödinger equation and they also determine the proper boundary conditions. We can define an operator $\hat{P}_{ij}$ that exchanges the particles $i$ and $j$. Acting $\hat{P}_{12}$ on each of the Jacobi coordinates $\rho_{i}$ gives

[TABLE]

The hyperspherical coordinates $(R,\theta,\phi)$ transform as

[TABLE]

Since particles $1$ and $2$ are identical fermions, the wavefunction $\Phi(R,\theta,\phi)$ must be antisymmetric upon exchanging particles 1 and 2. Thus,

[TABLE]

Similarly for $\hat{P}_{34}$,

[TABLE]

From the above two equations, the wavefunction defined on $\phi\in[-\pi,\pi)$ can be determined by knowing its value for $\phi\in[0,\pi/2]$, since based on Eq. (A4), the value of $\Phi(R,\theta,\phi)$ for $\phi\in[-\pi,0]$ can be determined by that for $\phi\in[0,\pi]$, and based on Eq. (A3), $\Phi(R,\theta,\phi)$ for $\phi\in[\pi/2,\pi]$ can be determined by that for $\phi\in[0,\pi/2]$. From Eq. (A3) and (A4), the boundary conditions on this reduced region are derived by letting $\phi\rightarrow 0$ or $\phi\rightarrow\pi/2$, resulting in $\Phi(R,\theta,\phi=0)=0$ and $\Phi(R,\theta,\phi=\pi/2)=0$.

Spatial inversion symmetry $\hat{\Pi}$ also reduces the coordinate region. This symmetry yields $\rho_{i}\to-\rho_{i}$ for $i=1,2,3$, which corresponds to $(\theta,\phi)\to(\pi-\theta,\pi+\phi)$. The wavefunction must be invariant up to a phase under this symmetry, and $\hat{\Pi}^{2}=1$, so

[TABLE]

where $p=-1$ for odd parity, and $p=1$ for even parity. By arguments similar to those above for the fermionic exchange, Eq. (A5) implies that the range of $\theta$ can be reduced from $[0,\pi]$ to $[0,\pi/2]$.

Lastly, an inversion of the system’s geometry, $\rho_{3}\rightarrow-\rho_{3}$, while fixing $\rho_{1}$ and $\rho_{2}$ (equivalent to applying $\hat{\Pi}\hat{P}_{12}\hat{P}_{34}$) should not alter the boundary conditions for the reduced region, i.e.,

[TABLE]

Setting $\theta=\pi/2$, we arrive at another boundary condition: $\Phi(R,\theta,\phi)|_{\theta=\pi/2}=0$ for odd parity and $\frac{\partial}{\partial\theta}\Phi(R,\theta,\phi)|_{\theta=\pi/2}=0$ for even parity.

Appendix B: Solving the Schrödinger Equation

The hyperangular Schrödinger equation

[TABLE]

is solved with B-spline bases Liu et al. (2003) for $R$, $\theta$, and $\phi$, by writing

[TABLE]

$u_{n}(\phi)$ and $v_{m}(\theta)$ are B-spline functions. First, we insert the B-spline expansion in Eq. (B2) into Eq. (Appendix B: Solving the Schrödinger Equation) and integrate both side with $\int{d\Omega}$ over the unit sphere, yielding

[TABLE]

where

[TABLE]

The expression of $t_{\theta}(m^{\prime},m)$ is obtained from integration by parts. Values in Eq. (B4) are numerically integrated.

The hyperradial equation

[TABLE]

is similarly solved on a B-spline basis. For all $R$, $\theta$, and $\phi$, B-splines are of order 5, and all numerical integrals such as those in Eq. (B4) are evaluated with a Gaussian quadrature with 10 nodes. We choose the numbers of B-splines in $\theta$ and $\phi$ directions such that the error of the adiabatic potential curves at the largest plotted hyperradius is $\leq 1\%$ of the average spacing; the number in each direction varies from 100 to 160 depending on the specific parameters of the two-atom interaction. We numerically calculate the eigenvalues and eigenvectors in Eq. (B3) and Eq. (B5) using the Arnoldi algorithm in ARPACK.