Ultracold chemistry and its reaction kinetics

A chemical reaction normally occurs at a few hundred kelvin between reagents involving large numbers of particles ($\sim {{10}^{23}}$). This is because reactions are usually activated by thermal fluctuations which are only significant for large concentrations of particles with high momenta. Therefore, the possibility of a chemical reaction taking place in the dilute and ultracold ($T\lt 1\;\mu {\rm K}$) regime is somewhat counterintuitive. However, the formation of ultracold molecules from atoms in a Bose–Einstein condensate (BEC) was observed almost 15 years ago [1]. This interaction of atoms and molecules close to the absolute zero of temperature has been referred to as ultracold chemistry [2].

A variety of experimental techniques such as the coupling of atoms and molecules via magnetic Feshbach resonance [3, 4] have been successfully employed to achieve chemical bonding in an ultracold environment. For example, ultracold Potassium-Rubidium molecules [5] have been investigated to analyse the quantum mechanical effects of particle statistics on molecular reactivity. The now ready experimental accessibility of chemical processes in the dilute ultracold regime strongly motivates us to develop a general physical understanding of their reaction kinetics.

A quantized description is required in order to study the dynamics of ultracold reactions, in order to fully account for the effects of quantum fluctuations and entanglement. Here we should replace the classical notion of a temperature-dependent reaction with a coherent reversible Hamiltonian evolution. The first phenomenological steps toward such a description were taken in [6], where a mean field ansatz was exploited to describe the coherent formation of diatomic molecules in a BEC. Since then a variety of extensions to this model, mainly focussed on adding quantum corrections to the original mean field ansatz, have been studied [7, 8]. Concomitantly, theoretical explanations, via two-body scattering processes, for the ultracold chemical reaction rates observed in recent experiments have been proposed [9]. In the case of fermionic particles with several internal states the observed reaction rates may be understood in terms of a simple quantum threshold model [10]. When there are different particle types, e.g., bosons and fermions, more sophisticated modelling in terms of multichannel quantum defect theory [11] is required. So far, a general and systematic investigation of the dynamics of ultracold reactions, analogous to the study of reaction kinetics for classical thermal reactions, has not yet been undertaken. Such an approach seems to be indispensable if one wants to study the role of quantum coherence and the production of quantum entanglement in ultracold chemical systems.

In this paper we propose a general scheme to study the kinetics of ultracold chemical reactions. Investigating the predicted dynamics, we find that entanglement replaces the role of thermal fluctuations in activating ultracold chemical reactions. This allows us to draw a parallel between ultracold chemical reactions and quenched dynamics. Thanks to the generality of our formulation it is possible to consider complex chemical reactions involving many reagents. We exemplify our approach by studying an experimentally accessible example of an ultracold reaction exhibiting rich quantum phenomena leading to a dynamical relaxation to local equilibrium.

The subject of reaction kinetics in classical chemistry is concerned with the temporal dynamics of a chemical reaction and its reaction rate. One usually studies a general reaction

$$\begin{eqnarray}&&\mathop{\sum }\limits_{i}{{\mu }_{i}}{{A}_{i}}\mathop{\mathop{\rightleftarrows }\limits_{{{k}_{ba}}}}\limits^{{{k}_{ab}}}\mathop{\sum }\limits_{j}{{\nu }_{j}}{{B}_{j}},\end{eqnarray} \tag{ 1 }$$

where the numbers ${{\mu }_{i}}$ and ${{\nu }_{j}}$ are referred to as stoichiometric coefficients and k_ab and k_ba are called reaction constants. If the reaction is elementary, i.e.there are no catalytic, intermediate, or concurrent reactions, the rate equation for the concentration of species A_i, denoted as $[{{A}_{i}}]$, can be inferred straightforwardly [12]:

$$\begin{eqnarray}&&\frac{1}{{{\mu }_{i}}}\frac{d[{{A}_{i}}]}{dt}=-{{k}_{ab}}\mathop{\prod }\limits_{j}{{[{{A}_{j}}]}^{{{\mu }_{j}}}}+{{k}_{ba}}\mathop{\prod }\limits_{k}{{[{{B}_{k}}]}^{{{\nu }_{k}}}}.\end{eqnarray} \tag{ 2 }$$

