Ultracold Chemistry and its Reaction Kinetics

We study the reaction kinetics of chemical processes occurring in the ultracold regime and systematically investigate their dynamics. Quantum entanglement is found to play a key role in driving an ultracold reaction towards a dynamical equilibrium. In case of multiple concurrent reactions Hamiltonian chaos dominates the phase space dynamics in the mean field approximation.


Introduction
A chemical reaction normally occurs at a few hundred kelvin between reagents involving large numbers of particles (∼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 1 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.

Classical reaction kinetics
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 In the terminology of irreversible reactions, i.e. when either k ab or k ba equals zero, the number 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: 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].

Proposed framework
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 highmomentum 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 4 . 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 ψ 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 ψ ψ = nˆÂ A A † . The dynamics of the system is induced by its k i n p o t 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 int to describe particle conversion, such as molecule formation. Taking guidance from the classical setting we propose the following Ĥ int for the elementary reaction (1): where ψ A i and ψ 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 = ∑ H Ĥî i int int, . Notice that within this framework we are restricted to the description of reversible reactions with a single reaction rate ≡ = k k k ab 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 i denotes the single-particle ground state wave-function of the respective species. Consequently, the density observable is replaced by the number operator = n a âˆÂ A A † i i i and Ĥ 0 simplifies to a harmonic oscillator: 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.
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 + H Ĥ0 i n t , 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 5 , 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 0 , this coupling results in a quadratic scaling of the average particle number of the coupled species 〈 〉 ∝ n t t ( ) A 2 . Notice, by way of contrast, that in the classical case we obtain a linear growth or decay of the particle concentration 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 int † 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 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 nonquadratically interacting field theories. These models are not generally exactly solvable and we must employ approximations, heuristics, and numerical methods to study their dynamics.

Diatomic molecule formation
Here we study the most elementary second-order reaction, diatomic molecule formation [20]: The stoichiometric coefficient for the atomic species μ A equals two and the molecular coefficient ν A 2 is one. According to our proposed framework this reaction is modelled by the Hamiltonian: To avoid experimentally unrealistic situations, we also restrict our considerations to states with finite particle number.
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][24][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 E E A A 2 , so that the effect of a change in the binding energy can be neglected. Additionally, the overall particle number operator = + N N Nˆ2t ot A 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 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 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 semiclassical regime, where the quantum and mean field dynamics coincide. At the breakdown time τ 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 N .
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 τ 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.
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 ( ) , 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 2 , even though our system is always in a global pure state. We obtained N and ΔN 2 via exact diagonalisation. Although the full Hamiltonian has degenerate eigenvalues, the dynamical problem, due to the conservation of N 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 [30,31] are the coefficients in the complete energy eigenbasis. Moreover, the time-averaged fluctuations around this mean value can be obtained via 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: . 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 Figure 2. Evolution of the quantum entanglement (in terms of the von Neumann entropy of the atomic reduced density operator ρ A ) between the atomic and molecular modes for the reaction (5) with respect to rescaled time τ =  tk . For an initially coherent atomic state the amount of entanglement is close to zero for short times. In the vicinity of the breakdown time t MF the entanglement rapidly increases and stays roughly at a constant level for later times. Figure 3. The time-averaged relative mean value of atoms and the fluctuations around it against the overall particle number for the reaction (5). As the particle number is increased the fraction of atoms in the asymptotic regime and the relative temporal fluctuations decrease.
through [33][34][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, with ω = 50 Hz.

Concurrent reaction
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 Applying the proposed rules and assuming that our particles are trapped in the ground state of some harmonic potential, we obtain the following Hamiltonian: 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 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 α α ( , ) labeling coherent states. Note that within this approximation we obtain the average particle number of a certain species 〈 〉 N i by considering the square of the modulus of the respective complex number α | | 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: To reduce the number of parameters in the system, we remove unnecessary degrees of freedom by replacing the dynamical variables with nondimensionalized quantities  with parameters Therefore, the dynamics of the system is completely determined by the choice of the two parameters ∈  c c , 1 2 and its initial conditions. Note that due to our choice of parameters = c 0 1 implies = k 0 1 , 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 1 which translates the constraint on the molecular 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 1 the system is also integrable due to the conservation of the overall particle number 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 α | | , 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) A 0 denote the square root of the initial nondimensionalized number of atoms, then we obtain the following analytical expressions for the amplitude A mod and frequency ω mod of this modulation:  . 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  (16)).
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 ( )  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. Note that in contrast to usual Poincaré sections in the literature, we plotted P A2 on the z-axis to get an better impression of the projected energy hypersurface. The sections are plotted for 25 long-time trajectories with arbitrary initial conditions: for a small perturbation pictured in (a) the system remains integrable. Increasing the perturbation in (b) leads to splitting of the first orbits into little islands according to the Poincaré-Birkhoff theorem [39]. In (c) the irregular trajectories begin to spread out and densely fill out the energy hypersurface. In (d) most of the energy hypersurface is covered by irregular trajectories. Figure (e) shows a significant deformation of the energy hypersurface. In (f) all of the sampled trajectories are again on integrable curves.

Summary and outlook
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.