Chemical Reactions Using a Non-Equilibrium Wigner Function Approach

Abstract

A three-dimensional model of binary chemical reactions is studied. We consider an ab initio quantum two-particle system subjected to an attractive interaction potential and to a heat bath at thermal equilibrium at absolute temperature $T>0$. Under the sole action of the attraction potential, the two particles can either be bound or unbound to each other. While at $T=0$, there is no transition between both states, such a transition is possible when $T>0$ (due to the heat bath) and plays a key role as $k_{\mathrm{B}}T$ approaches the magnitude of the attractive potential. We focus on a quantum regime, typical of chemical reactions, such that: (a) the thermal wavelength is shorter than the range of the attractive potential (lower limit on T) and (b) $(3/2)k_{\mathrm{B}}T$ does not exceed the magnitude of the attractive potential (upper limit on T). In this regime, we extend several methods previously applied to analyze the time duration of DNA thermal denaturation. The two-particle system is then described by a non-equilibrium Wigner function. Under Assumptions (a) and (b), and for sufficiently long times, defined by a characteristic time scale D that is subsequently estimated, the general dissipationless non-equilibrium equation for the Wigner function is approximated by a Smoluchowski-like equation displaying dissipation and quantum effects. A comparison with the standard chemical kinetic equations is made. The time τ required for the two particles to transition from the bound state to unbound configurations is studied by means of the mean first passage time formalism. An approximate formula for τ, in terms of D and exhibiting the Arrhenius exponential factor, is obtained. Recombination processes are also briefly studied within our framework and compared with previous well-known methods.

1. Introduction

While the foundations of equilibrium statistical mechanics are well established [1,2,3,4,5], those of non-equilibrium statistical mechanics continue to be an open and active important research area [1,6,7,8]. The monographs by [9,10,11] provide a wide and accessible perspective on classical equilibrium and irreversible thermodynamics. Within the realm of quantum thermodynamics (see [12] and the references therein), where the laws of thermodynamics from quantum mechanics and the theory of open quantum systems remain a mainstay [13], particular attention has been put to the Markovian assumptions and the dissipative dynamics proposed by Lindblad and Gorini–Kossakowski–Sudarshan [14,15]. One crucial hallmark is that the off-equilibrium evolution of statistical systems displays stochasticity: see [16,17,18] in the classical regime and [13,19,20,21] in the quantum one.

It is the goal of the present study to focus on non-equilibrium statistical mechanics in the quantum regime, which faces many open fundamental problems. A key question is how to formulate suitably approximate dynamical equations, displaying stochasticity, by starting out from a quantum mechanical framework and carrying out adequate long time approximations. Such a task turns out to be very difficult for isolated and large nonequilibrium quantum systems.

Fortunately, certain rewarding results have been obtained for a closed composite system consisting of the system $S_{\mathrm{int}}$ of actual interest and a larger “heat bath” (HB) in quantum mechanics, under various conditions. Those results required considering the quantum states of $S_{\mathrm{int}}$ and of the HB altogether and various conditions on Hamiltonians and that the states of $S_{\mathrm{int}}$ be spread over many different energies. In short, even if the composite is initially in a pure nonequilibrium state, the reduced state of $S_{\mathrm{int}}$ will tend to be canonical. Some of those results (which are related to and support our present work) are: (i) derivation of the canonical density operator for $S_{\mathrm{int}}$ for large times [22] and for a suitable overwhelming majority of wave functions (typicality) [23]; (ii) $S_{\mathrm{int}}$ will approach an equilibrium state and remain close to it for almost all times, independently both of the HB and, under suitable conditions, also of the initial state of $S_{\mathrm{int}}$ [24]. For equilibration of isolated systems, see [25]. For further strengthening results, which include a number of extensions, small $S_{\mathrm{int}}$ and equilibration for finite times, under wider conditions on energy spreads of states, as well as the interplay between quantum information and thermodynamics, see [26,27,28,29] and the references therein.

Chemical reactions constitute a very important and wide class of phenomena in which quantum mechanics and (equilibrium and non-equilibrium) statistical mechanics play essential roles. See [4,5] and the references therein for studies devoted to the equilibrium chemical constant. Much research has also been devoted to non-equilibrium chemical reactions: see, in particular, [9,11,17] for stochastic and thermodynamical approaches. The non-equilibrium rate constant is treated in [17] by means of a thorough analysis of the Kramers equation. Pulsed lasers allow one to probe, with very short temporal resolution, off-equilibrium dynamics in chemical reactions (see [30,31] and the references therein).

In biological systems, one prominent example is the thermal denaturation of DNA: a physical non-equilibrium process through which the two homologous strands, initially bound to each other, split off into two separate single ones, at a certain (melting) temperature [32,33]. In a previous work [34], and motivated by non-equilibrium DNA thermal denaturation, the present authors treated the dynamics of two interacting classical three-dimensional macromolecular chains. Each chain was assumed to be Gaussian (except for possible additional weak interactions along it). These chains were immersed in a medium at thermal equilibrium and initially bound to each other by potentials, which were attractive for medium and large distances (say, Morse-like). As temperature increased from room values up to about the melting one, the two strands became dissociated, and the mean first passage time of this process was calculated. Here, we shall generalize nontrivially the methods used in [34] to a different process, namely a non-equilibrium chemical reaction of two quantum particles (atoms and small molecules) in a medium, by resorting to a somewhat simplified model.

The strategy in our present study will parallel that of [34] very closely. Then, for a better understanding of the present work, it is adequate to start by summarizing below the main successive conceptual steps followed in [34]:

• The two macromolecular chains were described, from the outset, by a classical (non-negative) probability distribution, depending on all spatial positions of all atoms and on time and evolving through an irreversible Kramers-like master equation (stochasticity being thereby taken care of). Neither Planck’s constant ℏ nor thermal wavelengths were included. However, the formulation did include an important length scale, namely the bond length yielding the average distance between two successive atoms along each chain.
• The overall center of mass (CM) motion was factored out off-equilibrium from the remaining spatial variables (the set $\left[\mathit{z}\right]$ of all positions of the atoms along the two chains and the relative position y of the two centers of mass of the chains). The interest focused then on the non-equilibrium reduced probability distribution depending on $\left[\mathit{z}\right]$, y and time t. The inclusion of attractive potentials between the two chains in the corresponding Kramers-like master equation seemed to be, to the best of our knowledge, a novel feature.
• Upon integrating over all $\left[\mathit{z}\right]$ (leaving y unintegrated), the classical reduced Boltzmann equilibrium distribution for the two interacting chains was used as the generating function for an infinite family of orthogonal polynomials in $\left[\mathit{z}\right]$ (depending on y, parametrically). That generating function was not a Gaussian one, due to the various potentials.
• The orthogonal polynomials were used to define non-equilibrium moments (y- and t-dependent) for the non-equilibrium reduced probability distribution.
• The Kramers-like equation yielded a hierarchy for the non-equilibrium moments.
• At temperatures below, but close to, the melting one (i.e., DNA denaturation), the long-time dynamics was approximated by a Smoluchowski-like equation (containing a y-dependent effective potential) for the lowest non-equilibrium moment.
• The application of the mean first passage time formalism [16,17,35] to the Smoluchowski-like equation in step 6 enabled us to study approximately the time duration of non-equilibrium thermal denaturation of the two strands, initially bound to each other, towards configurations of two separate single strands.

We shall firstly outline the contents of the present work by emphasizing the correspondence with the above steps 1–7 in [34] and, secondly, underline the differences with [34]. Section 2 deals with the non-equilibrium evolution of two quantum particles in three spatial dimensions, subject to an attractive potential between them and in contact with a large HB at thermal equilibrium. Throughout our study, we assume that the state of the system, for sufficiently long times, approaches towards its own canonical density operator at thermal equilibrium with the HB, in agreement with the results in [22,23,24,25,26,27,28,29]. We shall simplify the treatment by omitting the detailed analysis of the states of the HB (just retaining the fact that it establishes the temperature), and in an effective way, we will focus exclusively on the quantum states of the two particles. The quantum Hamiltonian, nonequilibrium Wigner function and evolution equation (t-reversible and without ab initio dissipation) for the system, the separation (at both equilibrium and off-equilibrium) of the overall center of mass from the relative motion variables (position x and momentum q) of the two particles and certain spectral properties are treated in Section 2.1, Section 2.2 and Section 2.3 (counterparts of steps 1 and 2 in [34]). The construction of an infinite family of orthogonal polynomials in the relative momentum q, generated by the equilibrium relative Wigner function as the weight function, is carried out in Section 2.4 and Section 2.5 (counterpart of step 3 in [34]). The moments (x- and t-dependent) of the non-equilibrium relative Wigner function are constructed in Section 2.6 (counterpart of step 4 in [34]). Section 3 provides the dynamical equations for the non-equilibrium moments implied by those for the non-equilibrium relative Wigner function. The infinite hierarchy is presented and analyzed in Section 3.1 (counterpart of step 5) in [34]). The restrictions to the lowest non-equilibrium moment and to an irreversible Smoluchowski-like equation for it, under various approximations, are considered and justified in Section 3.2 and Section 3.3 (counterpart of step 6) in [34]). Section 4 deals with approximate kinetic equations. Section 5 compares approximately thermal and chemical equilibria. Section 6 and Section 7 are devoted to the application of the mean first passage time formalism to the Smoluchowski equation for the lowest non-equilibrium moment (counterpart of step 7) in [34]). Section 8 compares briefly our approach with other well-known methods and, in so doing, treats chemical recombination. Section 9 summarizes our conclusions.

Let us quote some key differences between the classical system considered in [34] and the quantum one to be treated here. In the quantum case: (a) even if, in general, Wigner functions could be negative in some domain [36,37], we shall argue that the actual equilibrium Wigner function could still be used as a weight function to generate orthogonal polynomials, by invoking a suitable mathematical framework (Section 2.4); (b) upon constructing the orthogonal polynomials, one integrates over the relative momentum q (while in [34], one integrated over all positions $\left[\mathit{z}\right]$ of the atoms along the two chains); (c) even if the system approaches thermal equilibrium for sufficiently long times [22,23,24,25,26,27,28,29], the nonequilibrium evolution equation for the Wigner function (our starting point in Section 2) does not display either irreversibility or a thermalizing evolution in general (while the Kramers-like master equation in [34] did display them from the outset). Then, approximations (for short thermal wavelengths and long times) are necessary here in order to arrive at the actual irreversible Smoluchowski-like equation for the lowest non-equilibrium moment in Section 3.3. In turn, for a simpler presentation, those approximations underlying the contents of Section 3.3 are carried out successively and separately, for a one-dimensional case, in Appendix C, Appendix D, Appendix E and Appendix F. One underlying basic question is why the starting point of our analysis is the non-equilibrium Wigner equation without explicit dissipation, instead of a master equation displaying the latter from the very outset: Appendix G will be devoted to discuss this issue, a posteriori.

2. Two Particles: Towards Non-Equilibrium Toy Chemical Reactions

2.1. General Features

We shall consider two non-relativistic spinless quantum particles of masses $m_j$ (with $j=1,2$) in three spatial dimensions, inside a large finite volume Ω (outside which, an infinitely repulsive potential is assumed). Those individual particles can be either atoms or small molecules. We emphasize that Ω is large only at the microscopic scale. In Cartesian coordinates, the position vector and the momentum operator of the j-th particle are ${\mathit{x}}_j=(x_{j,1},x_{j,2},x_{j,3})$ and ${\mathit{p}}_j=(p_{j,1},p_{j,2},p_{j,3})$, where $p_{j,\alpha }=-i\hslash (\partial /\partial x_{j,\alpha })$, $\alpha =1,2,3$ and ℏ is Planck’s constant. We shall denote $\mathbf{x}=({\mathit{x}}_1,{\mathit{x}}_2)$.

There are no external forces, and all interactions “seen by the system” are described, in principle, by a time-independent and velocity-independent two-body instantaneous potential V between the particles. The total (Hermitian) quantum Hamiltonian is:

$$H=-\sum _{j=1}^2\sum _{\alpha =1}^3\frac{{\hslash }^2}{2m_j}\frac{{\partial }^2}{\partial x_{j,\alpha }^2}+V\left(\right|x_1-x_2\left|\right).$$

(1)

Let $\mathit{P}={\sum }_{j=1}^2{\mathit{p}}_j$ be the total momentum. Then, one has: $[H,\mathit{P}]=0$ ($[A,B]=AB-BA$ being the commutator of the operators A and B).

We shall separate off the standard CM: its position vector being $\mathit{X}=M^{-1}\left[{\sum }_{j=1}^2{m}_j{x}_j\right]$ with $M={\sum }_{j=1}^2{m}_j$. Let us introduce the relative position vector $\mathit{x}={\mathit{x}}_2-{\mathit{x}}_1$ and the reduced mass $m^{-1}=m_1^{-1}+m_2^{-1}$. Then, one has:

$$H=H_{\mathrm{CM}}+H_{\mathrm{re}},H_{\mathrm{CM}}=-\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial {\mathit{X}}^2},\left[H_{\mathrm{CM}},H_{\mathrm{re}}\right]=0,$$

(2)

$$H_{\mathrm{re}}=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}+V\left(\right|\mathit{x}\left|\right).$$

(3)

Here, $H_{\mathrm{CM}}$ is the free CM Hamiltonian, and $H_{\mathrm{re}}$ is the Hamiltonian of the reduced-mass particle. We now express x in terms of the standard radial ($r=\left|\mathit{x}\right|$) and angular coordinates (θ, φ). Then:

$$H_{\mathrm{re}}=\frac{p_{\mathrm{re},\mathrm{r}}^2}{2m}+V\left(r\right)+\frac{1}{2m{r}^2}{\mathit{l}}^2,p_{\mathrm{re},r}^2=-{\hslash }^2\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial }{\partial r}\right),$$

(4)