In the terminology of irreversible reactions, i.e. when either k_ab or k_ba equals zero, the number ${\rm max} \{M,N\}$, where $M:= {{\sum }_{j}}{{\mu }_{j}}$ and $N := {{\sum }_{k}}{{\nu }_{k}}$, is called the order of the reaction. Systems of parallel reactions are treated in an analogous manner. An important feature of classical reaction kinetics is that the temperature dependence of the reaction constants is described by Arrhenius' law:

$$\begin{eqnarray}&&k(T)\propto {\rm exp} \left( \frac{-{{E}_{A}}}{\kappa T} \right),\end{eqnarray} \tag{ 3 }$$

where E_A represents the activation energy of the corresponding reaction. Microscopically this means the statistical fluctuations drive the reaction. This model of classical chemical reaction kinetics has been validated by its long term success in describing a wide range of astonishing phenomena from the oscillating Belousov–Zhabotinskii reaction [13, 14] to deterministic chaos [15]. However, classical reaction kinetics is still an active field, e.g., the question of the existence of an equilibrium was only recently answered in [16].

The description of the dynamics of chemical reactions at ultracold temperatures falls into the framework of quantum mechanics as the basic assumption of Arrhenius' Law, namely a disordered movement of many high-momentum particles, is no longer satisfied in the low temperature regime. Hence we must replace Arrhenius' law with a new concept. Here we propose that quantum entanglement provides the driving source of fluctuations in the ultracold regime. This is because when the system is in a non-equilibrium pure state quantum correlations are generically induced by unitary dynamics⁴ . An analogous situation occurs in the study of quench dynamics for quantum systems [17, 18], i.e. large many-body systems where the interaction undergoes a sudden change at some fixed time. Based on this analogy, we expect the dynamics of ultracold reactions to be dominated by entanglement-induced phenomena, such as the relaxation of subsystems to a local equilibrium.

A chemical reaction is inherently a many particle problem so that second quantization is the natural framework. Hence we describe every reacting species A_i with a corresponding quantum field ${{\hat{\psi }}_{{{A}_{i}}}}$. The underlying Hilbert space is a tensor product of the single-species Fock spaces. The role of the classical particle concentration $[A]$ is replaced by the particle density operator ${{\hat{n}}_{A}}=\hat{\psi }_{A}^{\dagger }{{\hat{\psi }}_{A}}$. The dynamics of the system is induced by its Hamiltonian $\hat{H}={{\hat{H}}_{0}}+{{\hat{H}}_{\operatorname{int}}}$, where ${{\hat{H}}_{0}}={{\hat{H}}_{{\rm kin}}}+{{\hat{H}}_{{\rm pot}}}$ represents the standard second-quantised kinetic and potential terms, respectively. The crucial part of our proposal is now the choice of a proper interaction term ${{H}_{\operatorname{int}}}$ to describe particle conversion, such as molecule formation. Taking guidance from the classical setting we propose the following ${{\hat{H}}_{\operatorname{int}}}$ for the elementary reaction (1):

$$\begin{eqnarray}&&{{\hat{H}}_{\operatorname{int}}}=k\int dx\mathop{\prod }\limits_{i,j}{{\left( \hat{\psi }_{{{A}_{i}}}^{\dagger }(x) \right)}^{{{\mu }_{i}}}}{{\left( {{{\hat{\psi }}}_{{{B}_{j}}}}(x) \right)}^{{{\nu }_{j}}}}+{\rm h}.{\rm c}.,\end{eqnarray} \tag{ 4 }$$

where ${{\hat{\psi }}_{{{A}_{i}}}}$ and ${{\hat{\psi }}_{{{B}_{j}}}}$ obey the canonical commutation relations (CCR) or canonical anti-commutation relations (CAR) according to particle type of the species. If there are multiple concurrent reactions one should take the sum of each of the interaction Hamiltonians ${{\hat{H}}_{\operatorname{int}}}={{\sum }_{i}}{{\hat{H}}_{\operatorname{int},i}}$. Notice that within this framework we are restricted to the description of reversible reactions with a single reaction rate $k\equiv {{k}_{ab}}={{k}_{ba}}$. However, with a little work, the interaction Hamiltonian can be naturally extended via the inclusion of ancillary baths to incorporate more general effects such as a self-interaction of particles, particle loss, and dissipation. For the rest of the paper, we neglect those effects. We make an additional crucial assumption, namely that the positional degrees of freedom are not excited during the reaction. Therefore, we perform a cutoff after the first term in the operator mode expansion, i.e. we replace field operators with annihilation operators according to ${{\hat{\psi }}_{{{A}_{i}}}}\to \phi _{0}^{{{A}_{i}}}(x){{\hat{a}}_{{{A}_{i}}}}$, where $\phi _{0}^{{{A}_{i}}}(x)$ denotes the single-particle ground state wave-function of the respective species. Consequently, the density observable is replaced by the number operator ${{\hat{n}}_{{{A}_{i}}}}=\hat{a}_{{{A}_{i}}}^{\dagger }{{\hat{a}}_{{{A}_{i}}}}$ and ${{\hat{H}}_{0}}$ simplifies to a harmonic oscillator: ${{\hat{H}}_{0}}={{E}_{{{A}_{i}}}}{{\hat{n}}_{{{A}_{i}}}}$. For simplicity we also restrict ourselves to the bosonic systems although the addition of fermionic particles is straightforward. Keeping these approximations in mind we summarise in 1 the proposed interaction Hamiltonians for the first few elementary low-order ultracold chemical reactions.

Table 1.  The proposed interaction Hamiltonians for low-order bosonic reactions.

Order	Reaction	Proposed ${{\hat{H}}_{\operatorname{int}}}$
0.	${\rm bath}\mathop{\rightleftarrows }\limits^{k}A$	$k(\hat{a}_{A}^{\dagger }+{\rm h}.{\rm c}.)$
1.	$A\mathop{\rightleftarrows }\limits^{k}B$	$k(\hat{a}_{A}^{\dagger }{{\hat{a}}_{B}}+{\rm h}.{\rm c}.)$
2.	$A+B\mathop{\rightleftarrows }\limits^{k}C$	$k(\hat{a}_{A}^{\dagger }\hat{a}_{B}^{\dagger }{{\hat{a}}_{C}}+{\rm h}.{\rm c}.)$
3.	$A+B\mathop{\rightleftarrows }\limits^{k}C+D$	$k(\hat{a}_{A}^{\dagger }\hat{a}_{B}^{\dagger }{{\hat{a}}_{C}}{{\hat{a}}_{D}}+{\rm h}.{\rm c}.)$

To complete our proposal, we need to consider suitable initial conditions for our dynamical system. It does not make sense to initialize the system in an equilibrium state, i.e. an eigenstate or thermal state of ${{\hat{H}}_{0}}+{{\hat{H}}_{\operatorname{int}}}$, if one is interested in dynamical effects. Rather, we assume that the system is in an eigenstate of the non-interacting part of the Hamiltionian, which corresponds to the dynamical isolation of the reacting species before the chemical reaction starts. In order to make the comparison to classical reaction kinetics as meaningful as possible, we consider the system to be initialized in a product of coherent states ⁵ , although the generalization to different initial conditions is obvious.

Considering the scheme outlined in 1, we see that a zeroth-order reaction is not a chemical reaction per se, due to the absence of any interaction of different species. It can be rather understood as a connection of the considered species to some reservoir. Neglecting H₀, this coupling ${{\hat{H}}_{\operatorname{int}}}=k(\hat{a}_{A}^{\dagger }+{{\hat{a}}_{A}})$ results in a quadratic scaling of the average particle number of the coupled species $\langle {{\hat{n}}_{A}}(t)\rangle \propto {{t}^{2}}$. Notice, by way of contrast, that in the classical case we obtain a linear growth or decay of the particle concentration $[A]\propto \pm t$ according to the choice of the sign of the coupling constant.

There are two types of elementary first-order reaction. The first describes the simplest of all ultracold chemical reactions, namely a simple particle conversion between two species. The Hamiltonian ${{H}_{\operatorname{int}}}=k(\hat{a}_{A}^{\dagger }{{\hat{a}}_{B}}+{\rm h}.{\rm c}.)$ models a linear interaction between two quantum fields (e.g. a beam-splitter in quantum optics [19]). It is exactly solvable in the sense that it can be linearly transformed into decoupled harmonic oscillators. The average particle number of each species therefore periodically oscillates, in contrast to the classical first-order reaction, which relaxes to a stationary state. The other type of first-order reaction is modelled by ${{H}_{\operatorname{int}}}=k(\hat{a}_{A}^{\dagger }\hat{a}_{B}^{\dagger }+{\rm h}.{\rm c}.)$ and describes the production of 'pairs' from a bath. This interaction is familiar in quantum optics where it models two-mode squeezing.