where l is the orbital angular momentum operator. Spatial integrations over X and x will be carried out inside large finite volumes Ω. The eigenfunctions corresponding to both $H_{\mathrm{CM}}$ and $H_{\mathrm{re}}$ will vanish, by assumption, at the surfaces enclosing Ω.

2.2. Non-Equilibrium Statistical Formulation and Equilibrium Distributions

Let the two-particle system be immersed in an HB at thermal equilibrium, which is at absolute temperature T. The temporal evolution for $t>0$ is given by the density operator $\rho =\rho \left(t\right)$ (a statistical mixture of quantum states). It fulfills the (t-reversible) operator equation $\partial \rho /\partial t={\left(i\hslash \right)}^{-1}[H,\rho ]$ with initial condition $\rho (t=0)={\rho }_{\mathrm{in}}$. The Hermitian and positive-definite linear operators $\rho \left(t\right)$ and ${\rho }_{\mathrm{in}}$ act on the Hilbert space spanned by the set of all eigenfunctions of H. Unless otherwise stated, we shall not impose that $\rho \left(t\right)$ be normalized.

Let $\beta ={\left(k_{\mathrm{B}}T\right)}^{-1}$, with $k_{\mathrm{B}}$ being Boltzmann’s constant. The canonical density operator describing the two-particle system at thermal equilibrium with the HB is ${\rho }_{\mathrm{eq}}=exp[-\beta H]$ (with $\partial {\rho }_{\mathrm{eq}}/\partial t=0$), and it factorizes as ${\rho }_{\mathrm{eq}}={\rho }_{\mathrm{CM},\mathrm{eq}}\otimes {\rho }_{\mathrm{re},\mathrm{eq}}$, where ${\rho }_{\mathrm{CM},\mathrm{eq}}=exp[-\beta H_{\mathrm{CM}}]$ and ${\rho }_{\mathrm{re},\mathrm{eq}}=exp[-\beta H_{\mathrm{re}}]$.

We now resort to the non-equilibrium Wigner function [1,8,38,39,40,41,42,43], which reads as:

$$W(\mathbf{x},\mathbf{q};t)=\frac{1}{{\left(\pi \hslash \right)}^6}\int {\mathbf{dx}}^{\prime }exp\left[\frac{2i}{\hslash }\sum _{j=1}^2\sum _{\alpha =1}^3{x}_{i,\alpha }^{\prime }{q}_{i,\alpha }\right]\langle \mathbf{x}-{\mathbf{x}}^{\prime }\left|\rho \left(t\right)\right|\mathbf{x}+{\mathbf{x}}^{\prime }\rangle ,$$

(5)

where $\mathbf{x}=(x_1,x_2)$, $\mathbf{q}=(q_1,q_2)$ and ${\mathbf{dx}}^{\prime }={\prod }_{i=1}^2{\prod }_{\alpha =1}^{3}d{x}_{i,\alpha }^{\prime }$.

The time evolution of the Wigner function is furnished by:

$$\frac{\partial W(\mathbf{x},\mathbf{q};t)}{\partial t}=-\sum _{j=1}^2\sum _{\alpha =1}^3\frac{q_{j,\alpha }}{m_j}\frac{\partial W(\mathbf{x},\mathbf{q};t)}{\partial x_{j,\alpha }}+M_{Q}W,$$

(6)

with:

$$\begin{array}{ccc}\hfill M_{Q}W& =& \frac{i}{\hslash {\left(\pi \hslash \right)}^6}\int {\mathbf{dq}}^{\prime }W(\mathbf{x},{\mathbf{q}}^{\prime };t)\int {\mathbf{dx}}^{\prime }exp\left[\frac{2i}{\hslash }{\displaystyle \sum _{j=1}^2}{\displaystyle \sum _{\alpha =1}^3}{x}_{j,\alpha }^{\prime }(q_{j,\alpha }-q_{j,\alpha }^{\prime })\right]\hfill \\ & \times & \left[V\left(\right|x_1+x_1^{\prime }-x_2-x_2^{\prime }\left|\right)-V\left(\right|x_1-x_1^{\prime }-x_2+x_2^{\prime }\left|\right)\right],\hfill \end{array}$$

(7)

where ${\mathbf{dq}}^{\prime }={\prod }_{i=1}^2{\prod }_{\alpha =1}^{3}d{q}_{i,\alpha }^{\prime }$. As Ω is large, we shall approximate spatial integrals by those for an infinite volume when such approximations are harmless, unless some specific discussion is required. Strictly speaking, as Ω is not infinite, the ${\mathbf{q}}^{\prime }$ are discretized momenta, and $\int {\mathbf{dq}}^{\prime }$ in Equation (7) should be interpreted as a six-fold series. However, as Ω is large, we shall disregard the small spacings in ${\mathbf{dq}}^{\prime }$ and approximate it as a six-fold integral (thus, varying continuously). A similar remark and interpretation will apply and be understood whenever integrations over momenta occur (namely, in Equations (10), (16), (18), (19), and so on). We shall accept that all integrals (or all series) over momenta converge for large values of the latter: explicit expressions for $W_{\mathrm{CM},\mathrm{eq}}$ below, together with (11), (12) and the computations in Appendix A will support this assumption. The initial condition $W_{\mathrm{in}}(\mathbf{x},\mathbf{q})$ is determined by ${\rho }_{\mathrm{in}}$ through Equation (5). The equilibrium Wigner function $W_{\mathrm{eq}}(\mathbf{x},\mathbf{q})$ is given by Equation (5) for the actual ${\rho }_{\mathrm{eq}}$. Let $\mathit{Q}={\mathit{q}}_1+{\mathit{q}}_2$ and $\mathit{q}=\left(m_1{\mathit{q}}_2-m_2{\mathit{q}}_1\right)/M$ be the CM and relative momentum vectors, respectively. The actual $W_{\mathrm{eq}}(\mathbf{x},\mathbf{q})$ factorizes as:

$$W_{\mathrm{eq}}(\mathbf{x},\mathbf{q})=W_{\mathrm{CM},\mathrm{eq}}{W}_{\mathrm{re},\mathrm{eq}},$$

(8)

with $W_{\mathrm{re},\mathrm{eq}}=W_{\mathrm{re},\mathrm{eq}}\left(\mathit{x},\mathit{q}\right)$ being the relative equilibrium Wigner function. $W_{\mathrm{CM},\mathrm{eq}}$ is proportional to $exp\left[-\beta ({\mathit{Q}}^2/2M)\right]$. Correspondingly, there is also factorization off-equilibrium: $W(\mathbf{x},\mathbf{q};t)=W_{\mathrm{CM}}(\mathit{X},\mathit{Q};t)W_{\mathrm{re}}\left(\mathit{x},q\mathit{;}t\right)$, with $\partial W_{\mathrm{CM}}/\partial t=-(\mathit{Q}/M)\partial W_{CM}/\partial \mathit{X}$, and the (t-reversible) Wigner function $W_{\mathrm{re}}$ of the relative particle fulfills:

$$\frac{\partial W_{\mathrm{re}}}{\partial t}=-\frac{\mathit{q}}{m}\frac{\partial W_{\mathrm{re}}}{\partial \mathit{x}}+M_{\mathrm{re},Q}{W}_{\mathrm{re}},$$

(9)

$$\begin{array}{ccc}\hfill M_{\mathrm{re},Q}{W}_{\mathrm{re}}& =& \int d^3{\mathit{q}}^{\prime }{W}_{\mathrm{re}}(\mathit{x},{\mathit{q}}^{\prime };t)\int d^3{\mathit{x}}^{\prime }\frac{i}{\hslash {\left(\pi \hslash \right)}^3}exp\left[\frac{i2(\mathit{q}-{\mathit{q}}^{\prime }){\mathit{x}}^{\prime }}{\hslash }\right]\hfill \\ & \times & \left[V\left(\right|\mathit{x}+{\mathit{x}}^{\prime })-V\left(\right|\mathit{x}-{\mathit{x}}^{\prime }\left|\right)\right].\hfill \end{array}$$

(10)

We shall introduce the fixed, physically relevant and x-independent momentum $q_{\mathrm{eq}}={(2m/\beta )}^{1/2}$. Recall that, in the classical case, the equilibrium (or Boltzmann’s) canonical distribution describing thermal equilibrium of two classical particles with an HB is (proportional to a) Gaussian in $\mathbf{q}$: that is, $exp\left[-\beta \left({\sum }_{j=1}^2{\sum }_{\alpha =1}^3{q}_{j,\alpha }^2/\left(2m_j\right)+V\right)\right]$, with $q_{j,\alpha }$ being now the classical momenta. In the high temperature or small β (quasiclassical) regime, Wigner [38] obtained successive approximations for the equilibrium quantum distribution. Based on [38], we shall give the leading terms in the small β expansion for $W_{\mathrm{re},\mathrm{eq}}\left(\mathit{x}=(x_1,x_2,x_3),\mathit{q}=(q_1,q_2,q_3)\right)$:

$$\begin{array}{ccc}\hfill W_{\mathrm{re},\mathrm{eq}}& =& c_0{f}_{\mathrm{re},\mathrm{eq}}\left[1+c_1+c_2\left[{\displaystyle \sum _{\alpha =1}^3}\frac{{\partial }^{2}V}{\partial x_{\alpha }^2}{H}_2(q_{\alpha }/q_{\mathrm{eq}})\right.\right.\hfill \\ & +& \left.\left.{\displaystyle \sum _{\alpha ,{\alpha }^{\prime }=1;\alpha \ne {\alpha }^{\prime }}^3}\frac{{\partial }^{2}V}{\partial x_{\alpha }\partial x_{{\alpha }^{\prime }}}{H}_1(q_{\alpha }/q_{\mathrm{eq}})H_1(q_{{\alpha }^{\prime }}/q_{\mathrm{eq}})\right]\right],\hfill \end{array}$$

(11)

where:

$$f_{\mathrm{re},\mathrm{eq}}=exp\left[-\beta (({\mathit{q}}^2/2m)+V)\right],$$

(12)

$$c_1={\hslash }^2\left[-\frac{{\beta }^2}{12m}\sum _{\alpha =1}^3\frac{{\partial }^{2}V}{\partial x_{\alpha }^2}+\frac{{\beta }^3}{24m}\sum _{\alpha =1}^3{\left(\frac{\partial V}{\partial x_{\alpha }}\right)}^2\right],$$

(13)

and:

$$c_2=\frac{{\beta }^2{\hslash }^2}{48m},$$

(14)

with $c_0={\left(2\pi \hslash \right)}^{-3}$. $H_2$ and $H_1$ denote standard Hermite polynomials. We shall refer to the contributions associated with $c_1$ and $c_2$ as quantum contributions in the quasiclassical regime.

It is a reasonable assumption that the initial conditions ${\rho }_{\mathrm{in}}$ and $W_{\mathrm{in}}(\mathbf{x},\mathbf{q})$ factorize into CM and relative ones. The initial condition for $W_{\mathrm{re}}$ is $W_{\mathrm{re},\mathrm{in}}$.

2.3. Assumptions on $V\left(r\right)$, the Spectrum of $H_{\mathrm{re}}$ and Application to $W_{\mathrm{re},\mathrm{eq}}$

To further advance with the present framework, we suppose that the two-body interaction potential V satisfies the following conditions:

• $V\left(r\right)$ is repulsive ($>0$) for $0\le r<r_0$ (“hard core”, with adequately small $r_0$), attractive ($<0$) in the interval $r_0<r<{\left(3\mathsf{\Omega }/4\pi \right)}^{1/3}$ and vanishes fast as $r\to {\left(3\mathsf{\Omega }/4\pi \right)}^{1/3}$.
• $V\left(r\right)$ is finite everywhere, and its magnitude $\left|V\right|$ is appreciable in $r_0<r<r_0+a<{\left(3\mathsf{\Omega }/4\pi \right)}^{1/3}$. a is understood to be the range of V.
• $V\left(r\right)$ and all $d^{n}V\left(r\right)/d{r}^n$, for $n=1,2,3,\dots$, are continuous for all $r>0$.
• $V\left(r\right)$ does give rise to only one bound state (bound spectrum). Thus, the relevance of the region where $V<0$ is larger than that of the hard core.

As we have mentioned, each individual particle could be either an atom or a small molecule. Hence, on physical grounds, we assume that $V\left(r\right)$ is an effective potential between them (which, in particular, includes and averages over Coulomb interactions). Let ${\phi }_j={\phi }_j\left(\mathit{x}\right)$ denote a suitably-normalized eigenfunction of $H_{\mathrm{re}}$ with corresponding eigenvalue $E_j$ and j denoting a set of labels. $H_{\mathrm{re}}$ has both a discrete spectrum (due to the above assumptions, just one bound state) and a denumerably-infinite number of discrete states. For the discrete spectrum $E_j=E_{\mathrm{d}}<0$ and ${\phi }_j={\phi }_{\mathrm{d}}\left(\mathit{x}\right)$. The denumerably-infinite discrete spectrum has a small spacing, and it becomes a continuous one (sweeping the continuous positive real axis) as ${\mathsf{\Omega }}^{-1}\to 0$. We shall denote it by CS, even if the small ${\mathsf{\Omega }}^{-1}$ remains positive. The eigenfunctions corresponding to the continuous spectrum of $H_{\mathrm{re}}$ are ${\phi }_j={\phi }_{\mathit{k}}\left(\mathit{x}\right)$, with $j\equiv \mathit{k}$ being an almost continuous wavevector, and eigenvalues $E_j=E_{\mathit{k}}\ge 0$. The CS eigenfunctions are normalized through: $({\phi }_{\mathit{k}},{\phi }_{{\mathit{k}}^{\prime }})=\int d^3\mathit{x}{\phi }_{\mathit{k}}^{*}{\phi }_{{\mathit{k}}^{\prime }}={\delta }_{\mathit{k}{\mathit{k}}^{\prime }}$ (a Kronecker delta). Furthermore, $({\phi }_{\mathrm{d}},{\phi }_{\mathrm{d}})=\int d^{3}x{\phi }_{\mathrm{d}}^{*}{\phi }_d=1$ (normalized) and $({\phi }_{\mathrm{d}},{\phi }_{\mathit{k}})=\int d^3\mathit{x}{\phi }_{\mathrm{d}}^{*}{\phi }_{\mathit{k}}=0$. Hence, ${\phi }_{\mathrm{d}}$, and all CS ${\phi }_{\mathit{k}}$ span two separate Hilbert subspaces ${\mathcal{H}}_{\mathrm{d}}$ and ${\mathcal{H}}_{\mathrm{CS}}$. Accordingly, ${\sum }_j$ will include the contribution of both the single discrete eigenfunction plus that of a three-fold infinite summation over the whole CS ones. For the CS only: ${\sum }_j\to (\mathsf{\Omega }/{(2\pi )}^3)\int d^3\mathit{k}$ as ${\mathsf{\Omega }}^{-1}\to 0$. Therefore,