In order to describe more complex chemical reactions, e.g. the formation of a molecule from two particles, we have to go beyond first-order reactions. The corresponding interaction Hamiltonians describe non-quadratically interacting field theories. These models are not generally exactly solvable and we must employ approximations, heuristics, and numerical methods to study their dynamics.

Here we study the most elementary second-order reaction, diatomic molecule formation [20]:

$$\begin{eqnarray}&&A+A\mathop{\rightleftarrows }\limits^{k}{{A}_{2}}\end{eqnarray} \tag{ 5 }$$

The stoichiometric coefficient for the atomic species ${{\mu }_{A}}$ equals two and the molecular coefficient ${{\nu }_{{{A}_{2}}}}$ is one. According to our proposed framework this reaction is modelled by the Hamiltonian:

$$\begin{eqnarray}&&\hat{H}={{E}_{A}}{{\hat{n}}_{A}}+{{E}_{{{A}_{2}}}}{{\hat{n}}_{{{A}_{2}}}}+k(\hat{a}_{A}^{\dagger }\hat{a}_{A}^{\dagger }{{\hat{a}}_{{{A}_{2}}}}+{\rm h}.{\rm c}.),\end{eqnarray} \tag{ 6 }$$

where E_A and ${{E}_{{{A}_{2}}}}$ label the corresponding ground-state energies. This Hamiltonian has been the subject of much recent research as it also models second harmonic generation [21] or two-photon down conversion [22], as well as molecule formation via coherent photoassociation in an atomic BEC [23–25]. Nonetheless, a comprehensive investigation of its dynamical regimes relevant for the reaction kinetics of ultracold chemistry is missing.

We are interested in studying the role that entanglement plays in the kinetics of the reaction (5). For this purpose, and to avoid having to explore a large parameter space, we consider a regime where the reaction rate is dominant, i.e. $k\gg |{{E}_{A}}|+|{{E}_{{{A}_{2}}}}|$, so that the effect of a change in the binding energy can be neglected.

Additionally, the overall particle number operator ${{\hat{N}}_{tot}}={{\hat{N}}_{A}}+2{{\hat{N}}_{{{A}_{2}}}}$ is a conserved quantity for the dynamics, hence we can restrict our study to the reduced dynamics of the atomic mode, writing $N\equiv {{N}_{A}}$. The behaviour of the molecular mode can be deduced straightforwardly.

In figure 1 we have plotted the full quantum time evolution of the expectation value of $\hat{N}$ (obtained via exact diagonalisation) together with the mean field prediction for the dynamics of the reaction. The system is initially in a product state of a coherent atomic state and a completely depleted molecular state. Mean field theory predicts a complete inversion of the population, where the system is driven to an unstable fixed point [26]. However, considering the full quantum solution, we identify three different dynamical regimes: First, a semi-classical regime, where the quantum and mean field dynamics coincide. At the breakdown time ${{\tau }_{{\rm MF}}}$ the full solution drifts away from the mean field approximation [27] and the semi-classical regime transitions to an intermediate evanescent regime, where the quantum trajectory oscillates with an increasingly damped amplitude. Eventually, the system reaches the asymptotic regime, where the expectation value of the population imbalance relaxes to a stationary value $\bar{N}$.

Zoom In Zoom Out Reset image size

We can understand the three different dynamical regimes by studying the time evolution of the quantum entanglement between the atomic and molecular modes. These results are shown in 2. In the semi-classical regime we see a rapid increase of the entanglement at the beginning of the reaction, which is necessary for the formation of molecules. It is initially rather surprising that the mean field approximation works as well as it does in the semi-classical regime given that the state rapidly becomes entangled and is not well-modelled via a product ansatz. A possible explanation is that the entanglement evolution in the semi-classical regime is typical of that produced by integrable interactions [28], at least until the breakdown time ${{\tau }_{{\rm MF}}}$. After the breakdown time, in the evanescent regime, the system rapidly reaches the maximum available entanglement and begins to explore the full Hilbert space. Soon after, it enters the asymptotic regime where it ergodically evolves through highly entangled states. It remains in the asymptotic regime until it experiences a quantum revival.

Zoom In Zoom Out Reset image size

The dynamical behavior exhibited by the reaction (5) is reminiscent of the local relaxation observed in quenched many particle quantum systems [29]. This hypothesis is supported by studying the time-averaged fluctuations $\Delta {{N}^{2}}={{{\rm lim} }_{\tau \to \infty }}\frac{1}{\tau }\int _{0}^{\tau }dt{{\left( \langle {{{\hat{N}}}_{(t)}}\rangle -\bar{N} \right)}^{2}}$ relative to the mean value $\bar{N}={{{\rm lim} }_{\tau \to \infty }}\frac{1}{\tau }\int _{0}^{\tau }dt\langle {{\hat{N}}_{(t)}}\rangle$, plotted in figure 3: we find that the fluctuations decrease as the particle number N is increased. However, the mechanism leading to the local relaxation observed in quenched dynamics is slightly different to that found here. In quenched many particle systems the incoherent interference of localised excitations traveling at different velocities leads to a cumulative effect of relaxation. However, in our case, we have an interaction between just two modes and the relaxation we observe here is directly related to the growth of entanglement between them, namely the loss of coherence, or purity, of the reduced density operators. Finally, we point out that the relaxation behavior in the asymptotic regime is remarkably similar to the classical high temperature kinetics of (5), which relaxes to the fixed point ${{[A]}^{2}}=[{{A}_{2}}]$, even though our system is always in a global pure state.

Zoom In Zoom Out Reset image size

We obtained $\bar{N}$ and $\Delta {{N}^{2}}$ via exact diagonalisation. Although the full Hamiltonian has degenerate eigenvalues, the dynamical problem, due to the conservation of ${{\hat{N}}_{{\rm tot}}}$, can be separated into finite-dimensional problems with no degeneracy. Exploiting this we find that the time-averaged expectation value of the atoms and the molecules coincides with the predictions given by the diagonal ensemble ${{\bar{N}}_{{\rm ens}}}:={{\sum }_{\alpha }}|{{c}_{\alpha }}{{|}^{2}}{{N}_{\alpha ,\alpha }}$, where ${{c}_{\alpha }}=\langle {{\psi }_{{\rm in}}}|\alpha \rangle$ and ${{N}_{\alpha ,\beta }}=\langle \alpha |\hat{N}|\beta \rangle$ [30, 31] are the coefficients in the complete energy eigenbasis. Moreover, the time-averaged fluctuations around this mean value can be obtained via $\Delta N_{t}^{2}={{\sum }_{\alpha \ne \beta }}|{{c}_{\alpha }}{{|}^{2}}|{{c}_{\beta }}{{|}^{2}}|{{N}_{\alpha ,\beta }}{{|}^{2}}$ [32]. However, the system does not thermalize as the predicted expectation values do not coincide with those of the microcanonical ensemble.

For an experimental implementation of (5) we consider ultracold Cesium atoms:

$$\begin{eqnarray}&&{\rm Cs}+{\rm Cs}\mathop{\rightleftarrows }\limits^{k}{\rm C}{{{\rm s}}_{2}}\end{eqnarray} \tag{ 7 }$$