$$\begin{array}{ccc}\hfill W_{\mathrm{re},\mathrm{eq}}\left(\mathit{x},\mathit{q}\right)& =& \frac{1}{{\left(\pi \hslash \right)}^3}\int d^3{\mathit{x}}^{\prime }exp\left(\frac{i2\mathit{q}{\mathit{x}}^{\mathit{\prime }}}{\hslash }\right)\langle \mathit{x}-{\mathit{x}}^{\prime }|exp[-\beta H_{\mathrm{re}}]|\mathit{x}+{\mathit{x}}^{\prime }\rangle \hfill \\ & =& \frac{1}{{\left(\pi \hslash \right)}^3}\int d^3{\mathit{x}}^{\prime }exp\left(\frac{i2\mathit{q}{\mathit{x}}^{\prime }}{\hslash }\right)\sum _{j}exp[-\beta E_j]{\phi }_j(\mathit{x}-{\mathit{x}}^{\prime }){\phi }_j^{*}(\mathit{x}+{\mathit{x}}^{\prime }).\hfill \end{array}$$

(15)

While at $T=0$, there is no transition between the bound state and the CS ones, such a transition is indeed possible for $T>0$ (due to the HB) and plays a key role as $k_{\mathrm{B}}T$ approaches $|E_d|$, as will be discussed below. Thus far, we have implemented the counterparts of steps 1 and 2 in [34]. We emphasize that the contributions of the bound states have disappeared in the high temperature (quasi-classical) regime corresponding to Equations (11)–(14), which should not be employed in regions where bound states be relevant. However, there seems to be no compelling reason for not using Equations (11)–(14) for rough or zeroth order approximations in regions where the contributions due to the bound states are negligible.

2.4. $W_{\mathrm{re},\mathrm{eq}}$ as a Quasi-Definite Functional of Momenta

As $W_{\mathrm{CM},\mathrm{eq}}$ is Gaussian in Q, it gives rise, as a weight function, to an infinite number of orthogonal polynomials in Q: the standard Hermite ones [44]. Furthermore, notice that neither $W_{\mathrm{re},\mathrm{eq}}$, nor $W_{\mathrm{re}}$ can be warranted to be nonnegative in general [38,39,42]. A necessary and sufficient condition for the Wigner function associated with a wave function to be nonnegative is that the latter be a Gaussian distribution [36]. However, the domain in which $W_{\mathrm{re},\mathrm{eq}}<0$ may occur cannot be large and has to be consistent with the fact that the marginals of both $W_{\mathrm{re},\mathrm{eq}}$ and $W_{\mathrm{re}}$ are nonnegative.

Let us now invoke, in outline, a mathematical framework, based on the theory of orthogonal polynomials [45], the physical interest of which will be appreciated below. For simplicity, we shall focus on the one-dimensional case, leaving out the direct extension to three dimensions. Let us consider a kernel $K=K\left(y\right)$ (which could be negative), a set of functions $f=f\left(y\right)$ and the following functionals $L_K$ determined by the kernel K: $f\to L_K\left[f\right]={\int }_{-\infty }^{+\infty }dyK\left(y\right)f\left(y\right)$. We will assume that all integrals over y are convergent. Let ${\mu }_n=L_K\left[y^n\right]$, $n=0,1,2,3,\dots$, and consider the set of all $(S+1)\times (S+1)$ matrices $M_S$ (where $S=0,1,2,3,\dots$) with the $(i,j)$-th element equal to ${\mu }_{i+j}$ (where $i,j=0,...,S$) and their respective determinants: $det\left(M_S\right)$. By definition, the functional $L_K$ is quasi-definite if $det\left(M_S\right)\ne 0$ for any $S=0,1,2,3,\dots$ [45]. If $L_K$ is a quasi-definite functional, we resort to a theorem in [45], which implies the following: (i) the existence of an infinite family of orthogonal polynomials, denoted here as $H_{K,n}=H_{K,n}\left(y\right)$, for $n=0,1,2,\dots$, with weight function K (even if $K<0$ in some domain of y); and (ii) the $H_{K,n}$’s satisfy a basic three-term linear recurrence relation, omitted here.

Let $\mathit{y}=\mathit{q}/q_{\mathrm{eq}}$. The argument in Appendix A supports (albeit, it does not prove) that $W_{\mathrm{re},\mathrm{eq}}$ determines a quasi-definite functional of y for any x. Based on that argument, we shall assume henceforth that, regarding their $\mathit{q}$-dependences, the Wigner function $W_{\mathrm{re},\mathrm{eq}}$ determines a quasi-definite functional $L_{W_{\mathrm{re},\mathrm{eq}}}$ (for any x and t). The interest of this assumption is obvious: if it holds (as we shall suppose), it implies the existence of an infinite family of orthogonal polynomials $H_{\left[n\right]}\left(\mathit{y}\right)$ generated by $W_{\mathrm{re},\mathrm{eq}}$ (Section 2.5). It is possible, but unnecessary, to extend the assumption to $W_{\mathrm{re}}$.

2.5. Orthogonal Polynomials $H_{\left[n\right]}$ Generated by $W_{\mathrm{re},\mathrm{eq}}$

Let us define $\left[n\right]=(n_1,n_2,n_3)$, with the $n_{\alpha }$ being non-negative integers ($\alpha =1,2,3$). Use will also be made of the additional shorthand notations: $\left[0\right]=(0,0,0)$, $\left[1_{\alpha }\right]$, which represents a vector with a one in the α position with the rest of the components being zero, and $\left[2_{\alpha }\right]$ similar to $\left[1_{\alpha }\right]$, but with a two in the α position. We introduce (unnormalized) orthogonal polynomials $H_{\left[n\right]}$ in y determined by $W_{\mathrm{re},\mathrm{eq}}$, which acts as a (in general, non-Gaussian) weight function. We choose $H_{\left[0\right]}\left(\mathit{y}\right)=1$. We shall also employ additional useful notations, like: $H_{\left[1_{\alpha }\right]}=H_{\left[1_{\alpha }\right]}\left(\mathit{y}\right)=y_{\alpha }$, $H_{\left[1_{\alpha }\right],\left[1_{\gamma }\right]}=H_{\left[1_{\alpha }\right]}.H_{\left[1_{\gamma }\right]}$ for $\alpha \ne \gamma$ and $H_{\left[2_{\alpha }\right]}=H_{\left[2_{\alpha }\right]}\left(\mathit{y}\right)=y_{\alpha }^2+{ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$. In general, the $H_{\left[n\right]}$’s are constructed recurrently as follows. We impose for $\left[n\right]\ne \left[n^{\prime }\right]$, and any x (left unintegrated), that:

$$\int d^3\mathit{y}{W}_{\mathrm{re},\mathrm{eq}}\left(\mathit{x},\mathit{y}\right)H_{\left[n\right]}\left(\mathit{y}\right)H_{\left[n^{\prime }\right]}\left(\mathit{y}\right)=0,$$

(16)

where:

$$H_{\left[n\right]}\left(\mathit{y}\right)=y_1^{n_1}{y}_2^{n_2}{y}_3^{n_3}+\sum _{\left[j\right]}{ϵ}_{\left[n\right],[n-j]}{y}_1^{n_1-j_1}{y}_2^{n_2-j_2}{y}_3^{n_3-j_3}+\cdots$$

(17)

$0\le j_{\alpha }\le n_{\alpha }$ and ${\sum }_{\alpha =1}^3(n_{\alpha }-j_{\alpha })=2,4,\dots$. Here, ${ϵ}_{\left[n\right],[n-j]}$ are dimensionless and y-independent (though x-dependent, in general). One has ${ϵ}_{\left[n\right],[n-j]}=0$ if ${\sum }_{\alpha =1}^3{j}_{\alpha }$ is odd, so that $H_{\left[n\right]}(-\mathit{y})={(-1)}^{{\sum }_{\alpha =1}^3{n}_{\alpha }}{H}_{\left[n\right]}\left(\mathit{y}\right)$. The ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$’s are given and estimated in Appendix B.

The orthonormalized polynomials are $H_{\left[n\right]}\left(\mathit{y}\right)/{\left(h_{\left[n\right]}\right)}^{1/2}$, with a (x-dependent) normalization factor $h_{\left[n\right]}$ defined through:

$$h_{\left[n\right]}\left(\mathit{x}\right)\equiv \int d^3\mathit{y}{W}_{\mathrm{re},\mathrm{eq}}\left(\mathit{x},\mathit{y}\right)H_{\left[n\right]}^2\left(\mathit{y}\right).$$

(18)

Expressing the $H_{\left[n\right]}$’s as sums of products of standard Hermite polynomials is also possible, but less convenient here. The orthogonal polynomials determined by the classical relative Boltzmann distribution $f_{\mathrm{re},\mathrm{eq}}$ in (12) are proportional to products of the standard Hermite polynomials $H_n\left(y_{\alpha }\right)$ [44]. Then, for the latter, the coefficients in the (classical) counterpart of (17) can be directly obtained and are x-independent. Section 2.4 and the present one have implemented the counterpart of step 3 in [34].

2.6. Moments: Off-Equilibrium and at Equilibrium

We now turn to the non-equilibrium moments. They are defined through:

$$\begin{array}{c}\hfill W_{\left[n\right]}\equiv W_{\left[n\right]}(\mathit{x};t)=\int d^3\mathit{y}{H}_{\left[n\right]}\left(\mathit{y}\right)W_{\mathrm{re}}\left(\mathit{x},\mathit{y};t\right).\end{array}$$

(19)

When $W_{\mathrm{re}}=W_{\mathrm{re},\mathrm{eq}}$, all equilibrium moments with $\left[n\right]\ne \left[0\right]$ vanish, except:

$$\begin{array}{ccc}\hfill W_{\mathrm{eq},\left[0\right]}\left(\mathit{x}\right)& =& \int d^3\mathit{y}{W}_{\mathrm{re},\mathrm{eq}}(\mathit{x},\mathit{y})=\frac{1}{q_{\mathrm{eq}}^3}{\displaystyle \sum _j}exp[-\beta E_j]{\phi }_j\left(\mathit{x}\right){\phi }_j^{*}\left(\mathit{x}\right)\hfill \\ \hfill & =& \frac{1}{q_{\mathrm{eq}}^3}exp\left(-\beta E_{\mathrm{d}}\right){\phi }_{\mathrm{d}}\left(\mathit{x}\right){\phi }_{\mathrm{d}}^{*}\left(\mathit{x}\right)+W_{\mathrm{eq},CS,\left[0\right]}.\hfill \end{array}$$

(20)

$W_{\mathrm{eq},CS,\left[0\right]}$ is the CS contribution. This justifies, a posteriori, the interest of having introduced the $H_{\left[n\right]}$’s and the $W_{\mathrm{eq},\left[n\right]}$’s. In the t-evolution, not far from thermal equilibrium, and recalling [22,23,24,25,26,27,28,29], one expects that $W_{\left[0\right]}$ would be dominant and that any $W_{\left[n\right]}$, for any $\left[n\right]\ne \left[0\right]$, be small (the more negligible the larger t and ${\sum }_{\alpha =1}^3{n}_{\alpha }$ are). The initial condition $W_{\mathrm{in},\left[n\right]}$ for $W_{\left[n\right]}$ is obtained by replacing $W_{\mathrm{re}}$ by $W_{\mathrm{re},\mathrm{in}}$ in Equation (19). $W_{\mathrm{re}}$ and any $W_{\left[n\right]}$ have the same dimension as ${\hslash }^{-3}$. This section has implemented the counterpart of step 4 in [34].

3. Hierarchy for the Non-Equilibrium Moments and Approximations

3.1. Non-Equilibrium Hierarchy

The transformation of Equations (9), (10) and (19) into an infinite linear hierarchy for the nonequilibrium moments $W_{\left[n\right]}$ can be carried out through direct computations and cancellations, increasingly cumbersome as ${\sum }_{\alpha =1}^3{n}_{\alpha }$ grows. The equations for the lowest order moments are:

$$\frac{\partial W_{\left[0\right]}}{\partial t}=-\frac{q_{\mathrm{eq}}}{m}\sum _{\alpha =1}^3\frac{\partial }{\partial x_{\alpha }}{W}_{\left[1_{\alpha }\right]},$$

(21)

$$\frac{\partial W_{\left[1_{\alpha }\right]}}{\partial t}=-\frac{q_{\mathrm{eq}}}{m}\frac{\partial }{\partial x_{\alpha }}{W}_{\left[2_{\alpha }\right]}-\frac{q_{\mathrm{eq}}}{m}\sum _{\gamma =1,\gamma \ne \alpha }^3\frac{\partial }{\partial x_{\gamma }}{W}_{\left[1_{\alpha }\right],\left[1_{\gamma }\right]}-M_{\left[1_{\alpha }\right],\left[0\right]}{W}_{\left[0\right]},$$

(22)

where:

$$M_{\left[1_{\alpha }\right],\left[0\right]}{W}_{\left[0\right]}=-\frac{q_{\mathrm{eq}}}{m}\frac{\partial }{\partial x_{\alpha }}\left({ϵ}_{\left[2_{\alpha }\right],\left[0\right]}{W}_{\left[0\right]}\right)+\frac{1}{q_{\mathrm{eq}}}\frac{\partial V}{\partial x_{\alpha }}{W}_{\left[0\right]}.$$

(23)

The general (t-reversible) hierarchy implied by Equations (9) and (10) for any $\left[n\right]$ is:

$$\frac{\partial W_{\left[n\right]}}{\partial t}=-\sum _{\gamma =1}^3{M}_{\left[n\right],\left[{(n+1)}_{\gamma }\right]}{W}_{\left[{(n+1)}_{\gamma }\right]}-\sum _{\left[n^{\prime }\right]=1}^{\left[n\right]}{M}_{\left[n\right],[n-n^{\prime }]}{W}_{[n-n^{\prime }]},$$

(24)

with:

$$M_{\left[n\right],\left[{(n+1)}_{\gamma }\right]}{W}_{{[n+1]}_{\gamma }}\equiv \frac{q_{\mathrm{eq}}}{m}\frac{\partial W_{{[n+1]}_{\gamma }}}{\partial x_{\gamma }}.$$

(25)

Having in mind that $W_{\left[n\right]}=W_{n_1,n_2,n_3}$, the subscript ${[n+1]}_{\gamma }$ in $W_{\left[{(n+1)}_{\gamma }\right]}$ denotes successively $W_{n_1+1,n_2,n_3}$, $W_{n_1,n_2+1,n_3}$ and $W_{n_1,n_2,n_3+1}$, and analogously for the corresponding interpretation of the sum ${\sum }_{\left[n^{\prime }\right]=1}^{\left[n\right]}{M}_{\left[n\right],[n-n^{\prime }]}{W}_{[n-n^{\prime }]}$ in (24) (which we do not make explicit to avoid unnecessary notational complications). Here, $M_{\left[n\right],\left[{(n+1)}_{\gamma }\right]}$, and any $M_{\left[n\right],[n-n^{\prime }]}$ come from the first and the second terms on the right-hand side of Equation (9), respectively. We shall omit, for simplicity, the $M_{\left[n\right],[n-n^{\prime }]}$’s. It is simpler to give and discuss the counterparts of the $M_{\left[n\right],[n-n^{\prime }]}$’s in the one-dimensional case (see Appendix C for details). The structure of the hierarchy (24) is a genuine consequence of quantum mechanics. All coefficients in (24) are expressed in terms of V and of quantities computed out of the equilibrium solution $W_{\mathrm{re},\mathrm{eq}}$.

Let $W_{\mathrm{re}}=W_{\mathrm{re},\mathrm{eq}}$. By using the expression for ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$, given later in Equations (B1) and (B2), together with Equations (20) and (23), we obtain that:

$$M_{\left[1_{\alpha }\right],\left[0\right]}{W}_{\mathrm{eq},\left[0\right]}=0.$$

(26)

Moreover, $W_{\mathrm{eq},\left[0\right]}(\ne 0)$ and $W_{\mathrm{eq},\left[n\right]}=0$ for $\left[n\right]\ne \left[0\right]$ solve the hierarchy (24), with all $\partial W_{\mathrm{eq},\left[n\right]}/\partial t=0$. We shall treat briefly the following stationary solution of Equation (24): $W_{\left[n\right]}=0$ for $\left[n\right]\ne \left[0\right]$, $n=1,2,3,...$ and $W_{\left[0\right]}\equiv W_{\left[0\right],\mathrm{st}}$, with $\partial W_{\left[0\right],\mathrm{st}}/\partial t=0$. With the ${ϵ}_{[2_{\alpha },\left[0\right]}$ (determined by $W_{\mathrm{eq}}$) used so far, it follows that: $M_{\left[1_{\alpha }\right],\left[0\right]}{W}_{\left[0\right],\mathrm{st}}=0$; hence, a direct computation by recalling Equation (26) shows that, consistently, $W_{\left[0\right],\mathrm{st}}=W_{\mathrm{eq},\left[0\right]}$. This section has implemented the counterpart of step 5 in [34].

3.2. Order of Magnitude Estimates

From now on, we shall consider values typical of microscopic scales and phenomena. The range a of V is a few Å (for instance, 3 up to 12 Å). We shall write the reduced mass as $m=n\times m_{\mathrm{ne}}$, with $m_{\mathrm{ne}}$ and n being, respectively, the neutron mass and a positive (approximate) integer (not to be confused with the n’s appearing as subscripts in polynomials and moments). If the individual quantum particles are atoms, then n is less than about one half of the largest atomic mass number (say, n is less than about 125). If the individual quantum particles are molecules, then n can be $\gg 125$. We shall introduce the length scale $\delta x=a/l$, where $l>2$, characterizing approximately the smallest scale of appreciable variations for V. For instance, $\delta x=(1/5)a$, so that 0.6 Å$<\delta x<2.5$ Å. Moreover, the magnitude $V_0$ of $\left|V\right|$, averaged over the domain in which $V<0$, lies between one and 10 electron volts. Notice that $V_0$ has a similar order of magnitude as the energy $|E_{\mathrm{d}}|$ for the bound state (see Section 2.3). $\delta V(<V_0)$ will denote the average variation of V within its range (where $V\ne 0$) in a scale $\delta x$. Similarly, ${\delta }^{n}V$, $n=2,3,...$, will denote the average variation of ${\delta }^{n-1}V$ in a scale $\delta x$ (${\delta }^{1}V\equiv \delta V$). For estimates, one could regard that $\left|\delta V\right|$ be about one order of magnitude smaller than $V_0$ and that $|{\delta }^{n}V|$ be about one order of magnitude smaller than $|{\delta }^{n-1}V|$, $n=2,3,...$.

Let ${\lambda }_{\mathrm{th}}=\hslash {[\beta /2m]}^{1/2}=\hslash /q_{\mathrm{eq}}$ be some suitable thermal wavelength. We shall assume that n and T are such that the following conditions are fulfilled: (a) ${\lambda }_{\mathrm{th}}$ is smaller than $\delta x$, which sets a lower limit on $nT$; and (b) $(3/2)k_{\mathrm{B}}T$ is smaller than or at most $V_0$, which sets an upper limit on T. Both (a) and (b) correspond to the regime of typical chemical reactions (that is, a quantum regime, but with small ${\lambda }_{\mathrm{th}}$). Thus, even if the relative particle is in the quantum regime (but not in the classical high-temperature one), ${\lambda }_{\mathrm{th}}$ is, on average, smaller than $\delta x$ and a.

By using computations in Appendix B, and as zeroth order approximations, we shall take ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}\simeq -(1/2)$ for $\left|\mathit{x}\right|>a+r_0$ (dominated by the CS and, in turn, approximated through the leading term $c_0{f}_{\mathrm{re},\mathrm{eq}}$ in (11)) and ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}\simeq -({\lambda }_{\mathrm{th}}^2/{\left(\delta x\right)}^2)$ for $r_0<\left|\mathit{x}\right|<a+r_0$ (dominated by the bound state), both of them independent of the index α and of the angular coordinates. We shall disregard ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$ and approximations for it within the hard core domain $\left|\mathit{x}\right|<r_0$ (and, more generally, contributions from the latter domain), because they will be neither required nor relevant.