with ${\rm C}{{{\rm s}}_{2}}={{|}^{3}}\Sigma _{u}^{+},v=0,J=0\rangle$ and ${\rm Cs}=|F=4,{{m}_{F}}=4\rangle$. The chosen atomic state has pure triplet spin character, and we neglect spin mixing effects due to spin-orbit coupling, therefore excluding any vibrational relaxation to the singlet ground state. Potential energy curves of the states under consideration are available through [33–35]. E_A and ${{E}_{{{A}_{2}}}}$ are given by an optical trap confinement and will be different and tunable due to the different dynamic polarizabilities of the molecular and atomic state [36]. For our calculation, ${{E}_{{{A}^{2}}}}\approx 1.05{{E}_{A}}=1.05(\frac{1}{2}m\omega _{A}^{2})$ with ${{\omega }_{A}}=50\;{\rm Hz}.$ The reaction is realized through a 2-color photoassociation in continuous Raman configuration. For the molecule conversion, one has $k=\frac{{{k}_{1}}{{k}_{2}}}{2\Delta }$, where k_i (i = 1,2) are the individual Rabi frequencies and Δ the one-photon detuning. We evaluate the corresponding wave functions and transition dipole matrix elements to and from the intermediate states $|{{b}^{3}}{{\Pi }_{u}},v^{\prime} ,J^{\prime} =0\rangle$ and $|{{c}^{3}}\Sigma _{g}^{+},v^{\prime} ,J^{\prime} =0\rangle$. As an example, for a transition to the $|{{b}^{3}}{{\Pi }_{u}},v^{\prime} ,J^{\prime} =0\rangle$ manifold via lasers operating around 1275 nm and 1232 nm respectively, our calculation yields Rabi frequencies $2\pi \times 1.17{\rm MHz}[0.4{\rm kHz}]\sqrt{I/{\rm mW}\;{\rm c}{{{\rm m}}^{-2}}}$ for the molecular [atomic] coupling respectively. Therefore for a typical Raman setup ($\Delta =750\;{\rm MHz},{{I}_{1}}={{I}_{2}}=(10\;{\rm mW}/0.25\;{\rm m}{{{\rm m}}^{2}})$, hence $k=1.25\;{\rm kHz}$), $k\gt |{{E}_{A}}|+|{{E}_{{{A}_{2}}}}|$ is safely fulfilled.

In classical high-temperature kinetics we need to consider complex reactions involving numerous reactants in order to obtain oscillating or irregular dynamics. However, we will see that when we add to the ultracold diatomic molecule formation a simple zeroth order reaction, we already encounter Hamiltonian chaos in the mean field regime. For this purpose consider the chemical reaction

$$\begin{eqnarray}&&A+A\mathop{\rightleftarrows }\limits^{{{k}_{1}}}{{A}_{2}}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\rm bath}\mathop{\rightleftarrows }\limits^{{{{\rm k}}_{2}}}{\rm A}\end{eqnarray} \tag{ 9 }$$

Applying the proposed rules and assuming that our particles are trapped in the ground state of some harmonic potential, we obtain the following Hamiltonian:

$$\begin{eqnarray}&&\hat{H}={{E}_{A}}{{\hat{n}}_{A}}+{{E}_{{{A}_{2}}}}{{\hat{n}}_{{{A}_{2}}}}+{{k}_{1}}(\hat{a}_{A}^{\dagger }\hat{a}_{A}^{\dagger }{{\hat{a}}_{{{A}_{2}}}}+{{\hat{a}}_{A}}{{\hat{a}}_{A}}\hat{a}_{{{A}_{2}}}^{\dagger })+{{k}_{2}}(\hat{a}_{A}^{\dagger }+{{\hat{a}}_{A}}),\end{eqnarray} \tag{ 10 }$$

where ${{E}_{A}}$ and ${{E}_{{{A}_{2}}}}$ denote the ground state energy of the respective molecular or atomic species. The concurrent zeroth order reaction (9) breaks the conservation of the overall particle number ${{\hat{N}}_{tot}}$ and therefore impedes a full quantum mechanical treatment. Consequently, we investigate the dynamics of the system in the mean field approximation by replacing the creation and annihilation operators in (10) with complex numbers $({{\alpha }_{A}},{{\alpha }_{{{A}_{2}}}})$ labeling coherent states. Note that within this approximation we obtain the average particle number of a certain species $\langle \widehat{{{N}_{i}}}\rangle$ by considering the square of the modulus of the respective complex number $|{{\alpha }_{i}}{{|}^{2}}$.

In the previous section we have seen that the diatomic molecule formation amounts to deviations from mean field dynamics because of the occurrence of quantum effects. Therefore, we need to keep the coupling parameter ${{k}_{1}}$ as small as possible compared to some relevant energies to consider the mean field limit as an appropriate description of the actual dynamics. Keeping this in mind, we obtain the equations of motion from the variational principle:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} i{{{\dot{\alpha }}}_{A}} & = & {{E}_{A}}{{\alpha }_{A}}+2{{k}_{1}}{{{\bar{\alpha }}}_{A}}{{\alpha }_{{{A}_{2}}}}+{{k}_{2}} \\ i{{{\dot{\alpha }}}_{{{A}_{2}}}} & = & {{E}_{{{A}_{2}}}}{{\alpha }_{{{A}_{2}}}}+{{k}_{1}}\alpha _{A}^{2}. \\ \end{array}\end{eqnarray} \tag{ 11 }$$

To reduce the number of parameters in the system, we remove unnecessary degrees of freedom by replacing the dynamical variables with nondimensionalized quantities

$$\begin{eqnarray}&&{{\tilde{\alpha }}_{A}}=\frac{{{\alpha }_{A}}}{{{\alpha }_{0}}},\qquad {{\tilde{\alpha }}_{{{A}_{2}}}}=\frac{{{\alpha }_{{{A}_{2}}}}}{{{\alpha }_{0}}},\qquad \tau =\frac{t}{{{t}_{0}}}.\end{eqnarray} \tag{ 12 }$$

This amounts to the coupled equations

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} i\frac{d{{\tilde{\alpha }}_{A}}}{d\tau }-{{\tilde{\alpha }}_{A}} & = & 1+{{c}_{1}}2{{\overline{\tilde{\alpha }}}_{A}}{{\tilde{\alpha }}_{{{A}_{2}}}} \\ i\frac{d{{\tilde{\alpha }}_{{{A}_{2}}}}}{d\tau }-{{c}_{2}}{{\tilde{\alpha }}_{{{A}_{2}}}} & = & {{c}_{1}}\tilde{\alpha }_{A}^{2}, \\ \end{array}\end{eqnarray} \tag{ 13 }$$

and energy

$$\begin{eqnarray}&&\tilde{H}:=\frac{H}{k_{2}^{2}/{{E}_{A}}}=|{{\tilde{\alpha }}_{A}}{{|}^{2}}+{{c}_{2}}|{{\tilde{\alpha }}_{{{A}_{2}}}}{{|}^{2}}+{{c}_{1}}\left( {{\overline{\tilde{\alpha }}}_{A}}{{\overline{\tilde{\alpha }}}_{A}}{{\tilde{\alpha }}_{{{A}^{2}}}}+{\rm c}.{\rm c}. \right)+({{\overline{\tilde{\alpha }}}_{A}}+{\rm c}.{\rm c}.),\end{eqnarray} \tag{ 14 }$$

with parameters

$$\begin{eqnarray}&&{{t}_{0}}=\frac{1}{{{E}_{A}}},\ {{\alpha }_{0}}=\frac{{{k}_{2}}}{{{E}_{A}}},\ {{c}_{1}}=\frac{{{k}_{1}}{{k}_{2}}}{E_{A}^{2}},\ {{c}_{2}}=\frac{{{E}_{{{A}_{2}}}}}{{{E}_{A}}}.\end{eqnarray} \tag{ 15 }$$

Therefore, the dynamics of the system is completely determined by the choice of the two parameters ${{c}_{1}},{{c}_{2}}\in \mathbb{R}$ and its initial conditions. Note that due to our choice of parameters ${{c}_{1}}=0$ implies ${{k}_{1}}=0$, i.e. no molecule formation, and arbitrary ${{k}_{2}}$. Moreover, the constraint stemming from the validity of the mean field approximation can now precisely be expressed as ${{c}_{1}}\ll 1$ which translates the constraint on the molecular reaction constant to ${{k}_{1}}\ll \frac{E_{A}^{2}}{{{k}_{2}}}$.