All operators to be considered (namely, $M_{\left[n\right],\left[{(n+1)}_{\alpha }\right]}$ and $M_{\left[n\right],[n-n^{\prime }]}$) below in this section, which act on the $W_{\left[n\right]}$’s in Equation (24), have dimension (time)^‒1. Their orders of magnitude (denoted as ${\left({\tau }^{*}\right)}^{-1}$) will now be estimated. Their inverses (namely, the ${\tau }^{*}$’s) can be regarded as effective evolution times, and henceforth, we will refer to them as such.

The order of magnitude involved upon applying $(q_{\mathrm{eq}}/m)(\partial /\partial x_{\alpha })$ (that is, $M_{\left[n\right],\left[{(n+1)}_{\alpha }\right]}$), is of the order of $\hslash /\left({\lambda }_{\mathrm{th}}m\delta x\right)$. The operator $M_{\left[1_{\alpha }\right],\left[0\right]}$ in (23) gives rise to: (i) $(1/q_{\mathrm{eq}})(\partial V/\partial x_{\alpha })$, which has an order of magnitude about $({\lambda }_{\mathrm{th}}/\hslash \delta x)\delta V$ for $r_0<\left|\mathit{x}\right|<r_0+a$ and fully negligible for $\left|\mathit{x}\right|>r_0+a$; (ii) $(q_{\mathrm{eq}}/m){ϵ}_{\left[2_{\alpha }\right],\left[0\right]}(\partial /\partial x_{\alpha })$, which has an order of magnitude about ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}\hslash /\left({\lambda }_{\mathrm{th}}m\delta x\right)$, where, in turn, ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$ has been estimated above for $r_0<\left|\mathit{x}\right|<r_0+a$ and for $\left|\mathit{x}\right|>r_0+a$; (iii) $(q_{\mathrm{eq}}/m)(\partial {ϵ}_{\left[2_{\alpha }\right],\left[0\right]}/\partial x_{\alpha })$, which has an order of magnitude of (at most) $\hslash /\left({\lambda }_{\mathrm{th}}m\delta x\right){ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$ for $r_0<\left|\mathit{x}\right|<r_0+a$ and negligible for $\left|\mathit{x}\right|>r_0+a$. In practice, the order of magnitude involved upon applying $M_{\left[1_{\alpha }\right],\left[0\right]}$ is the sum of (i) plus (ii) (since that of (iii) will be neglected). Notice that the above estimate of $M_{\left[1_{\alpha }\right],\left[0\right]}$ depends on $\left|\mathit{x}\right|$.

We now turn to the orders of magnitude involved upon applying $M_{\left[n\right],[n-n^{\prime }]}$ in (24). The orders of magnitude of $M_{\left[n\right],[n-1]}$ can be expected to be respectively similar to the various contributions of $M_{\left[1_{\alpha }\right],\left[0\right]}$. On a similar basis, they can be justified by considering explicit expressions in Equations (C5)–(C9) for the counterparts of $M_{\left[1_{\alpha }\right],\left[0\right]}$ and of the $M_{\left[n\right],[n-1]}$’s in the one-dimensional case (Appendix C). On the other hand, the specific quantum contributions $M_{\left[n\right],[n-n^{\prime }]}$, $n^{\prime }>1$ have structures proportional to $({\hslash }^{2n^{\prime \prime }}/q_{\mathrm{eq}}^{2n^{\prime \prime }+1})$ times spatial partial derivatives of V of order $2n^{\prime \prime }+1$ (for various $n^{\prime \prime }$, $2n^{\prime \prime }+1$ growing as $2^{-1}{\sum }_{\alpha =1}^3{n}_{\alpha }-1$), as exemplified through Equations (C2) and (C9). Such structures yield contributions smaller than that of $(1/q_{\mathrm{eq}})(\partial V/\partial x_{\alpha })$ by factors ${({\lambda }_{\mathrm{th}}/\delta x)}^{2n^{\prime \prime }}\times ({\delta }^{2n^{\prime \prime }}V/\delta V)$, by virtue of Assumption (a) above. Moreover, one could expect $\left|\right({\delta }^{2n^{\prime \prime }}V/\delta V\left)\right|$ to be smaller than unity and to decrease as n grows. We anticipate that the quantum contributions $M_{\left[n\right],[n-n^{\prime }]}$, $n^{\prime }>1$, will be neglected from now on (see the analysis in Appendix C and Appendix D, in one spatial dimension). The main conclusion is that, under the above Assumption (a), the overall order of magnitude involved upon applying $M_{\left[n\right],[n-1]}$ in (24) is about the same as for $M_{\left[1_{\alpha }\right],\left[0\right]}$. Namely, the sum of (i) plus (ii) plus (iii).

There are two other effective evolution times respectively associated with the operator $D[\left[1_{\alpha }\right];ϵ]$ (Section 3.3 and Appendix D and Appendix E) and to the mean first passage time (Section 6). After having introduced the latter two, the estimates for all those effective evolution times will be compared and discussed a posteriori in Section 7.

The quantum contributions due to (13) and (14) in the quasiclassical regime, for $\left|\mathit{x}\right|$ outside the region where the bound state is concentrated, can be regarded as negligible compared to unity in Equation (11), as V decreases very fast there.

3.3. Approximations: Small Thermal Wavelength and Long Time

By using the assumption and the order of magnitude estimates in Section 3.2, one analyzes the non-equilibrium hierarchy (24) for the $W_{\left[n\right]}$’s and proceeds to the small thermal wavelength quantum regime (STWQR). In the latter, x varies by units of order $\delta x>{\lambda }_{\mathrm{th}}$. Any $M_{\left[n\right],{[n-n^{\prime }]}_{\gamma }}$, with $n^{\prime }>1$, yields a contribution smaller than that from $M_{\left[n\right],{[n-1]}_{\gamma }}$, as explained in Section 3.2. Then, by neglecting all $M_{\left[n\right],{[n-n^{\prime }]}_{\gamma }}$’s with $n^{\prime }>1$, Equation (24) becomes the approximate (t-reversible) three-term hierarchy:

$$\frac{\partial W_{\left[n\right]}}{\partial t}=-\sum _{\gamma =1}^3{M}_{\left[n\right],\left[{(n+1)}_{\gamma }\right]}{W}_{\left[{(n+1)}_{\gamma }\right]}-\sum _{\gamma =1}^3{M}_{\left[n\right],[n-1]}{W}_{[n-1]}.$$

(27)

See Appendix D for explanations of the argument yielding (27) in the one-dimensional case. We emphasize that (27) is still a quantum hierarchy, since it contains all $M_{\left[n\right],[n-1]}$, which are quantum-mechanical. For simplicity (and without essential loss of generality), we shall assume the initial condition $W_{\mathrm{in},\left[0\right]}\ne 0$, $W_{\mathrm{in},\left[0\right]}\ne W_{\mathrm{eq},\left[0\right]}$ and $W_{\mathrm{in},\left[n\right]}=0$ for $\left[n\right]\ne \left[0\right]$. The solution of the resulting approximate non-equilibrium hierarchy (27) is given, through a Laplace transform from t to the variable s, in terms of products of certain s-dependent generalized operator continued fractions $D\left[\right[n];s]$, which satisfy the three-dimensional generalization of the one-dimensional Equation (D3). Next, we proceed to the long-time approximation (thus introducing irreversibility) for t longer than a certain largest effective evolution time ${\tau }_{\max }^{*}$, to be estimated in Section 7, and also by following the one-dimensional case of Appendix D and Appendix E (in particular, Equation (E1)). One gets (for fixed and small $s=ϵ>0$) the following irreversible approximate quantum (Smoluchowski-like) equation for the lowest non-equilibrium moment $W_{\left[0\right]}(x;t)$:

$$\frac{\partial W_{\left[0\right]}}{\partial t}=\frac{q_{\mathrm{eq}}}{m}\sum _{\alpha =1}^3\frac{\partial }{\partial x_{\alpha }}D[\left[1_{\alpha }\right];ϵ]\left[M_{\left[1_{\alpha }\right],\left[0\right]}\right]W_{\left[0\right]},$$

(28)

with $D[\left[1_{\alpha }\right];ϵ]$ being a quantum operator and assuming approximately the initial condition $W_{\mathrm{in},\left[0\right]}$. Notice that ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$, which is contained in $M_{\left[1_{\alpha }\right],\left[0\right]}]$, also displays quantum effects (compare with Equation (E1)). Our Equation (28) embodies stochasticity and displays a structure typical of diffusion-convection-reaction equations: it is linear; convection is determined by the V-dependent term; and external sources are absent in it. Section 3.2 and the present section have implemented the counterpart of step (6) in [34].