What is the physical meaning of those two parameters? The ODE system (13) is equivalent to a pair of nonlinearly coupled harmonic oscillators. Whereas the oscillator describing the atoms has eigenfrequency one, ${{c}_{2}}$ determines the eigenfrequency of the molecular oscillator. Moreover, ${{c}_{1}}$ is the coupling strength between the two oscillators and is the only nonlinear term in the system. We therefore expect a regular behaviour for small values of ${{c}_{1}}$. Moreover, we see from (14) that if ${{c}_{1}}\gg 1$ the system is also integrable due to the conservation of the overall particle number ${{N}_{tot}}=|{{\tilde{\alpha }}_{A}}{{|}^{2}}+2|{{\tilde{\alpha }}_{{{A}_{2}}}}{{|}^{2}}$. In what follows, we investigate the different dynamical regimes determined by the choice of ${{c}_{1}}$ and fix ${{c}_{2}}$ to an experimentally realistic value. Applying perturbation theory for anharmonic oscillators [37] provides us analytic solutions of the occurring dynamical phenomena for sufficiently small coupling constants ${{c}_{1}}$: the trajectories depicted in figure 4 show the expectation value of atoms $|{{\tilde{\alpha }}_{A}}{{|}^{2}}$ and molecules $|{{\tilde{\alpha }}_{{{A}_{2}}}}{{|}^{2}}$ for an initially depleted molecular mode and 40000 atoms. In case of no coupling at all, i.e. ${{c}_{1}}=0$, the molecular mode remains completely depleted, whereas the atomic mode oscillates due to the coupling to the bath. However, as we increase ${{c}_{1}}$ the system progressively enters a modulational regime, in which the molecular site regularly oscillates and the amplitude of the free oscillation on the atomic site is modulated. Let ${{A}_{0}}={{\tilde{\alpha }}_{A}}(0)$ denote the square root of the initial nondimensionalized number of atoms, then we obtain the following analytical expressions for the amplitude ${{A}_{{\rm mod} }}$ and frequency ${{\omega }_{{\rm mod} }}$ of this modulation:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{A}_{{\rm mod} }} & = & \frac{4({{A}_{0}}+1)A_{0}^{3}c_{1}^{2}}{{{({{c}_{2}}-2)}^{2}}}+\mathcal{O}\left( c_{1}^{4} \right), \\ {{\omega }_{{\rm mod} }} & = & {{c}_{2}}-1+\mathcal{O}\left( A_{0}^{2}c_{1}^{2} \right). \\ \end{array}\end{eqnarray} \tag{ 16 }$$

This means increasing ${{c}_{1}}$ causes an quadratic increase of ${{A}_{{\rm mod} }}$ whereas the frequency of modulation ${{\omega }_{{\rm mod} }}$ remains approximately unchanged. The restriction for the parameters for perturbation theory to be valid are $c_{1}^{2}A_{0}^{2}\ll 1\ll {{A}_{0}}$ and ${{c}_{2}}\in (1,2)$.

Zoom In Zoom Out Reset image size

What happens to the system, if we increase ${{c}_{1}}$ beyond the regime of perturbation theory? We already mentioned that ${{c}_{1}}$ interpolates between integrable systems. But does the system remain integrable for all choices of ${{c}_{1}}$? A well-known tool to characterize irregular behaviour of a system is to consider its Poincaré sections [38]. In our case, we choose the quadratures ${{X}_{i}}=\frac{1}{2}({{\alpha }_{i}}+\overline{{{\alpha }_{i}}})$ and ${{P}_{i}}=\frac{1}{2i}({{\alpha }_{i}}-\overline{{{\alpha }_{i}}})$ with $i\in \{A,{{A}_{2}}\}$ as dynamical variables and the surface ${{X}_{{{A}_{2}}}}=0$ as intersection surface. Illustrative examples of these sections for different ${{c}_{1}}$ are shown in figure 5. We find that for a certain range of ${{c}_{1}}$ the system shows behaviour, which is typical for Hamiltonian chaos: as long as the system remains integrable, the Poincaré section consists of closed curves corresponding to sections of two-dimensional tori. However, increasing ${{c}_{1}}$ deforms and finally destroys some of the closed curves. Some of the sampled trajectories start to densely fill out parts of the energy hypersurface. We call this the chaotic regime of the reaction. Finally, further increase of ${{c}_{1}}$ leads to deformation of the energy hypersurface and eventually restores the integrability of the system.

Zoom In Zoom Out Reset image size

We have presented a proposal to systematically investigate the kinetics of ultracold chemical reactions. The dynamical analysis for the most elementary bosonic examples already implies that entanglement plays a major role in the formation of ultracold molecules. This leads to a relaxation of the atomic and molecular subsystems resulting in the stationarity of local observables similar to quenched systems. Moreover, considering two concurrent elementary reactions amounts to complex dynamical behaviour like Hamiltonian chaos in the mean field approximation. Extending our analysis to fermionic or mixed systems and validating the predicted phenomena against experimental data are the next natural steps to be taken.

Helpful discussions with J Eisert, R Werner, E Schreiber, S Schmidt, L C Venuti, B Neukirchen, A Werner and A Milsted are most gratefully acknowledged. This work was supported by the ERC grants QFTCMPS, SIQS, and POLAR and by the cluster of excellence EXC 201 Quantum Engineering and Space-Time Research. We acknowledge support by Deutsche Forschungsgemeinschaft and Open Access Publishing Fund of Leibniz Universität Hannover.