As a side remark, Equation (28) can be recast, by using (26), into a form in which only $D[\left[1_{\alpha }\right];ϵ]$, statistical equilibrium averages and $W_{\mathrm{eq},\left[0\right]})$ appear:

$$\frac{\partial W_{\left[0\right]}}{\partial t}={\left(\frac{q_{\mathrm{eq}}}{m}\right)}^2\sum _{\alpha =1}^3\frac{\partial }{\partial x_{\alpha }}D[\left[1_{\alpha }\right];ϵ]\left(\frac{\partial \left({\langle y_{\alpha }^2\rangle }_{\mathrm{eq}}{W}_{\left[0\right]}\right)}{\partial x_{\alpha }}-\frac{\partial ln\left({\langle y_{\alpha }^2\rangle }_{\mathrm{eq}}{W}_{\mathrm{eq},\left[0\right]}\right)}{\partial x_{\alpha }}{\langle y_{\alpha }^2\rangle }_{\mathrm{eq}}{W}_{\left[0\right]}\right).$$

(29)

The computation of ${\langle y_{\alpha }^2\rangle }_{\mathrm{eq}}$ and $D[\left[1_{\alpha }\right];ϵ]$ in Equation (29) does require quantum mechanics.

4. Towards Kinetic Equations

In what follows, we will consider and interpret Equation (28) by replacing approximately the operator $D[\left[1_{\alpha }\right];ϵ]$ by a positive and α-independent constant D, by assuming that ${ϵ}_{\left[2_{\alpha }\right],0}={ϵ}_{\left[2\right],0}\left(r\right)$ does not depend on α nor on the angular variables and by supposing that $W_{\left[0\right]}$ equals $W_{\left[0\right]}(r,t)$ (the angular variables are absent). Such Equation (28), with the latter approximation, will be equivalent to Equation (46) in Section 6. Let $r_1<{\mathsf{\Omega }}^{1/3}$ be some radial coordinate just slightly larger than $r_0+a$, so that V is negligible at $r_1$ and, of course, for $r>r_1$. Let $S\left(r_1\right)$ be a spherical surface with center at $r=0$ and radius $r_1$. Let $\mathit{n}=(n_1,n_2,n_3)$ be the outer normal unit vector at a generic point on $S\left(r_1\right)$, and let $\mathsf{\Omega }\left(r_1\right)$ be the sphere enclosed by $S\left(r_1\right)$. First, we shall integrate Equation (28) over $\mathsf{\Omega }\left(r_1\right)$ (at fixed t) and apply Green’s theorem so as to transform the volume integral on the resulting right-hand side into a surface integral over $S\left(r_1\right)$ (with differential surface element $dS\left(r_1\right)$). The t-dependence will not be written explicitly. We have:

$$\begin{array}{ccc}\hfill \frac{\partial }{\partial t}{\int }_{\mathsf{\Omega }\left(r_1\right)}{d}^3\mathit{x}{W}_{\left[0\right]}& =& \frac{q_{\mathrm{eq}}D}{m}{\int }_{S\left(r_1\right)}dS\left(r_1\right){\displaystyle \sum _{\alpha =1}^3}{n}_{\alpha }\left[M_{\left[1_{\alpha }\right],\left[0\right]}\right]W_{\left[0\right]}\hfill \\ & \simeq & -\frac{q_{\mathrm{eq}}^{2}D}{m^2}{\int }_{S\left(r_1\right)}dS\left(r_1\right){\left[\frac{\partial }{\partial r}\left[{ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right]\right]}_{r_1}\hfill \\ & \simeq & -\frac{q_{\mathrm{eq}}^{2}D}{m^2\mathsf{\Delta }r}{\int }_{S\left(r_1\right)}dS\left(r_1\right)\left[\left[{ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right]\left(r_1+\frac{\mathsf{\Delta }r}{2}\right)-\left[{ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right]\left(r_1-\frac{\mathsf{\Delta }r}{2}\right)\right],\hfill \end{array}$$

(30)

where ${\sum }_{\alpha =1}^3{n}_{\alpha }[\partial /\partial x_{\alpha }]=\partial /\partial r$, the subscript $r_1$ denotes the derivative at $r=r_1$, and we have neglected the contribution of V. In the rightmost side of Equation (30), we have replaced $\partial /\partial r{\left[{ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right]}_{r_1}$ by ${\left(\mathsf{\Delta }r\right)}^{-1}[\left[{ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right](r_1+2^{-1}\mathsf{\Delta }r)-\left[{ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right](r_1-2^{-1}\mathsf{\Delta }r)]$ with small $\mathsf{\Delta }r$ (such that $r_1-2^{-1}\mathsf{\Delta }r<r_0+a$). Second, we shall integrate Equation (28) over the volume (denoted as $\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right)$) corresponding to $r_1\le r\le {[\left(3\mathsf{\Omega }\right)/\left(4\pi \right)]}^{1/3}$. The latter volume is enclosed by $S\left(r_1\right)$ and by the spherical surface (denoted as $S\left(\mathsf{\Omega }\right)$) with center at $r=0$ and radius ${[\left(3\mathsf{\Omega }\right)/\left(4\pi \right)]}^{1/3}$ (also at fixed t). We shall also apply Green’s theorem so as to transform the latter volume integral into two surface integrals over $S\left(r_1\right)$ and $S\left(\mathsf{\Omega }\right)$. Furthermore, let $\mathit{n}\left(\mathsf{\Omega }\right)=(n_1\left(\mathsf{\Omega }\right),n_2\left(\mathsf{\Omega }\right),n_3\left(\mathsf{\Omega }\right))$ be the outer normal unit vector at a generic point of $S\left(\mathsf{\Omega }\right)$. Then, by proceeding as in the derivation of Equation (30), we get:

$$\frac{\partial }{\partial t}{\int }_{\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right)}{d}^3\mathit{x}{W}_{\left[0\right]}\simeq \frac{q_{\mathrm{eq}}^{2}D}{m^2}{\int }_{S\left(r_1\right)}dS\left(r_1\right){\left[\frac{\partial }{\partial r}\left({ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right)\right]}_{S\left(r_1\right)}+I\left(\mathsf{\Omega }\right),$$

(31)

$$I\left(\mathsf{\Omega }\right)=-\frac{q_{\mathrm{eq}}^{2}D}{m^2}{\int }_{S\left(\mathsf{\Omega }\right)}dS\left(\mathsf{\Omega }\right){\left[\frac{\partial }{\partial r}\left({ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right)\right]}_{S\left(\mathsf{\Omega }\right)}.$$

(32)

For large Ω, by neglecting $I\left(\mathsf{\Omega }\right)$ and approximating like in Equation (30), Equation (31) becomes:

$$\begin{array}{c}\frac{\partial }{\partial t}{\int }_{\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right)}{d}^3\mathit{x}{W}_{\left[0\right]}\simeq \hfill \\ \frac{q_{\mathrm{eq}}^{2}D}{m^2\mathsf{\Delta }r}{\int }_{S\left(r_1\right)}dS\left(r_1\right)\left[\left[{ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right](r_1+2^{-1}\mathsf{\Delta }r)-\left[{ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right](r_1-2^{-1}\mathsf{\Delta }r)\right].\hfill \end{array}$$

(33)

Notice that, consistently:

$$\frac{\partial }{\partial t}{\int }_{\mathsf{\Omega }\left(r_1\right)}{d}^3\mathit{x}{W}_{\left[0\right]}+\frac{\partial }{\partial t}{\int }_{\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right)}{d}^3\mathit{x}{W}_{\left[0\right]}\simeq 0,$$

(34)

which supports a posteriori disregarding $I\left(\mathsf{\Omega }\right)$. For adequately large t, equilibrium sets in and the right-hand sides of the approximate Equations (33) and (30) vanish. We shall analyze the resulting equilibrium equation in Section 5.

Next, assume for the time being that, for adequately large t and as rough approximations: (i) for x inside $\mathsf{\Omega }\left(r_1\right)$, spatial variations of $W_{\left[0\right]}$ are not significant, and the latter is dominated by the contribution of the bound state; so that ${\int }_{\mathsf{\Omega }\left(r_1\right)}{d}^3\mathit{x}{W}_{\left[0\right]}\simeq \mathsf{\Omega }\left(r_1\right)W_{\left[0\right]}(r_1-2^{-1}\mathsf{\Delta }r)$; (ii) for x inside $\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right)$, spatial variations of $W_{\left[0\right]}$ are not important, and the latter is dominated by the contribution of the CS, so that ${\int }_{\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right)}{d}^3\mathit{x}{W}_{\left[0\right]}\simeq (\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right))W_{\left[0\right]}(r_1+2^{-1}\mathsf{\Delta }r)$. The volume of the sphere $\mathsf{\Omega }\left(r_1\right)$ has been denoted as $\mathsf{\Omega }\left(r_1\right)$, and so on, for $\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right)$. Then, Equations (33) and (30) become two approximate linearly-coupled kinetic equations for $W_{\left[0\right]}(r_1-2^{-1}\mathsf{\Delta }r)$ and $W_{\left[0\right]}(r_1+2^{-1}\mathsf{\Delta }r)$,

$$\begin{array}{ccc}\hfill \frac{\partial W_{\left[0\right]}(r_1-2^{-1}\mathsf{\Delta }r)}{\partial t}& \simeq & -\frac{q_{\mathrm{eq}}^{2}D}{m^2\mathsf{\Delta }r\mathsf{\Omega }\left(r_1\right)}{\int }_{S\left(r_1\right)}dS\left(r_1\right)\left[{ϵ}_{\left[2\right],\left[0\right]}\left(r_1+\frac{\mathsf{\Delta }r}{2}\right)W_{\left[0\right]}\left(r_1+\frac{\mathsf{\Delta }r}{2}\right)\right.\hfill \\ & -& \left.{ϵ}_{\left[2\right],\left[0\right]}\left(r_1-\frac{\mathsf{\Delta }r}{2}\right)W_{\left[0\right]}\left(r_1-\frac{\mathsf{\Delta }r}{2}\right)\right],\hfill \end{array}$$

(35)

and:

$$\frac{\partial W_{\left[0\right]}(r_1-2^{-1}\mathsf{\Delta }r)}{\partial t}=-\frac{\partial W_{\left[0\right]}(r_1+2^{-1}\mathsf{\Delta }r)}{\partial t}\frac{\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right)}{\mathsf{\Omega }\left(r_1\right)}.$$

(36)

For illustrative purposes, let us compare the above expressions with standard chemical kinetics [11]. Let us suppose that particles 1 and 2 are similar (or identical, ignoring possible effects due to quantum-mechanical symmetrization). For a chemical reaction involving two chemical species $A_1$, $A_1$ and $A_{\mathrm{d}}$: $A_{\mathrm{d}}\to A_1+A_1$, and its inverse reaction $A_1+A_2\to A_{\mathrm{d}}$, the phenomenological kinetic equations read as [11]:

$$\frac{d{C}_{\mathrm{d}}}{dt}=k^{+}{C}_1^2-k^{-}{C}_{\mathrm{d}},$$

(37)

$$\frac{d{C}_1}{dt}=-k^{+}{C}_1^2+k^{-}{C}_{\mathrm{d}},$$

(38)

where $C_1$ and $C_{\mathrm{d}}$ are the t-dependent concentrations of $A_1$ and $A_{\mathrm{d}}$, with $k^{+}$ and $k^{-}$ being suitable forward and backward rate constants. We stress that (37)–(38) are non-linearly-coupled equations for $C_1$ and $C_{\mathrm{d}}$. In principle, the $W_{\left[0\right]}$’s in different regions arising in the approximate linear Equations (35) and (36) (with the CM evolution factored out) do not appear to coincide with the concentrations in (37) and (38), so that, strictly speaking, there should be no essential contradiction or problem of principle in the simultaneous approximate validity of both pairs of equations. In spite of that, one could still ask whether there could be some analogy pointing out towards the approximate compatibility between the linear (35) and (36) and the non-linear (37) and (38). We shall outline an argument indicating that such a consistency appears to hold, at least in a semi-quantitative sense and adequately close to equilibrium. Recall that the $W_{\left[0\right]}$’s have the same dimension as ${\hslash }^{-3}$. We shall introduce the t-dependent quantities $a_{\mathrm{d}}$ and $a_{\mathrm{CS}}$ through $W_{\left[0\right]}(r_1-2^{-1}\mathsf{\Delta }r))=k^{\prime }{a}_{\mathrm{d}}$ and $W_{\left[0\right]}(r_1+2^{-1}\mathsf{\Delta }r))=k^{\prime }{k}^{\prime \prime }{a}_{\mathrm{CS}}^2$. Here, $k^{\prime }$ is an arbitrary constant, to be fixed from the outset, while $k^{\prime \prime }$ is another constant to be determined below. The underlying physical purpose is to interpret $a_{\mathrm{d}}$ and $a_{\mathrm{CS}}$ as proportional to concentrations inside $\mathsf{\Omega }\left(r_1\right)$ and $\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right)$, respectively. For large t, we argue that $a_{\mathrm{CS}}$ equals the constant equilibrium value $a_{\mathrm{CS},\mathrm{eq}}$ plus a smaller (decreasing) t-dependent contribution, so that: $a_{\mathrm{CS}}(d{a}_{\mathrm{CS}}/dt)\simeq a_{\mathrm{CS},\mathrm{eq}}(d{a}_{\mathrm{CS}}/dt)$. Then, Equations (35) and (36) become:

$$\frac{\partial a_{\mathrm{d}}}{\partial t}\simeq -\frac{q_{\mathrm{eq}}^{2}D}{m^2\mathsf{\Delta }r\mathsf{\Omega }\left(r_1\right)}{\int }_{S\left(r_1\right)}dS\left(r_1\right)\left[{ϵ}_{\left[2\right],\left[0\right]}\left(r_1+\frac{\mathsf{\Delta }r}{2}\right)k^{\prime \prime }{a}_{\mathrm{CS}}^2-{ϵ}_{\left[2\right],\left[0\right]}\left(r_1-\frac{\mathsf{\Delta }r}{2}\right)a_{\mathrm{d}}\right],$$

(39)

$$\begin{array}{ccc}\hfill \frac{\partial a_{\mathrm{CS}}}{\partial t}& \simeq & \frac{q_{\mathrm{eq}}^{2}D}{m^2\mathsf{\Delta }r(2k^{\prime \prime }{a}_{CS,eq})(\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right))}{\int }_{S\left(r_1\right)}dS\left(r_1\right)\left[{ϵ}_{\left[2\right],\left[0\right]}\left(r_1+\frac{\mathsf{\Delta }r}{2}\right)k^{\prime \prime }{a}_{\mathrm{CS}}^2\right.\hfill \\ & -& \left.{ϵ}_{\left[2\right],\left[0\right]}\left(r_1-\frac{\mathsf{\Delta }r}{2}\right)a_{\mathrm{d}}\right].\hfill \end{array}$$

(40)

Notice that $k^{\prime }$ has factored out. Next, we compare the phenomenological Equations (37)–(38) to Equations (39) and (40), interpreting $a_{\mathrm{d}}$ as $C_{\mathrm{d}}$ and $a_{\mathrm{CS}}$ as $C_1$. The two sets are compatible with each other if $k^{\prime \prime }{a}_{\mathrm{CS},\mathrm{eq}}$ fulfills: $\mathsf{\Omega }\left(r_1\right)=(\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right))2k^{\prime \prime }{a}_{\mathrm{CS},\mathrm{eq}}$. The latter equation and the restriction of $W_{\left[0\right]}(r_1+2^{-1}\mathsf{\Delta }r))=k^{\prime }{k}^{\prime \prime }{a}_{\mathrm{CS}}^2$ at equilibrium determine $k^{\prime \prime }$ and $a_{\mathrm{CS},\mathrm{eq}}$ (once $k^{\prime }$ has been fixed). Then, for large t, the non-linear Equations (37) and (38) and the linear Equations (39) and (40) would be compatible.

5. Thermal and Chemical Equilibria

Let us now return to Equations (15) and (20). Upon integration over the large finite volume Ω:

$$\int d^3\mathit{x}\int d^3\mathit{q}{W}_{\mathrm{re},\mathrm{eq}}\left(\mathit{x},\mathit{q}\right)=\frac{1}{q_{\mathrm{eq}}^3}\sum _{j}exp[-\beta E_j]\int d^3\mathit{x}{\phi }_j\left(\mathit{x}\right){\phi }_j^{*}\left(\mathit{x}\right).$$

(41)

Then, the right-hand side of Equation (41) contains the contribution of the discrete spectrum plus the one due to the whole, almost continuous spectrum, which is finite (as long as Ω is finite, but diverging like Ω if ${\mathsf{\Omega }}^{1/3}\to \infty$). If, for a finite (adequately large) Ω, T decreases within the range specified in Section 3.2, then the discrete contribution to $W_{\mathrm{re},\mathrm{eq}}$ dominates: the two particles remain bound to each other. As T increases (for fixed large Ω), the almost continuous spectrum contribution to $W_{\mathrm{re},\mathrm{eq}}$ dominates. Those two extreme situations and the intermediate ones could correspond to thermal equilibrium, but there may be, in general, an imbalance between the discrete and the CS contributions to (41), unless some further condition be fulfilled. To grasp the latter, notice that there exist key ranges of intermediate T and ${\mathsf{\Omega }}^{1/3}$, such that the orders of magnitude of the discrete and the CS contributions are not significantly different from each other. From Equation (20), the discrete contribution to (41) is:

$$\frac{1}{q_{\mathrm{eq}}^3}exp\left(-\beta E_{\mathrm{d}}\right)\int d^3\mathit{x}{\phi }_{\mathrm{d}}\left(\mathit{x}\right){\phi }_{\mathrm{d}}^{*}\left(x\right)=\frac{1}{q_{\mathrm{eq}}^3}exp\left(-\beta E_{\mathrm{d}}\right).$$

(42)

The order of magnitude of the CS contribution in Equation (41) can be estimated by disregarding the bound state term and approximating $W_{\mathrm{re},\mathrm{eq}}\simeq W_{\mathrm{eq},CS,\left[0\right]}$ (recall (20)) through (11) and (12) (disregarding $c_1$ and $c_2$). Therefore, with $\mathsf{\Omega }-\mathsf{\Omega }\left(r_1\right)\simeq \mathsf{\Omega }$, a zeroth-order approximation to the CS contribution in (41) is:

$$\frac{1}{q_{\mathrm{eq}}^3}\frac{\mathsf{\Omega }}{8{\pi }^{3/2}{\lambda }_{\mathrm{th}}^3}.$$

(43)

There would be a comparative similarity of the orders of magnitude of (42) and (43) provided that $exp\left(-\beta E_{\mathrm{d}}\right)$ be about $\mathsf{\Omega }/\left(8{\pi }^{3/2}{\lambda }_{\mathrm{th}}^3\right)$, which holds for suitable Ω, n and T. We shall now discuss the size of the microscopically large Ω. One could argue physically that our two-particle system is contained in some macroscopic gaseous system at T, containing a very large number N of similarly interacting pairs. The gaseous system occupies a macroscopically large volume ${\mathsf{\Omega }}^{\prime }\gg \mathsf{\Omega }$. Then, one could choose $\mathsf{\Omega }\simeq {\mathsf{\Omega }}^{\prime }/N$. For typical diatomic gases at adequately low pressures and T’s not smaller than room temperature, one has that Ω is about ${10}^4$–${10}^5$ ų per pair. Such an order of magnitude for Ω could allow for the compatibility between (42) and (43).

Based on the above comparison, we now turn to the possibility of another kind of equilibrium for large t, namely chemical equilibrium. In fact, at equilibrium, the vanishing of the right-hand sides of the approximate Equations (33) and (30) can be recast as:

$$\frac{W_{\left[0\right]}(r_1+2^{-1}\mathsf{\Delta }r)}{W_{\left[0\right]}(r_1-2^{-1}\mathsf{\Delta }r)}\simeq \frac{{ϵ}_{\left[2\right],\left[0\right]}(r_1-2^{-1}\mathsf{\Delta }r)}{{ϵ}_{\left[2\right],\left[0\right]}(r_1+2^{-1}\mathsf{\Delta }r)}.$$

(44)

At equilibrium, we approximate $W_{\left[0\right]}(r_1-2^{-1}\mathsf{\Delta }r)$ by $q_{\mathrm{eq}}^{-3}exp\left(-\beta E_{\mathrm{d}}\right)\mathsf{\Omega }{\left(r_1\right)}^{-1}$ and $W_{\left[0\right]}(r_1+2^{-1}\mathsf{\Delta }r)$ by $q_{\mathrm{eq}}^{-3}{\left(8{\pi }^{3/2}{\lambda }_{\mathrm{th}}^3\right)}^{-1}$, consistently with (42) and (43), and we replace the ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$’s by their zeroth order approximations in Section 3.2. Then, (44) becomes:

$$\frac{\mathsf{\Omega }\left(r_1\right)exp\left(\beta E_{\mathrm{d}}\right)}{8{\pi }^{3/2}{\lambda }_{\mathrm{th}}^3}\simeq \frac{2{\left({\lambda }_{\mathrm{th}}\right)}^2}{{\left(\delta x\right)}^2},$$

(45)

where $\mathsf{\Omega }\left(r_1\right)$ is of order $a^3$ (see Section 3.2). Equation (44) is, within our rough approximations (say, about one order of magnitude more or less), compatible with the approximate similarity of the orders of magnitude of Equations (42) and (43) when the latter occurs (which, in turn, holds if $\mathsf{\Omega }\left(r_1\right)/\mathsf{\Omega }$ is about $2{\left({\lambda }_{\mathrm{th}}\right)}^2/{\left(\delta x\right)}^2)$). Hence, there is a simultaneous coexistence of thermal and chemical equilibria in the actual simple framework. Moreover, it appears adequate to regard Equation (44) (in which no chemical potentials have been introduced) as a counterpart of the law of mass action for equilibrium chemical reactions. Notice also that $W_{\left[0\right]}(r_1+2^{-1}\mathsf{\Delta }r)$ is then smaller than $W_{\left[0\right]}(r_1-2^{-1}\mathsf{\Delta }r)$ by a factor of about $\mathsf{\Omega }\left(r_1\right)/\mathsf{\Omega }$: as a side remark, this justifies the neglect of the CS contribution compared to the bound state one, when estimating ${ϵ}_{\left[2_{\alpha }\right],0}$ for $r_0<\left|\mathit{x}\right|<r_0+a$ in Appendix B.

6. Transitions among a Bound and Continuous Spectrum: Mean First Passage Time

In a similar fashion as in Section 4, let us replace approximately the quantum operator $D[\left[1_{\alpha }\right];ϵ]$ by a positive and both α- and x-independent constant D, and suppose that ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$ depends solely on the radial coordinate r and that it is negative. Let us restrict to $W_{\left[0\right]}=W_{\left[0\right]}(r,t)$, independent by assumption on the angular coordinates. Implicitly, we are accepting that $W_{\left[0\right]}(r,t)$ is an adequate (not to say, dominant) approximation to the full $W_{\left[0\right]}$. Then, Equation (28) becomes:

$$\frac{\partial W_{\left[0\right]}}{\partial t}=\frac{q_{\mathrm{eq}}D}{m}\left(\frac{\partial }{\partial r}+\frac{2}{r}\right)\left(-\frac{q_{\mathrm{eq}}}{m}\frac{\partial }{\partial r}\left({ϵ}_{\left[2\right],\left[0\right]}{W}_{\left[0\right]}\right)+\frac{1}{q_{\mathrm{eq}}}\frac{\partial V}{\partial r}{W}_{\left[0\right]}\right).$$

(46)

It will be useful to replace $W_{\left[0\right]}(r,t)$ by a new distribution $f=f(r,t)$ defined through $W_{\left[0\right]}(r,t)=r^{-2}f(r,t)$. Then, Equation (46) reduces to:

$$\frac{\partial f}{\partial t}=\frac{q_{\mathrm{eq}}D}{m}\frac{\partial }{\partial r}\left[-\frac{q_{\mathrm{eq}}}{m}\frac{\partial }{\partial r}\left({ϵ}_{\left[2\right],\left[0\right]}f\right)+\frac{1}{q_{\mathrm{eq}}}\frac{\partial V}{\partial r}f+\frac{2}{r}\frac{q_{\mathrm{eq}}}{m}{ϵ}_{\left[2\right],\left[0\right]}f\right].$$

(47)

We are interested in finding approximately the time required for the relative particle to proceed from the domain in which the bound state is concentrated to a continuous spectrum configuration (that is, large r). For that purpose, we shall apply the mean first passage time formalism (MFPT) [16,17,35]. We shall extend the treatment in [34]. An MFPT $\tau \left(r\right)$ (which can be regarded as another effective evolution time) is the solution of the so-called adjoint equation associated with Equation (47), namely of:

$$\frac{q_{\mathrm{eq}}^{2}D(-{ϵ}_{\left[2\right],\left[0\right]})}{m^2}\frac{{\partial }^2\tau \left(r\right)}{\partial r^2}-\frac{q_{\mathrm{eq}}D}{m}\left[\frac{1}{q_{\mathrm{eq}}}\frac{\partial V}{\partial r}+\frac{2q_{\mathrm{eq}}{ϵ}_{\left[2\right],\left[0\right]}}{mr}\right]\frac{\partial \tau \left(r\right)}{\partial r}=-1.$$

(48)

Given that $r_0$ is an adequately small radius, such that $V>0$ for $0\le r<r_0$ (i.e., the hard core), and that for the appreciably larger interval $r_0<r<a+r_0$, $V<0$ and $\left|V\right|$ takes on its largest values, let $r_2$ be considerably larger than $r_0+a$ (and so, than the $r_1$ in Section 4 and Section 5), so that V is negligible at $r_2$. Since $r_2$ is smaller than ${[\left(3\mathsf{\Omega }\right)/\left(4\pi \right)]}^{1/3}$, we shall look for the solution $\tau \left(r\right)$ of Equation (48), such that $\partial \tau \left(r\right)/\partial r=0$ at $r=r_0$ (“reflecting boundary condition”) and $\tau \left(r\right)=0$ at $r=r_2$ (“absorbing boundary condition”). One finds readily for $r_0<r<r_2$ that:

$$\begin{array}{ccc}\hfill \tau \left(r\right)& =& {\int }_r^{r_2}\frac{d{r}^{\prime }}{{(q_{\mathrm{eq}}/m)}^{2}D{r}^{\prime 2}}J\left(r^{\prime }\right),\hfill \end{array}$$

(49)

$$\begin{array}{ccc}\hfill J\left(r^{\prime }\right)& =& -{\int }_{r_0}^{r^{\prime }}d{r}^{\prime \prime }\frac{r^{\prime \prime 2}}{{ϵ}_{\left[2\right],\left[0\right]}\left(r^{\prime \prime }\right)}exp\left[-\frac{\beta }{2}{\int }_{r^{\prime \prime }}^{r^{\prime }}d{r}^{\prime \prime \prime }\frac{(\partial V/\partial r^{\prime \prime \prime })}{{ϵ}_{\left[2\right],\left[0\right]}\left(r^{\prime \prime \prime }\right)}\right].\hfill \end{array}$$

(50)

In Section 7 and Appendix F, we shall be able to approximate the operator $D[\left[1_{\alpha }\right];ϵ]$ by a positive and α-independent function $D\left(r\right)$, instead of by a constant D directly. One would like to be able to control the further approximation of the function $D\left(r\right)$ by a suitable constant D upon evaluating the MFPT, so as to perform estimates with Equation (49). For that purpose, the following remark regarding Equation (49) will be useful. Let us come back to Equation (28), replace in it the operator $D[\left[1_{\alpha }\right];ϵ]$ by a function $D\left(r\right)$ and rewrite the resulting Equation (28) as $\partial W_{\left[0\right]}/\partial t=L_S{W}_{\left[0\right]}$. Here, $L_S$ is the operator defined by the right-hand side of Equation (28) when the operator $D[\left[1_{\alpha }\right];ϵ]$ is replaced by the function $D\left(r\right)$. Let $L_S^{+}$ be the adjoint of $L_S$. An MFPT ${\tau }_L\left(r\right)$ for $L_S^{+}$ can be introduced, by extending [16,17,35], as the solution of the differential equation $L_S^{+}{\tau }_L\left(r\right)=-1$ in $r_0<r<r_2$, with similar boundary conditions $\partial {\tau }_L\left(r\right)/\partial r=0$ at $r=r_0$ (“reflecting boundary condition”) and ${\tau }_L\left(r\right)=0$ at $r=r_2$ (“absorbing boundary condition”). A calculation similar to that yielding Equation (49) leads to:

$${\tau }_L\left(r\right)={\int }_r^{r_2}\frac{d{r}^{\prime }}{{(q_{\mathrm{eq}}/m)}^{2}D\left(r^{\prime }\right)r^{\prime 2}}J\left(r^{\prime }\right),$$

(51)

with the same $J\left(r^{\prime }\right)$ as in Equation (49). The interest of Equation (51) is two-fold: either $D\left(r^{\prime }\right)$ inside the integral could now be approximated by some constant D (some average of $D\left(r^{\prime }\right)$ in $r<r^{\prime }<r_2$) and, so, allow for the use of $\tau \left(r\right)$ and Equation (49) (as done below in this section) or variations of $D\left(r^{\prime }\right)$ inside the integral could be considered in certain cases (see Section 7). We rewrite Equation (50) as:

$$\begin{array}{ccc}\hfill J\left(r^{\prime }\right)& =& {\int }_{r_0}^{r_0+a}d{r}^{\prime \prime }\frac{r^{\prime \prime 2}}{(-{ϵ}_{\left[2\right],\left[0\right]}\left(r^{\prime \prime }\right))}exp\left[\frac{\beta }{2}{\int }_{r^{\prime \prime }}^{a+r_0}d{r}^{\prime \prime \prime }\frac{(\partial V/\partial r^{\prime \prime \prime })}{(-{ϵ}_{\left[2\right],\left[0\right]}\left(r^{\prime \prime \prime }\right))}+\frac{\beta }{2}{\int }_{r_0+a}^{r^{\prime }}d{r}^{\prime \prime \prime }\frac{(\partial V/\partial r^{\prime \prime \prime })}{(-{ϵ}_{\left[2\right],\left[0\right]}\left(r^{\prime \prime \prime }\right))}\right]\hfill \\ & +& {\int }_{r_0+a}^{r^{\prime }}d{r}^{\prime \prime }\frac{r^{\prime \prime 2}}{(-{ϵ}_{\left[2\right],\left[0\right]}\left(r^{\prime \prime }\right))}exp\left[\frac{\beta }{2}{\int }_{r^{\prime \prime }}^{r^{\prime }}d{r}^{\prime \prime \prime }\frac{(\partial V/\partial r^{\prime \prime \prime })}{(-{ϵ}_{\left[2\right],\left[0\right]}\left(r^{\prime \prime \prime }\right))}\right].\hfill \end{array}$$

(52)

Recall that it was supposed that $V\left(r\right)$ has no discontinuities and that it becomes negligible as $r>r_0+a$. As stated in Section 3.2, $V_0$ lies between 1 and 10 eV, and T satisfies Assumptions (a) and (b) there. Let r be larger than $r_0+a$ (with $r<r_2$). Then, it is permissible to replace ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$ by its zeroth order approximation, given in Section 3.2. Notice that ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$ for $r<r_0$ is not needed. After that approximation, one performs the three integrations over $r^{\prime \prime \prime }$ in (52). We take $V\left(r^{\prime }\right)\simeq 0$ for $r^{\prime }>r_0+a$, $V(r_0+a)\simeq -V_0$ and $V\left(r^{\prime \prime }\right)\simeq 0$ for $r^{\prime }>r^{\prime \prime }>r_0+a$. Hence, Equation (52) simplifies to:

$$J\left(r^{\prime }\right)\simeq \frac{{\left(\delta x\right)}^2}{{\left({\lambda }_{\mathrm{th}}\right)}^2}exp\left(\beta V_0\right){\int }_{r_0}^{r_0+a}d{r}^{\prime \prime }{r}^{\prime \prime 2}exp\left[\frac{m{\left(\delta x\right)}^2}{{\hslash }^2}\left(V(r_0+a)-V\left(r^{\prime \prime }\right)\right)\right]+2{\int }_{r_0+a}^{r^{\prime }}d{r}^{\prime \prime }{r}^{\prime \prime 2}.$$

(53)

Then, approximating $D\left(r^{\prime }\right)$ by some (average) constant D and using $\tau \left(r\right)$ together with Equation (49), we get:

$$\begin{array}{ccc}\hfill \tau \left(r\right)& \simeq & \frac{m^2}{q_{\mathrm{eq}}^{2}D}\left[\frac{{\left(\delta x\right)}^2}{{\left({\lambda }_{\mathrm{th}}\right)}^2}exp\left(\beta V_0\right)\left(\frac{1}{r}-\frac{1}{r_2}\right)J_1+\frac{2}{3}\left(\frac{r_2^2-r^2}{2}-{(r_0+a)}^3\left(\frac{1}{r}-\frac{1}{r_2}\right)\right)\right]\hfill \\ & \simeq & \frac{m^2}{q_{\mathrm{eq}}^{2}D}\frac{{\left(\delta x\right)}^2}{{\left({\lambda }_{\mathrm{th}}\right)}^2}exp\left(\beta V_0\right)\left(\frac{1}{r}-\frac{1}{r_2}\right)J_1,\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill J_1& =& {\int }_{r_0}^{r_0+a}d{r}^{\prime \prime }{r}^{\prime \prime 2}exp\left[\frac{m{\left(\delta x\right)}^2}{{\hslash }^2}\left(V(r_0+a)-V\left(r^{\prime \prime }\right)\right)\right].\hfill \end{array}$$

(55)

The last contribution in the rightmost side of (54) has been neglected by invoking that $exp\left(\beta V_0\right)$ is adequately larger than unity (unless $k_{\mathrm{B}}T$ is about $V_0$) and the approximate relationships in Section 5. As a rough approximation, we replace $exp\left[\frac{m{\left(\delta x\right)}^2}{{\hslash }^2}\left(V(r_0+a)-V\left(r^{\prime \prime }\right)\right)\right]$ for any $r_0<r^{\prime \prime }<r_0+a$ inside the integral in (the β- independent) $J_1$ by some ($r^{\prime \prime }$ and T independent) κ. The magnitudes of $J_1$, the constant κ and related quantities will be discussed briefly in Section 8. Then, Equation (54) reduces to:

$$\tau \left(r\right)\simeq \frac{exp\left(\beta V_0\right)m^2\kappa [{(r_0+a)}^3-r_0^3]{\left(\delta x\right)}^2}{3D{\hslash }^2}\left(\frac{1}{r}-\frac{1}{r_2}\right)\simeq \frac{exp\left(\beta V_0\right)m^2\kappa a^3{\left(\delta x\right)}^2}{3D{\hslash }^{2}r}.$$

(56)

Consistently with the assumptions in Section 2.3 and to fix the ideas, the shape of V resembles that of the Morse potential. Moreover, $\tau {\left(r\right)}^{-1}$ could be interpreted as some sort of rate constant for the transition between the bound state and the continuous spectrum states. At this point, we recall the exponential-like features of the rate constant in the classical Arrhenius formula and in the transition state theory [35,46]. We see that our quantum-mechanical approach has led approximately to the exponential Arrhenius factor, namely $exp\left(\beta V_0\right)$. Notice that Equation (51), by extending the approximations for $\tau \left(r\right)$, also leads to the exponential Arrhenius factor for ${\tau }_L\left(r\right)$. However, our approach fails so far to predict $\tau \left(r\right)$ (or ${\tau }_L\left(r\right)$) strictly, as long as it involves an unknown prefactor, namely either the function $D\left(r\right)$ or the average constant D (approximating roughly the unknown operator $D[\left[1_{\alpha }\right];ϵ]$). Through further approximations, estimates of D and of $\tau \left(r\right)$ will be provided in Section 7.

7. Estimating Effective Evolution Times

We disregard the fact that $D[\left[1_{\alpha }\right];ϵ]$ is an operator and approximate it by the function $D\left(r\right)$. The latter can be estimated by extending directly to three dimensions the one-dimensional analysis in Appendix F. This yields the following rough estimate for $D\left(r\right)$ (in seconds):

$$D\left(r\right)\simeq \frac{1}{{(-M_{\left[0\right],\left[1_{\alpha }\right]}{M}_{\left[1_{\alpha }\right],\left[0\right]})}^{1/2}}.$$

(57)

In turn, by recalling Section 3.2, we shall interpret $-M_{\left[0\right],\left[1_{\alpha }\right]}{M}_{\left[1_{\alpha }\right],\left[0\right]}$ as the product of the constant $\hslash /\left({\lambda }_{\mathrm{th}}m\delta x\right)$ (associated with $M_{\left[0\right],[\left(1_{\alpha }\right]}$, that is with $(q_{\mathrm{eq}}/m)(\partial /\partial x_{\alpha })$) times the overall order of magnitude estimate for $M_{\left[1_{\alpha }\right],\left[0\right]}$ (namely, the sum of the estimates (i) plus (ii) in Section 3.2, for $r_0<\left|\mathit{x}\right|<r_0+a$ and $\left|\mathit{x}\right|>r_0+a$, respectively), which depends on $r=\left|x\right|$. Then, Equation (57) does depend on r. For $r_0<r<r_0+a$, the dominant contributions give ($\delta V$ in eV, ${\lambda }_{\mathrm{th}}$, n and $\delta x$ being discussed in Section 3.2):

$$D\left(r\right)\simeq {\left(\frac{m{\left(\delta x\right)}^2}{\delta V}\right)}^{1/2}\simeq \frac{n^{1/2}}{3{\left(10\right)}^{1/2}{(\delta V/1\mathrm{eV})}^{1/2}}\left(\frac{\delta x}{1\mathsf{\AA }}\right)\times {10}^{-13}\mathrm{s},$$

(58)

which is T-independent.

For $r>r_0+a$, the overall order of magnitude estimate for $M_{\left[1_{\alpha }\right],\left[0\right]}$ is about $(1/2)\hslash /\left({\lambda }_{\mathrm{th}}m\delta x\right)$. Then:

$$D\left(r\right)\simeq \frac{{\lambda }_{\mathrm{th}}m\delta x}{2^{1/2}\hslash }\simeq \frac{n}{92}\left(\frac{\delta x}{1\mathsf{\AA }}\right)\left(\frac{{\lambda }_{\mathrm{th}}}{1\mathsf{\AA }}\right)\times {10}^{-11}\mathrm{s},$$

(59)

which depends on T like $T^{-1/2}$.

Therefore, for adequately small ${\lambda }_{\mathrm{th}}$ (fulfilling Assumptions (a) and (b) in Section 3.2, the estimate in (59) would decrease and become closer to that in (58). Then, provided that the orders of magnitude of the latter two are not too different from each other, the constant D considered in Section 6 and, in particular, in Equations (49), (54) and (56) is some average of (58) and (59). Recall that the average time for ironing out local thermal inhomogeneities over one nanometer in a gas under normal conditions is about ${10}^{-11}$ s [3].

Then, returning to Section 3.3 and as in Appendix F, it is reasonable to interpret ${\tau }_{\max }^{*}\simeq D/2$ and to expect that the long-time approximation be adequate for $t>D/2$.

Turning to the estimate for $\tau \left(r\right)$ in the rightmost side of (56), we get ($r_0+a<r<10a$):

$$\tau \left(r\right)D\simeq \left(\frac{n^2}{{10}^4}\right){\left(\frac{\delta x}{1\mathsf{\AA }}\right)}^2{\left(\frac{a}{1\mathsf{\AA }}\right)}^2\left(\frac{a}{r}\right)\kappa exp\left(\beta V_0\right)\times {10}^{-21}{\mathrm{s}}^2,$$

(60)

in which D is an average of (58) and (59) when they do not differ much.

Notice that if (58) is appreciably smaller than (59), then the use of Equation (51) would be more adequate than that of Equation (49). For $r_0<r<r_0+a$, the integral in (51) could then be restricted to ${\int }_r^{r_0+a}\frac{d{r}^{\prime }}{{(q_{\mathrm{eq}}/m)}^{2}D\left(r^{\prime }\right)r^{\prime 2}}J\left(r^{\prime }\right)$, where $D\left(r^{\prime }\right)$ should be replaced by (58). For $r>r_0+a$, $D\left(r^{\prime }\right)$ should be replaced by (59) in the integral in (51). Both resulting integrals could be estimated as in Section 6, so as to yield further values of the MFPT. This section and the preceding one have implemented the counterparts of step 7 in [34].

8. Comparison with Other Approaches: Recombination and Possible Extensions

Notice that, through Condition 1 in Section 2.3, we have restricted to a potential $V\left(r\right)$, which, in particular, is attractive ($<0$) in the interval $r_0<r<{\left(3\mathsf{\Omega }/4\pi \right)}^{1/3}$ (say, with a global structure resembling that of a Morse-like one). From the conditions in Section 2.3 and through the developments and approximations up to and including Section 3, we have arrived at the quantum Equation (28), which, in particular, provided the basis for the MFPT analysis in Section 6. Such a study can be compared to that in Subsection VII.C.4 in [35] for a radial classical Smoluchowski equation with a damping or friction parameter (${\gamma }_{fr}$, given from the outset) and a similar spherically-symmetric potential (Morse-like), in connection with diffusive problems and, specifically, with a recombination chemical reaction (see, in particular, Equation (7.30) and Figure 24 in [35]). The corresponding MFPT analysis yields Equations (7.31)–(7.33) in [35]. One key difference between the classical Equations (7.30)–(7.31) in [35] and our Equations (46)–(48) is that the latter contain the quantum contribution ${ϵ}_{\left[2\right],\left[0\right]}$. Another difference is that our approach allows one to estimate, after a number of approximations, the quantum counterpart of ${\gamma }_{fr}$, namely the values of $D\left(r\right)$ (or D): see Appendix F and Section 7. Notice that the values of ${ϵ}_{\left[2\right],\left[0\right]}$ in $r_0<r<a+r_0$ (dominated by the bound state), estimated in Section 3.2, are responsible for the Arrhenius factor $exp\left(\beta V_0\right)$ in (56). Our (quantum) solution (49) and (50) is also applicable to similar diffusive problems and is the counterpart of (with the same boundary conditions as) the classical solution in Equation (7.33) in [35].

It is worth emphasizing that the MFPT solution in Equation (7.32) in [35] corresponds to the opposite boundary conditions (absorption at $r_0$ and reflection far outside the range of the potential), and it has been applied to recombination processes in chemical reactions: see [35] and the references therein for generalizations. This suggests the interest of extending succinctly our study in Section 6 above so as to find the MFPT solution ${\tau }_{\mathrm{rec}}\left(r\right)$ for Equations (47) and (48) for boundary conditions $\partial {\tau }_{\mathrm{rec}}\left(r\right)/\partial r=0$ at $r=r_2$ (reflection) and ${\tau }_{\mathrm{rec}}\left(r\right)=0$ at $r=r_0$ (absorption boundary) and to compare it with [35]. The operator $D[\left[1_{\alpha }\right];ϵ]$ is approximated by the same positive and α-independent function $D\left(r\right)$ as in Section 6. A new computation, similar to that yielding Equations (50) and (51), gives readily for $r_0<r<r_2$:

$${\tau }_{\mathrm{rec}}\left(r\right)={\int }_{r_0}^r\frac{d{r}^{\prime }}{{(q_{\mathrm{eq}}/m)}^{2}D\left(r^{\prime }\right)r^{\prime 2}}{J}_{\mathrm{rec}}\left(r^{\prime }\right),$$

(61)

$$J_{\mathrm{rec}}\left(r^{\prime }\right)=-{\int }_{r^{\prime }}^{r_2}d{r}^{\prime \prime }\frac{r^{\prime \prime 2}}{{ϵ}_{\left[2\right],\left[0\right]}\left(r^{\prime \prime }\right)}exp\left[\frac{\beta }{2}{\int }_{r^{\prime }}^{r^{\prime \prime }}d{r}^{\prime \prime \prime }\frac{(\partial V/\partial r^{\prime \prime \prime })}{{ϵ}_{\left[2\right],\left[0\right]}\left(r^{\prime \prime \prime }\right)}\right].$$

(62)

which are the actual quantum counterparts of Equations (7.32), (7.34) and (7.35) in [35]. Here, Equations (61) and (62) would give approximately the (recombination) time required for the relative particle to proceed from the domain corresponding to the continuous spectrum configuration ($r_0+a<r$) to the one in which the bound state is concentrated (that is, $r_0<r<r_0+a$). Equations (61) and (62) can be estimated, by approximating $D\left(r^{\prime }\right)$ by a constant (D) and extending readily the approximations leading from (49) and (50) to (54). One finds:

$${\tau }_{\mathrm{rec}}\left(r\right)\simeq {\tau }_{\mathrm{rec},1}\left(r\right)+{\tau }_{\mathrm{rec},2}\left(r\right),$$

(63)

where:

$${\tau }_{\mathrm{rec},1}\left(r\right)=\frac{m\beta }{4D}\frac{exp\left(-\beta V_0\right)}{2}\frac{\left(r_3^2-{(r_0+a)}^3\right)}{3}{\int }_{r_0}^r\frac{d{r}^{\prime }}{r^{\prime 2}}exp\left[\frac{m{\left(\delta x\right)}^2}{{\hslash }^2}\left(V\left(r^{\prime }\right)-V(r_0+a)\right)\right],$$

(64)

$${\tau }_{\mathrm{rec},2}\left(r\right)=\frac{m^2{\left(\delta x\right)}^2}{{\hslash }^{2}D}{\int }_{r_0}^r\frac{d{r}^{\prime }}{r^{\prime 2}}{\int }_{r^{\prime }}^{r_0+a}d{r}^{\prime \prime }{r}^{\prime \prime 2}exp\left[\frac{m{\left(\delta x\right)}^2}{{\hslash }^2}\left(V\left(r^{\prime }\right)-V\left(r^{\prime \prime }\right)\right)\right],$$

(65)

with $exp\left(-\beta V_0\right)$ being an Arrhenius-like factor with a negative exponent, decreasing as T decreases. Notice that $D{\tau }_{\mathrm{rec},1}\left(r\right)$ depends on temperature only through $\beta exp\left(-\beta V_0\right)$, while $D{\tau }_{\mathrm{rec},2}\left(r\right)$ is temperature-independent. The possible magnitudes of the exponentials $exp\left[\frac{m{\left(\delta x\right)}^2}{{\hslash }^2}\left(V\left(r^{\prime }\right)-V\left(r^{\prime \prime }\right)\right)\right]$ for various $r^{\prime }$ and $r^{\prime \prime }$, which appear in $J_1$ in Equation (55) (and so, the related constant κ in Equation (56)), in ${\tau }_{\mathrm{rec},1}\left(r\right)$ and in ${\tau }_{\mathrm{rec},2}\left(r\right)$ depend rather strongly on whether V varies rapidly or smoothly and on the sign of $V\left(r^{\prime }\right)-V\left(r^{\prime \prime }\right)$. Let $V\left(r^{\prime }\right)-V\left(r^{\prime \prime }\right)>0$ (as, otherwise, those exponentials could be neglected). For suitably smooth V (say, rather small $\delta V$), those exponentials would not be much greater than unity (for instance, κ could not exceed ${10}^2$ or ${10}^3$), but they could be certainly large for rapidly-varying V. In any case, regardless of the magnitudes of those exponentials, we stress that $J_1$ (and hence, κ), together with the integrals in ${\tau }_{\mathrm{rec},1}\left(r\right)$ and ${\tau }_{\mathrm{rec},2}\left(r\right)$ are temperature independent. A numerical study of those integrals lies outside our scope here. The quantum Equations (63)–(65) can be compared with the classical Equations (74) and (75) in [35], which appear to yield a different temperature dependence. Recall that our methods rely on Assumption (a) in Section 3.2, imposing a lower limit on temperature.

The activated barrier crossing problem (say, for transitions from a bound state to continuous state configurations) has attracted enormous research attention in connection with chemical reactions [17,35,46]: in one simple version, it corresponds to a potential $V\left(r\right)$, which presents a “hard core” for $0\le r<r_0$; it is attractive in the interval $r_0<r<r_3$, repulsive in the interval $r_3<r<r_4$ and vanishes fast as $r\to +\infty$. Our Equation (28) and the arguments up to and including Section 3 would also hold if Condition 1 in Section 2.3 is replaced by: 1*. $V\left(r\right)$ is repulsive ($>0$) for $0\le r<r_0$ (“hard core”, with adequately small $r_0$), attractive ($<0$) in the interval $r_0<r<r_3$, repulsive ($>0$) in the interval $r_3<r<{\left(3\mathsf{\Omega }/4\pi \right)}^{1/3}$ and vanishes fast as $r\to {\left(3\mathsf{\Omega }/4\pi \right)}^{1/3}$. The estimate of ${ϵ}_{\left[2\right],\left[0\right]}$ in $r_0<r<a+r_0$ (dominated by the bound state) in Section 3.2 should now be revised, but it would not appear to spoil the validity of the corresponding Equation (28). The MFPT analysis in Section 6 up to and including Equations (46)–(51) would also apply (as would that leading to Equations (61) and (62)). However, the analysis in Section 4 and Section 5 and the approximate estimates in Section 6 beyond (51) should be revised. Therefore, in order to apply our methods to the activated barrier crossing problem, those parts of our analysis should have to be extended, which lie outside our scope here.

9. Conclusions and Discussions

In this work, we have put forward a nontrivial generalization of the methods previously used in the study of non-equilibrium DNA thermal denaturation, presented in [34], to provide a simple model for a binary chemical reaction starting from a quantum mechanical framework. The close correspondence between the methods and successive steps in [34] and the present work, as well as the main differences between them, have been highlighted. We have treated a quantum two-particle system in three spatial dimensions, subject to an attractive potential V and to a “heat bath” (HB) at thermal equilibrium at absolute temperature $T>0$. We have focused on a quantum regime typical of chemical reactions, such that: (a) the thermal wavelength ${\lambda }_{\mathrm{th}}$ is shorter than the range of the attractive potential; (b) the energy for the discrete bound state $|E_{\mathrm{d}}|$ (about the magnitude of the attractive potential within the negative region) is, in magnitude, $\ge (3/2)k_{\mathrm{B}}T$. Our starting point has been the non-equilibrium time-reversible Wigner equation without explicit dissipation (as justified in Appendix G). We have separated the overall center of mass (CM) of the relative motion. We have focused on the dynamics of the non-equilibrium Wigner function $W_{\mathrm{re}}$ for the relative particle motion. We have considered an infinite family of orthogonal polynomials $H_{\left[n\right]}$ (depending on the relative momentum), generated by the relative equilibrium Wigner function $W_{\mathrm{re},\mathrm{eq}}$ (non-Gaussian and non-positive, in general). These orthogonal polynomials have generated (by integrating over the relative momentum) non-equilibrium moments $W_{\left[n\right]}$ of $W_{\mathrm{re}}$. In turn, the general non-equilibrium Wigner equation for $W_{\mathrm{re}}$ has implied a linear nonequilibrium hierarchy for the $W_{\left[n\right]}$’s (time-reversible, as well, and depending on the relative position). Under the above Assumptions (a) and (b), and for sufficiently long times $t\ge {\tau }_{\max }^{*}$, the linear nonequilibrium hierarchy has been approximated by the irreversible Smoluchowski-like Equation (28) for the lowest moment $W_{\left[0\right]}$. Our Equation (28) contains quantum effects; namely the operator $D[\left[1_{\alpha }\right];ϵ]$ and, in particular, the function ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$. At a later stage, the operator $D[\left[1_{\alpha }\right];ϵ]$ has been approximated by the constant D, so that ${\tau }_{\max }^{*}\simeq D/2$. This has led to effective linear kinetic equations allowing, in turn, for approximate comparisons (with standard phenomenological non-linear kinetic equations for large t) and between thermal versus chemical equilibria. Then, we proceed to the simpler Equation (46), which depends only on t and on the radial separation r between both particles and still preserves a quantum signature via the term ${ϵ}_{\left[2_{\alpha }\right],\left[0\right]}$. The order of magnitude of D has been estimated. Equation (46) has allowed us to apply the mean first passage time (MFPT) formalism [16,17,35] which, in turn, has given rise to an explicit representation for the time $\tau =\tau \left(r\right)$ required for the two particles to proceed from the bound state to the continuous spectrum configurations, as a function of r. An explicit formula for $\tau \left(r\right)$, in terms of D and displaying the genuine Arrhenius exponential factor $exp\left[V_0/\left(k_{\mathrm{B}}T\right)\right]$, has been obtained and discussed briefly. A consistent extension, in which the operator $D[\left[1_{\alpha }\right];ϵ]$ is replaced by a suitable function of the radial distance, has also been outlined. The MFPT has also been applied briefly to chemical recombination. Our MFPT studies are compared also briefly to other MFPT studies by other authors [35]. For the sake of a simpler presentation without loss of generality, the various approximations leading to the Smoluchowski-like equation for $W_{\left[0\right]}$ have been carried out in detail for a one-dimensional case in Appendix C, Appendix D, Appendix E and Appendix F. We point out that no use of non-equilibrium Wigner functions, moments and hierarchies was made in the works [22,23,24,25,26,27,28,29]. Among several open problems (and hence, possible future uses of the methods in the present paper), we quote: (a2) extensions to include more than one bound state between the two particles, thereby allowing for transitions among those states (eventually, for instance, transitions resembling vibrational ones approximately) and (b2) generalizations to more general two-particle chemical reactions (like $A+B\to C+D$).