Non-Equilibrium Liouville and Wigner Equations: Classical Statistical Mechanics and Chemical Reactions for Long Times

Abstract

We review and improve previous work on non-equilibrium classical and quantum statistical systems, subject to potentials, without ab initio dissipation. We treat classical closed three-dimensional many-particle interacting systems without any “heat bath” ($hb$), evolving through the Liouville equation for the non-equilibrium classical distribution $W_c$, with initial states describing thermal equilibrium at large distances but non-equilibrium at finite distances. We use Boltzmann’s Gaussian classical equilibrium distribution $W_{c,eq}$, as weight function to generate orthogonal polynomials ($H_n$’s) in momenta. The moments of $W_c$, implied by the $H_n$’s, fulfill a non-equilibrium hierarchy. Under long-term approximations, the lowest moment dominates the evolution towards thermal equilibrium. A non-increasing Liapunov function characterizes the long-term evolution towards equilibrium. Non-equilibrium chemical reactions involving two and three particles in a $hb$ are studied classically and quantum-mechanically (by using Wigner functions W). Difficulties related to the non-positivity of W are bypassed. Equilibrium Wigner functions $W_{eq}$ generate orthogonal polynomials, which yield non-equilibrium moments of W and hierarchies. In regimes typical of chemical reactions (short thermal wavelength and long times), non-equilibrium hierarchies yield approximate Smoluchowski-like equations displaying dissipation and quantum effects. The study of three-particle chemical reactions is new.

1. Introduction

Research leading from equilibrium statistical mechanics [1,2,3,4,5,6] to non-equilibrium statistical mechanics [2,7,8,9] in both classical and quantum regimes faces many open fundamental difficulties.

One crucial hallmark is that the off-equilibrium evolution of statistical systems displays stochasticity: see [10,11,12] in the classical regime and [13,14,15,16,17,18,19] in the quantum one, specifically, in the framework of the theory of open quantum systems. For complementary background in classical (thermodynamical) frameworks, see [20,21,22].

Chemical reactions constitute a very important and wide class of phenomena in which equilibrium and non-equilibrium statistical mechanics play essential roles: such phenomena are genuinely quantum-mechanical, although classical approaches can provide useful information under various approximations. See [5,6] and references therein for studies devoted to the equilibrium chemical constant. Much research has also been devoted to non-equilibrium chemical reactions: see in particular [11,20,22], for stochastic and thermodynamical approaches. The non-equilibrium rate constant is treated in [11] by means of a thorough analysis of the Kramers equation.

Use has been made, by previous authors, of infinite hierarchies in the analysis of several stochastic equations [23,24]. The present author has carried out several detailed analyses of the classical Liouville and quantum Wigner equation, through infinite hierarchies for suitable non-equilibrium moments, in order to analyze approaches to thermal equilibrium in long-term approximations [25,26,27,28,29,30,31,32]. Motivated by previous work on non-equilibrium DNA thermal denaturation (the dynamics of two interacting classical three-dimensional macromolecular chains) [33], binary chemical reactions have been studied quantum-mechanically, based on the non-equilibrium Wigner function [34].

The intention of this work is to review and extend the work mentioned above in [25,26,27,28,29,30,31,32] and in [34]. To outline, we shall focus on the dynamics of non-equilibrium non-relativistic systems implied by the classical Liouville equation and the quantum Wigner one [2]. Throughout our study, we shall exclude by assumption, from the outset and in all cases, ab initio dissipation; there will be vanishing friction due to any other system. We shall omit throughout any detailed analysis of the states of the “heat bath” $hb$ (just retaining the fact that it establishes the temperature) and, in an effective way, focus exclusively on the states of the systems under study. The corresponding equilibrium distributions will be used to generate an infinite family of orthogonal polynomials. In turn, the latter will allow construction of non-equilibrium moments which, through either the Liouville equation or the Wigner one, will imply infinite linear hierarchies. These hierarchies, under suitable assumptions and approximations, will yield an approach to thermal equilibrium for the long term. The formalism will enable a theoretical analysis of chemical reactions.

The above program faces various difficulties, if it is to be treated successively. We shall begin with simpler one-dimensional settings (namely, one particle subject to an external potential and to a $hb$ at thermal equilibrium) in both classical (Section 2) and quantum (Section 4) cases, to illustrate how difficulties will be bypassed. Section 3 deals with classical closed three-dimensional many-particle interacting systems without any $hb$, evolving through the Liouville equation for the non-equilibrium classical probability distribution $W_c$ and with initial states describing thermal equilibrium at large distances but non-equilibrium at finite ones. Section 5 outlines useful features of the non-equilibrium dynamics of $N>1$ quantum particles subject to a $hb$ using Wigner functions, which will be useful in Section 6 and Section 7. Section 6 studies non-equilibrium chemical reactions involving two particles. Section 7 deals with a new short study of non-equilibrium chemical reactions involving three particles. Section 8 summarizes the contents of this article and distinguishes review material from new material, along with some comparative discussions. Some open problems are highlighted.

Throughout the text, successive assumptions, simplifications, and approximations are indicated separately. Successive results are explicitly indicated, for compactness, within the main text.

2. Open Classical One Particle Systems

This section will introduce techniques to be used directly in Section 3, dealing with many-particle systems.

2.1. One-Dimensional Case: Some General Aspects

Let a classical particle, with mass m, position x and momentum q, be subject to a real potential $V=V(x)$, in the presence of a “heat bath” ($hb$) at thermal equilibrium at temperature T. We shall employ the standard variable $\beta ={(k_{B}T)}^{-1}$ ($k_B$ being Boltzmann’s constant).

Assumptions: The potential is repulsive: $V(x)\ge 0$: either $V(x)\to 0$ as $\mid x\mid \to +\infty$ or $V(x)\equiv 0$.

The classical Hamiltonian of the particle is: $H_c=q^2/(2m)+V$. Let the classical particle be, at the initial time $t=0$, out of thermal equilibrium with the $hb$, and have a probability distribution function $W_{c,in}=W_{c,in}(x,q)(\ge 0)$ to be at the position x with momentum q. Then, the non-equilibrium particle could be, at time t($>0$), at the position x with momentum q, with probability distribution function $W_c=W_c(x,q;t)(\ge 0)$. For instance, the $hb$ could be air (at rest and at thermal equilibrium, at $T\simeq 300$ ${}^{\circ }$K) in a room, and the classical particle could be a virus or a grain of pollen, performing Brownian motion in air.

The time evolution of $W_c$ is given by the reversible Liouville equation:

$$\begin{array}{c}\frac{\partial W_c}{\partial t}+\frac{q}{m}\frac{\partial W_c}{\partial x}-\frac{\partial V}{\partial x}\frac{\partial W_c}{\partial q}=0\hfill \end{array}$$

(1)

with the initial condition $W_{c,in}$. The equilibrium (Boltzmann’s) distribution, the t-independent solution of Equation (1) describing thermal equilibrium of the particle with the $hb$, is: $W_c=W_{c,eq}=exp[-\beta (q^2/(2m)+V)]$, Gaussian in q.

Let $q_{eq}={[2m/\beta ]}^{1/2}$ be the physically natural fixed (scaling) momentum determined by T. We shall define $y=q/q_{eq}$. We shall introduce the denumerably infinite family of all (unnormalized) polynomials in y: $H_n=H_n(y)$ ($n=0,1,2,3,\dots$), the standard n-th Hermite polynomial [35], orthogonalized in y (for fixed x) by using $W_{c,eq}$ as (Gaussian) weight function. With $H_0(y)=1$, one has, for $n\ne n^{\prime }$ and any x (left unintegrated):

$$\begin{array}{c}{\int }_{-\infty }^{+\infty }dy{W}_{c,eq}(x,q)H_n(y)H_{n^{\prime }}(y)=0\hfill \end{array}$$

(2)

The orthonormalized polynomials are $H_n(y)/{({\pi }^{1/2}{2}^{n}n!)}^{1/2}$. We shall introduce the (normalized) non-equilibrium classical moments $W_{c,n}=W_{c,n}(x;t)$ ($n=0,1,2,\dots$) of $W_c$ [10,23,25,26]:

$$\begin{array}{c}{W}_{c,n}=W_{c,n}(x,t)=\int d(q/q_{eq})\frac{H_n(q/q_{eq})}{{({\pi }^{1/2}{2}^{n}n!)}^{1/2}}{W}_c(x,q;t)\hfill \end{array}$$

(3)

$W_{c,0}$ is the marginal probability distribution for x. If $W_c=W_{c,eq}$, then $W_{c,eq,0}={\pi }^{1/4}exp[-\beta V]$ and $W_{c,eq,n}=0$, $n=1,2,\dots$ Equation (3) can also be applied to the initial off-equilibrium distribution $W_{c,in}$ and gives the initial moments, $W_{c,in,n}$. One has the following (formal) expansion for $W_c$:

$$\begin{array}{c}{W}_c=W_{c,eq}(x,q)\sum _{n=0}^{+\infty }{W}_{c,n}(x;t)\frac{H_n(y)}{{({\pi }^{1/2}{2}^{n}n!)}^{1/2}}\hfill \end{array}$$

(4)

2.2. 3-Term Hierarchy and Operator-Continued Fractions

Equations (1)–(3) yield an exact three-term non-equilibrium hierarchy for all $W_{c,n}$’s. It is convenient to work with the symmetrized moments $g_n=g_n(x,t)=W_{c,eq,0}{(x)}^{-1/2}{W}_{c,n}(x,t)$. The hierarchy for the $W_{c,n}$’s becomes an exact three-term hierarchy for the $g_n$’s, for any n, as a result:

$$\begin{array}{c}\frac{\partial g_n}{\partial t}=-M_{c,n,n+1}{g}_{n+1}-M_{c,n,n-1}{g}_{n-1}\hfill \end{array}$$

(5)

$$\begin{array}{c}{M}_{c,n,n\pm 1}{g}_{n\pm 1}\equiv {[\frac{(n+(1/2)(1\pm 1))k_{B}T}{m}]}^{1/2}[\frac{\partial g_{n\pm 1}}{\partial x}-\frac{(\pm 1)g_{n\pm 1}}{2k_{B}T}\frac{\partial V}{\partial x}]\hfill \end{array}$$

(6)

with initial condition $g_{in,n}=W_{c,in,n}/W_{c,eq,0}^{1/2}$. One key fact is that $M_{c,n,n+1}$ and $-M_{c,n+1,n}$ are the adjoint of each other (Equations (5) and (6) providing an anti-Hermitian hierarchy).

Let us consider the Laplace transforms ${\tilde{g}}_n={\tilde{g}}_n(s)={\int }_0^{+\infty }dt{g}_{n}exp(-st)$, with inverse $g_n={\int }_{c-i\infty }^{c+i\infty }(ds/2\pi i){\tilde{g}}_{n}exp(st)$ (c being real and such that ${\tilde{g}}_n(s)$ is analytic in the half-plane $Res>c$ of the complex s-plane). This definition and Equation (5) yield the three-term hierarchy for ${\tilde{g}}_n$:

$$\begin{array}{c}s{\tilde{g}}_n=g_{in,n}-M_{c,n,n+1}{\tilde{g}}_{n+1}-M_{c,n,n-1}{\tilde{g}}_{n-1}\hfill \end{array}$$

(7)

The hierarchy Equation (7) can be solved formally by extending to it standard procedures for solving numerical three-term linear recurrence relations in terms of continued fractions (see, for instance, [10]). That formal procedure yields all ${\tilde{g}}_n(s)$, for any $n=1,\dots$, in terms of sums of products of certain s-dependent linear operators (actually, integral operators) $D[n^{\prime };s]$, $n^{\prime }\ge n$, acting upon ${\tilde{g}}_{n-1}(s)$ and upon all $g_{in,n^{\prime }}$’s, with $n^{\prime }\ge n$. The linear operators $D[n;s]$’s are defined recurrently through:

$$\begin{array}{c}D[n;s]={[sI-M_{c,n,n+1}D[n+1;s]M_{c,n+1,n}]}^{-1}\hfill \end{array}$$

(8)

I is the unit operator. By iteration of Equation (8), $D[n;s]$ becomes a formal infinite continued fraction of products of the linear operators $M_{c,n,n+1}$ and $M_{c,n+1,n}$ (which do not commute, in general). That formal infinite continued fraction of operators reads:

$$\begin{array}{c}D[n;s]=\frac{I}{sI-M_{c,n,n+1}\frac{I}{sI-M_{c,n+1,n+2}\frac{I}{sI-M_{c,n+2,n+3}\frac{I}{sI-\dots }{M}_{c,n+3,n+2}}{M}_{c,n+2,n+1}}{M}_{c,n+1,n}}\hfill \end{array}$$

(9)

Simplification—for a simpler hierarchy and a clearer exposition, without loss of generality, let us assume that $W_{c.in,n^{\prime }}=0$ for $n^{\prime }\ge 1$, with $W_{c,in,0}\ne 0$.

By solving Equation (7) for $n(\ge 1)$ (through Equation (8)), we get the following result, without approximations, so far:

$$\begin{array}{c}s{\tilde{g}}_0=W_{c,eq,0}^{-1/2}{W}_{c,in,0}-M_{c,0,1}{\tilde{g}}_1\hfill \end{array}$$

(10)

$$\begin{array}{c}{\tilde{g}}_n(s)=-D[n;s]M_{c,n,n-1}{\tilde{g}}_{n-1}(s),n\ge 1\hfill \end{array}$$

(11)

2.3. Properties of $D[n;s]$

Let us choose $n(\ge 1)$ and fix $s=ϵ>0$ (real and suitably small) in any $D[n;s]$. Then, the following crucial properties appear to hold (in general, for either $V\ne 0$ or $V=0$) [25,26]: if $D[n+1;ϵ]$ were Hermitian and if all its eigenvalues (which would be real) were non-negative, then the same would hold true for $D[n;ϵ]$. It is easy to confirm the validity of that property if $D[n+1;ϵ]$ is replaced by a $2\times 2$ matrix with the above (non-negativity and Hermiticity) properties and $M_{n,n+1}$ and $M_{n+1,n}$ are chosen as $2\times 2$ matrices such that $M_{n,n+1}$ be the adjoint of $-M_{n+1,n}$. In particular, a confirming example with a $2\times 2$ matrix is given in Appendix D in [26]. Then, through iterative arguments, as a result, $D[n;ϵ]$’s, $n=0,1,2,3\dots$, turn out to be Hermitian, and all their eigenvalues are non-negative, for $V\ne 0$, in general.

To understand the structure and convergence of the operator-continued fractions $D[n;s]$, we shall take $V\equiv 0$ and let $ϵ>0$ first, so as to allow for $ϵ\to 0$ later [25]. Let us perform a spatial Fourier transformation from configuration space (x) to wavevector space (k), by applying ${(2\pi )}^{-1/2}\int dxexp(-ikx)$. Let $e(k)\equiv {(2m)}^{-1}{k}_{B}T{k}^2$. Then, the Fourier transform of the operator-continued fraction in Equation (8), for Re$s>0$ and $n\ge 0$ is the following continued fraction involving ordinary functions $D_1[k;n;s]$:

$$\begin{array}{c}{D}_1[k;n;s]=[s+2e(k)(n+1)D_1{[k;n+1;s]}^{-1}\hfill \end{array}$$

(12)

By iteration, $D_1[k;n;s]$ becomes the following ordinary continued fraction (involving no non-commuting operators):

$$\begin{array}{c}{D}_1[k;n;s]=\frac{1}{s+2e(k)(n+1)\frac{1}{s+2e(k)(n+2)\frac{1}{s+2e(k)(n+3)\frac{1}{s+2e(k)(n+4)\frac{1}{s+\dots \dots ..}}}}}\hfill \end{array}$$

(13)

to be compared to Equation (9) (involving non-commuting operators).

It will be more convenient to make use of (shorter) standard notations [36] for ordinary infinite continued fractions (involving only commuting quantities). Accordingly, Equation (13) reads, in standard notations:

$$\begin{array}{c}{D}_1[k;n;s]=\frac{1}{s+}\frac{2e(k)(n+1)}{s+}\frac{2e(k)(n+2)}{s+}\dots =\frac{1}{e{(k)}^{1/2}}.[\frac{2^{-1}}{z+}\frac{2^{-1}(n+1)}{z+}\frac{2^{-1}(n+2)}{z+}\dots ].\hfill \end{array}$$

(14)

with $z=s/(2e{(k)}^{1/2})$. On the other hand, the continued fraction can be expressed as:

$$\begin{array}{c}\frac{2^{-1}}{z+}\frac{2^{-1}(n+1)}{z+}\frac{2^{-1}(n+2)}{z+}\dots =\frac{i^{n}erfc(z)}{i^{n-1}erfc(z)}.\hfill \end{array}$$

(15)

$i^{n}erfc(z)$ being the n-th repeated integral of the complementary error function ($i^{0}erfc(z)=erfc(z)$) [36]. As z varies in $0\le z<+\infty$, the values given by Equation (15) vary in a finite interval. With $ϵ>0$ and by using [36], one gets: $D_1[k;1;ϵ]\to {ϵ}^{-1}$ as $k\to 0$, while $D_1[k;1;ϵ]\to {(\pi e(k))}^{-1/2}$ as $\mid k\mid \to \infty$. For $s=ϵ=0$, the behavior of $D_1$ is different. Equation (14) gives:

$$\begin{array}{c}{D}_1[k;n;s=ϵ=0]=D_1[k;n;0]={[2e{(k)}^{1/2}]}^{-1}\frac{1}{0+}\frac{2^{-1}(n+1)}{0+}\frac{2^{-1}(n+2)}{0+}\dots \hfill \end{array}$$

(16)

The continued fraction in Equation (16) can, in turn, be evaluated in terms of the standard Gamma function $\mathsf{\Gamma }$ [36]:

$$\begin{array}{c}\frac{2^{-1}}{0+}\frac{2^{-1}(n+1)}{0+}\frac{2^{-1}(n+2)}{0+}\dots .=\frac{\mathsf{\Gamma }((n/2)+1/2)}{2\mathsf{\Gamma }((n/2)+1)}\hfill \end{array}$$

(17)

The ratio in Equation (17) behaves as $n^{-1/2}$, as $n\to +\infty$. We state that the various classical operator-continued fractions will also behave, on average, as $n^{-1/2}$ for large n: we shall omit details [25]. For $k\to 0$, $D_1[k;n;0]$ diverges as $k^{-1}$ (due to $e{(k)}^{-1/2}$). Then, $\int dk{D}_1[k;n;0]$ also diverges near $k=0$.

2.4. An Approximate Ansatz for Equation (8), for Large n

Using the developments in Section 2.3 as insights, we shall now outline an approximate ansatz, in order to simplify Equation (8) under the following conditions.

Approximation—Let $s=ϵ>0$ be real and suitably small and $n_0^{\prime }$ be a suitably large positive integer. Then, for $n,n^{\prime },n^{″}>n_0^{\prime }$ and as an approximation, we disregard the non-commutativities of all $M_{c,n,n+1}$, $M_{c,n+1,n}$ and $D[n^{″};ϵ]$ among themselves and regard all of them as commuting operators.

Notice also that the non-negative Hermitean operator $(-M_{c,n,n+1}{M}_{c,n+1,n})/(n+1)$ is independent of n and that the spectrum of all eigenvalues of $(-M_{c,n,n+1}{M}_{c,n+1,n})$ varies in $(0,+\infty )$. Then, after a direct formal manipulation inspired on the structures of Equations (14) and (16), Equation (8) becomes:

$$\begin{array}{c}D[n;ϵ]={[2(-M_{c,n,n+1}{M}_{c,n+1,n})/(n+1)]}^{-1/2}\frac{1}{z_1+}\frac{2^{-1}(n+1)}{z_1+}\frac{2^{-1}(n+2)}{z_1+}\dots ].\hfill \end{array}$$

(18)

with $z_1=ϵ{[2(-M_{c,n,n+1}{M}_{c,n+1,n})/(n+1)]}^{-1/2}$. A comparison with Equation (14) shows that Equation (18) is exact if $V\equiv 0$.

Approximation—At this point, for suitably small $ϵ$, we assume that we can neglect the contribution of the non-negative Hermitean operator $z_1$.

By using Equations (16) and (17), we get the approximate ($ϵ$-independent) operator, as a result:

$$\begin{array}{c}D[n;ϵ]={[2(-M_{c,n,n+1}{M}_{c,n+1,n})/(n+1)]}^{-1/2}\frac{\mathsf{\Gamma }((n/2)+1/2)}{\mathsf{\Gamma }((n/2)+1)}\hfill \end{array}$$

(19)

One expects that the dominant contribution to this approximate $D[n;ϵ]$ is determined by the set of lowest (non-negative) eigenvalues of $(-M_{c,n,n+1}{M}_{c,n+1,n})$.

2.5. Long-Term Approximation

Approximations—The long-term approximation for $n\ge n_0\ge 1$ reads in general (that is, without imposing $W_{c.in,n^{\prime }}=0$ for $n^{\prime }\ge 1$, with $W_{c,in,0}\ne 0$) as follows. One replaces for $n\ge n_0\ge 1$ any $D[n;s]$ yielding ${\tilde{g}}_n$ in terms of ${\tilde{g}}_{n-1}$ by $D[n;ϵ]$: this approximation is not done for $n<n_0$, which will be crucial, and is the better, the larger $n_0$. We regard $D[n_0;ϵ]$ for $n=n_0$ as a fixed (s-independent) operator. Some ansatz or approximation should be provided directly for $D[n_0;ϵ]$, like the one in Section 2.4. For simplicity, let us continue with the simplification: $W_{c.in,n^{\prime }}=0$ for $n^{\prime }\ge 1$, with $W_{c,in,0}\ne 0$. Moreover, after the above long-term approximation, we shall continue with the same initial condition $W_{c,in,0}$ at $t=0$: it may well be that this amounts to another kind of approximation.

Then, for small s, we approximate for $n=n_0$ as: ${\tilde{g}}_{n_0}(s)\simeq -D[n_0;ϵ]M_{n_0,n_0-1}{\tilde{g}}_{n_0-1}(s)$. The resulting hierarchy for ${\tilde{g}}_n$’s ($n=0,\dots ,n_0-1$), through the inverse Laplace transform, yields a closed approximate irreversible hierarchy for $g_n$, $n=0,1,\dots ,n_0-1$. The solutions of the last closed hierarchy for $g_n$ relax irreversibly, for large t and reasonable $W_{c,in,0}$, towards $W_{c,eq,0}^{1/2}\ne 0$ and $W_{c,eq,n}^{1/2}=0$, $n=1,\dots ,n_0-1$ (thermal equilibrium) [25,26]. Then, for the long term, the dominant moment is $g_0$, while any $g_n$ with $n>0$ being the smaller, the larger n and $t(>0)$ are (due to the behaviors of $M_{c,n,n\pm 1}$ and of $D[n;ϵ]$ with n). Similar behaviors hold for $W_{c,0}$ and $W_{c,n}$ with $n>0$.

Simplification—We shall illustrate the above facts by taking, for simplicity, $n_0=1$.

Then, Equation (11) for small s and $n=1$ becomes: ${\tilde{g}}_1(s)\simeq -D[1;ϵ]M_{c,1,0}{\tilde{g}}_0(s)$ (to be compared to Equation (11) before the approximation). Then, Equations (10) and (11) yield, by taking inverse Laplace transforms:

$$\begin{array}{c}\frac{\partial g_0}{\partial t}=-M_{c,0,1}{g}_1\hfill \end{array}$$

(20)

$$\begin{array}{c}{g}_1\simeq -D[1;ϵ]M_{c,1,0}{g}_0\hfill \end{array}$$

(21)

One finds directly the irreversible Smoluchowski-like equation for the $n=0$ moment, as a result:

$$\begin{array}{c}\partial g_0/\partial t=M_{c,0,1}D[1;ϵ]M_{c,1,0}{g}_0\hfill \end{array}$$

(22)

with initial condition $W_{c,eq,0}^{-1/2}{W}_{c,in,0}$.

The right-hand side of Equation (22) should be interpreted as

${\int }_{-\infty }^{+\infty }d{x}^{\prime }(M_{c,0,1}D[1;ϵ]M_{c,1,0})(x,x^{\prime })g_0(x^{\prime };t)$. Let: $(f_1,f_2)={\int }_{-\infty }^{+\infty }dx{f}_1{(x)}^{*}{f}_2(x)$ for suitable functions $f_1$ and $f_2$. Due to the Hermiticity of $D[1;ϵ]$: $(f_1,M_{c,0,1}D[1;ϵ]M_{c,1,0}{f}_2)=(M_{c,0,1}D[1;ϵ]M_{c,1,0}{f}_1,f_2)$, thereby checking that $M_{c,0,1}D[1;ϵ]M_{c,1,0}$ is Hermitian. Moreover:

$(f_1,M_{c,0,1}D[1;ϵ]M_{c,1,0}{f}_1)=-(M_{c,1,0}{f}_1,D[1;ϵ]M_{c,1,0}{f}_1)\le 0$ for arbitrary functions $f_1$, as all eigenvalues of $D[1;ϵ]$ are $\ge 0$. Let $f_{\lambda }(x)$ be an eigenfunction of the integral operator $M_{c,0,1}D[1;ϵ]M_{c,1,0}$ with eigenvalue $\lambda (\ge 0)$. Then, $(M_{c,0,1}D[1;ϵ]M_{c,1,0})(x,x^{\prime })={\sum }_{\lambda }\lambda f_{\lambda }(x)f_{\lambda }{(x^{\prime })}^{*}$. ${\sum }_{\lambda }$ is a short-hand notation denoting integration and summation over the whole spectrum of $M_{c,0,1}D[1;ϵ]M_{c,1,0}$. By expanding $W_{c,eq,0}^{-1/2}{W}_{c,in,0}={\sum }_{\lambda }{g}_{in,\lambda }{f}_{\lambda }(x)$, with x-independent $g_{in,\lambda }$, the solution of (22) with the above initial condition is $g_0={\sum }_{\lambda }{g}_{in,\lambda }{f}_{\lambda }(x)exp(-\lambda t)$, which relaxes irreversibly as $t\to +\infty$ towards $g_{in,0}{f}_0(x)$, corresponding to $\lambda =0$. At equilibrium, one has: $g_0=W_{c,eq,0}^{1/2}$ (proportional to $f_0$), $M_{c,1,0}{g}_0=0$ and $g_n=0$, $n=1,2,3\dots$, consistently. Then, Equation (22) is (at least, with $ϵ>0$) as irreversible as the standard heat equation: for long $t(>0)$ the dominant moment is $g_0$, while any $g_n$ with $n>0$ is negligible, $g_n$ being the smaller, the larger n and $t(>0)$, and so on for the $W_{c,n}$’s.

Irreversible thermalization does not happen in the absence of long-term approximations. Then, the above approximations ($D[n;s]\simeq D[n;ϵ]$ for $n\ge n_0$, but not for $n<n_0$), giving rise to the thermalization with the $hb$ implemented in Equation (22), is an alternative way of introducing irreversibility out of the reversible Equation (1).

3. Closed Classical Many-Particle Systems: Long-Term Approximation and Arrow of Time

3.1. Initial State Motivated by Fluid Dynamics, Hierarchy, and Continued Fractions

We shall outline the main developments (details being given in [25]). We treat a closed large system of many ($N\gg 1$) classical non-relativistic particles, in three spatial dimensions, with spatial coordinates ${\mathbf{x}}_1$,...,${\mathbf{x}}_N$ ($\equiv (\mathbf{x})$) and momenta ${\mathbf{q}}_1$,...,${\mathbf{q}}_N$ ($\equiv (\mathbf{q})$). Let $x_{i,\alpha }$ and $q_{i,\alpha }$ be the Cartesian components of ${\mathbf{x}}_i$ and ${\mathbf{q}}_i$, respectively ($i=1,\dots ,N$, $\alpha =1,2,3$). The non-equilibrium classical probability distribution function is: $W_c=W_c((\mathbf{x}),(\mathbf{q});t)$.

Assumptions—(a) All particles, which are identical, have mass m. Neither a $hb$ nor external friction mechanisms nor external forces are assumed. The physical (qualitative) idea is that the system contains a very large set ($s_1$) of degrees of freedom at large distances at thermal equilibrium at absolute temperature T with one another: $s_1$ plays the role of an (internal) $hb$. However, by assumption the system is in a non-equilibrium state, because it also contains a large set of degrees of freedom ($s_2$) at finite distances, which are off-equilibrium with the previous set $s_1$, and among themselves. We assume that the set $s_1$ is larger than the set $s_2$. Then, the initial non-equilibrium distribution $W_{c,in}=W_{c,in}((\mathbf{x}),(\mathbf{q}))$, which is known by assumption, describes thermal equilibrium at absolute temperature T for the set $s_1$ located at large distances, but off-equilibrium for the set $s_2$ located at finite distances [25]. The above qualitative idea will become clearer through $W_{c,in}$ and its properties, discussed below in (c).

(b) The interaction potential among the particles is: $V={\mathsf{\Sigma }}_{i,j=1,i<j}^N{v}_{ij}(\mid {\mathbf{x}}_i-{\mathbf{x}}_j\mid )$ and we suppose that all $v_{ij}(\mid {\mathbf{x}}_i-{\mathbf{x}}_j\mid )$ are repulsive ($\ge 0$) and tend quickly to zero for large $\mid {\mathbf{x}}_i-{\mathbf{x}}_j\mid$. $H_{c,N}={(2m)}^{-1}{\mathsf{\Sigma }}_{i=1}^N{\mathsf{\Sigma }}_{\alpha =1}^3{q}_{i,\alpha }^2+V$ is the classical N-particle Hamiltonian. Then, Boltzmann’s equilibrium (canonical) distribution at temperature T is: $W_{c,eq}=exp[-\beta H_{c,N}]$.

(c) The initial distribution function $W_{c,in}$, at $t=0$, will be chosen to be: (i) qualitatively consistent with the idea [3,37], typical of Information Theory that one should employ only distribution functions compatible with the limited information available which, in turn, refers to expectation values of only a subset of observable dynamical variables (these ideas being imposed for $t=0$ only, but not for $t>0$); (ii) also consistent with standard variables employed in equilibrium statistical mechanics and fluid dynamics [2,7,8,9]. Then, our ansatz for $W_{c,in}$ will depend on a finite number (actually, $2+3=5$) of functions of one single position vector, $\mathbf{x}$: ${\lambda }_k={\lambda }_k(\mathbf{x})$, $k=0,2$, and ${\lambda }_{1,\alpha }={\lambda }_{1,\alpha }(\mathbf{x})$, $\alpha =1,2,3$ (all of them being independent of time and on momenta). The expression for $W_{c,in}$ in terms of ${\lambda }_k$ and ${\lambda }_{1,\alpha }$ has appeared previously [2,9] (and is related to the Massieu-Planck function [9]). $W_{c,in}$, which describes thermal equilibrium with homogeneous temperature T for large distances (large $\mid \mathbf{x}\mid$) but non-equilibrium for intermediate and short distances $\mid \mathbf{x}\mid$ (with spatial inhomogeneities), reads:

$$\begin{array}{c}{W}_{c,in}={(N!)}^{-1}exp[-\int d^3\mathbf{x}({\lambda }_0(\mathbf{x})\sum _{i=1}^N{\delta }^{(3)}({\mathbf{x}}_i-\mathbf{x})+\sum _{\alpha =1}^3{\lambda }_{1,\alpha }(\mathbf{x})\sum _{i=1}^N\frac{q_{i,\alpha }}{m}\hfill \\ \times {\delta }^{(3)}({\mathbf{x}}_i-\mathbf{x})+{\lambda }_2(\mathbf{x})\sum _{i=1}^N(\frac{{\mathbf{q}}_i^2}{2m}+\frac{1}{2}\sum _{j=1,j\ne i}^N{v}_{ij}(\mid {\mathbf{q}}_i-{\mathbf{q}}_j\mid )){\delta }^{(3)}({\mathbf{x}}_i-\mathbf{x})].\hfill \end{array}$$

(23)

${\delta }^{(3)}$ denoting the three-dimensional Dirac delta function. Consistently with (i)–(ii) above, the $\lambda$’s will be uniquely determined in terms of $2+3=5$ $\mathbf{x}$-dependent observables (also independent of time and on momenta) typically employed in Fluid Dynamics, which, by assumption, are known at $t=0$: mass density, fluid velocity and some suitable energy density [9,25]. How does the equilibrium temperature T appear in this closed system, without $hb$? We accept that ${\lambda }_2(\mathbf{x})$ approaches quickly a non-vanishing constant, ${\lambda }_2(\infty )$, as $\mid \mathbf{x}\mid$ tends to ∞ along any direction and that the same holds for ${\lambda }_0(\mathbf{x})$. A similar statement holds for ${\lambda }_{1,\alpha }(\mathbf{x})$, the corresponding (large-$\mid \mathbf{x}\mid$) limiting value being zero. At finite $\mathbf{x}$, the off-equilibrium ${\lambda }_2(\mathbf{x})$, ${\lambda }_0(\mathbf{x})$ and ${\lambda }_{1,\alpha }(\mathbf{x})$ do depend on $\mathbf{x}$ and, so, differ from their respective constant (large-$\mid \mathbf{x}\mid$) limiting values, which describe equilibrium. Then, consistency is achieved (T being thereby introduced) if, in the thermodynamical limit, ${\lambda }_2(\infty )$ tends to ${(k_{B}T)}^{-1}$(plus corrections which approach zero in that limit). The dominant contributions to various statistical averages at $t=0$ come from large $\mathbf{x}$, up to corrections which decrease as N increases (and tend to vanish as $N\to \infty$). For a detailed analysis, see [25].

The reversible Liouville equation reads:

$$\begin{array}{c}\frac{\partial W_c}{\partial t}={\mathsf{\Sigma }}_{i=1}^N{\mathsf{\Sigma }}_{\alpha =1}^3[\frac{\partial V}{\partial x_{i,\alpha }}\frac{\partial W_c}{\partial q_{i,\alpha }}-\frac{q_{i,\alpha }}{m}\frac{\partial W_c}{\partial x_{i,\alpha }}]\hfill \end{array}$$

(24)

Let $[n]$ denote a set of non-negative integers $(n(i=1,\alpha =1),\dots ,n(i=N,\alpha =3))$ and let $n={\mathsf{\Sigma }}_{l=1}^N{\mathsf{\Sigma }}_{\alpha =1}^{3}n(l,\alpha )$. Let $\left[dq\right]={\prod }_{i=1}^N{\prod }_{\alpha =1}^{3}d{q}_{i,\alpha }$. We introduce non-equilibrium moments $W_c([n])$ of W (using products of Hermite polynomials, by generalizing Equation (3)):

$$\begin{array}{c}\int \left[dq\right]\prod _{i=1}^N\prod _{\alpha =1}^3\frac{H_{n(i,\alpha )}(q_{i,\alpha }/{(2m{k}_{B}T)}^{1/2})}{{({\pi }^{1/2}{2}^{n(i,\alpha )}n(i,\alpha )!)}^{1/2}}{W}_c((\mathbf{x}),(\mathbf{q}),t)\equiv W_c((\mathbf{x});[n];t)=W_c([n]),\hfill \end{array}$$

(25)

If $W_c=W_{c,eq}$, then $W_{c,eq}([0])$ ($[0]=(0,0,,,0)$) is proportional to $exp[-\beta V]$ and $W_{c,eq}([n])=0$, $[n]\ne [0]$ (say, $n\ne 0$). Equation (25) can also be applied to $W_{c,in}$ and gives the corresponding initial moments, $W_{c,in}([n])$. We shall work with the symmetrized moments $g([n])=W_{c,eq}{([0])}^{-1/2}{W}_c([n])$.

One gets an infinite reversible three-term linear recurrence for $g([n])$’s, generalizing Equations (5) and (6). It reads, as a result:

$$\begin{array}{c}\frac{\partial g(n(1,1),\dots ,n(j,\beta ),\dots ,n(N,3))}{\partial t}=-{\mathsf{\Sigma }}_{l=1}^N{\mathsf{\Sigma }}_{\alpha =1}^3[M_{l,\alpha ;n(l,\alpha );+}g(n(1,1),\dots ,n(l,\alpha )+1,\dots ,n(N,3))\\ +M_{l,\alpha ;n(l,\alpha );-}g(n(1,1),\dots ,n(l,\alpha )-1,\dots ,n(N,3))]\end{array}$$

(26)

$$\begin{array}{c}{M}_{l,\alpha ;n(l,\alpha );+}={[\frac{(n(l,\alpha )+1)k_{B}T}{m}]}^{1/2}[\frac{\partial }{\partial x_{l,\alpha }}-\frac{1}{2k_{B}T}\frac{\partial V}{\partial x_{l,\alpha }}]\hfill \end{array}$$

(27)

$$\begin{array}{c}{M}_{l,\alpha ;n(l,\alpha );-}={[\frac{n(l,\alpha )k_{B}T}{m}]}^{1/2}[\frac{\partial }{\partial x_{l,\alpha }}+\frac{1}{2k_{B}T}\frac{\partial V}{\partial x_{l,\alpha }}]\hfill \end{array}$$

(28)

The Laplace transform of the hierarchy Equation (26), which is the actual many-particle counterpart of Equation (7), can be formally solved in terms of linear operators $D[[n];s]$, which generalize the previous $D[n;s]$. For details, see [25]. All $D[[n];s]$ are square matrices, due to the indices $i=1,\dots ,N$ and $\alpha =1,\dots ,3$. In turn, each matrix element in those square matrices is an integral operator, arising from the partial differential operators $M_{l,\alpha ;n(l,\alpha );+}$ and $M_{l,\alpha ;n(l,\alpha );-}$, as l, $\alpha$ and $n(l,\alpha )$ vary. The $D[[n];s]$ fulfill the following formal hierarchy (which generalizes Equation (8)):

$$\begin{array}{c}D[[n];s]={[s-M_{+,[n+1]}D[[n+1];s]M_{-,[n]}]}^{-1}\hfill \end{array}$$

(29)

The linear operators $M_{\pm ,[n]}$ are rectangular matrices, the elements of which are formed out of the partial differential operators $M_{l,\alpha ;n(l,\alpha );+}$ and $M_{l,\alpha ;n(l,\alpha );-}$. $M_{+,[n+1]}$ can be shown to be the adjoint of $-M_{-,[n]}$ [25]. By iterating Equation (29) indefinitely, one can express formally the linear operator $D[[n];s]$ as an operator-continued fraction, which depends on all partial differential operators $M_{l,\alpha ;n(l,\alpha );+}$ and $M_{l,\alpha ;n(l,\alpha );-}$ and generalizes the operator-continued fraction for $D[n;s]$.

By generalizing the iterative arguments in Section 2.3, it follows that $D[[n];ϵ]$ for $ϵ>0$ is a Hermitian operator with non-negative eigenvalues for $n\ge n_0\ge 1$ [25] (in general, for either $V\ne 0$ or $V=0$), as a result. No approximation has been performed so far.

Simplification—A few remarks for the case $V\equiv 0$ may be clarifying. In the Laplace transform of Equation (26), let us perform a Fourier transformation from configuration space (${\mathbf{x}}_1$, ..., ${\mathbf{x}}_N$) to wavevector space (${\mathbf{k}}_1$, ..., ${\mathbf{k}}_N$)$\equiv [k]$. Let $e(k;N)\equiv k_{B}T{\sum }_{j=1}^N{(2m)}^{-1}{\mathbf{k}}_j^2$. Then, the Fourier transform $D_1[[k];[n];s]$ of $D[[n];s]$ for Re$s>0$ is an ordinary continued fraction, given in Equations (12)–(14) with $e(k)$ replaced by $e(k;N)$. Then, with such a replacement, the properties of $D_1(k;n;s)$ given in Section 2.4 also hold for $D_1[[k];[n];s]$. Notice that $D_1[[k];[n];s=0]$ diverges as $e{(k;N)}^{-1/2}$ if $e(k;N)\to 0$. On the other hand, and contrary to what happened for one particle in one spatial dimension (recall the comment after Equation (17)), $\int \left[dk\right]D_1[[k];[n];s=0]$ converges near $e(k;N)=0$. It is open whether actual counterpart of Equation (11) (containing $D[[n];ϵ]$ acting on various $M_{l,\alpha ;n(l,\alpha );-}g(n(1,1),\dots ,n(l,\alpha )-1,\dots ,n(N,d))$’s), when integrated over $[k]$’s to come back to $(\mathbf{x})$-space, would converge at small $[k]$’s as $ϵ\to 0$, when $V\ne 0$ and all $v_{ij}(\mid {\mathbf{x}}_i-{\mathbf{x}}_j\mid )$ are repulsive and tend quickly to zero for large distances. The approximate ansatz in Section 2.4 can be extended readily to $D[[n_0];ϵ]$.

3.2. Long-Term Approximation, Irreversibility, Liapunov Function and BBGKY Hierarchy

Approximations—Despite the very involved structure of Equation (26), we argue that a simple long-term approximation can be performed in it, for $v_{ij}\ge 0$ (and vanishing quickly at large distances) and very large N (eventually, in the thermodynamical limit), which generalizes directly the one in Section 2.5. This approximation consists in fixing $s=ϵ>0$ ($ϵ$ being small) in the whole hierarchy of operators $D[[n];s]$, for any $n(={\mathsf{\Sigma }}_{l=1}^N{\mathsf{\Sigma }}_{\alpha =1}^{3}n(l,\alpha ))\ge n_0(>0)$ which, then, become Hermitian operators $D[[n];ϵ]$ with no negative eigenvalues. It is crucial that s-dependences be kept in $D[[n];s]$, for $n<n_0$. Notice that the non-vanishing factors $n{(l,\alpha )}^{1/2}$ and ${((n(l,\alpha )+1))}^{1/2}$ in $M_{l,\alpha ;n(l,\alpha );-}$ and $M_{l,\alpha ;n(l,\alpha );+}$, respectively, tend to reduce, as the $n(l,\alpha )$’s increase, the relative importance of having fixed $s=ϵ$ and the contribution of the latter in the $D[[n];s=ϵ]$’s with $n\ge n_0$. This is a genuine feature of the $D[[n];s=ϵ]$’s. See [25], where it was shown that by imposing $n_0>2$, the long-term approximation is still exactly consistent with all (five) hydrodynamical balance equations. For simplicity, we discard all the initial moments $W_{c,in}([n])$ for $n\ge n_0$. We regard $D[[n_0];ϵ]$ as a fixed (s-independent) operator, yielding all $g([n_0])$ in terms of all $g([n_0-1])$. Moreover, after the above long-term approximation, we shall continue with the same initial condition $W_{c,in}([0])$ at $t=0$: it may amount to another kind of approximation.

All that leads to a closed approximate hierarchy for $g([n])$’s (with initial moments $W_{c,in}([n])$), with $n<n_0$, which appears to yield an approximate irreversible evolution towards thermal equilibrium at T. $g([n])$’s and, then, $W_c([n])$ relax the quicker the larger n, provided that $n<n_0$. $W_c([0])$ would dominate the approach towards equilibrium for $t\to +\infty$. See [25]. All that appears to work for fixed $D[[n_0];ϵ]$. An arrow of time would follow approximately in the present case, as discussed below.

Simplification—As an extreme example, let $n_0=1$ (which is strictly consistent only with the hydrodynamical balance equations for mass): see [25] for $n_0=2,3$. By making the above long-term approximation and taking inverse Laplace transforms, one finds directly the irreversible Smoluchowski-like equation for the $[n=0]$ moment, which generalizes (22) ($[n=0]$ meaning $n(1,1)=0,\dots ,n(j,\beta )=0,\dots ,n(N,3)=0$):

$$\begin{array}{c}\frac{\partial g([n=0])}{\partial t}={\mathsf{\Sigma }}_{l=1}^N{\mathsf{\Sigma }}_{\alpha =1}^3{M}_{l,\alpha ;n(l,\alpha )=0;+}\times \hfill \\ ({\mathsf{\Sigma }}_{l^{\prime }=1}^N{\mathsf{\Sigma }}_{{\alpha }^{\prime }=1}^3{[D[[n=1];ϵ]]}_{l,\alpha ;l^{\prime },{\alpha }^{\prime }}{M}_{l^{\prime },{\alpha }^{\prime };n(l^{\prime },{\alpha }^{\prime })=1;-})g([n=0])\hfill \end{array}$$

(30)

as a result. The operator $D[[n=1];ϵ]$ (Hermitian, with non-negative eigenvalues) has, as a square matrix, the matrix elements ${[D[[n=1];ϵ]]}_{l,\alpha ;l^{\prime },{\alpha }^{\prime }}$. The initial condition is $W_{c,eq}{([0])}^{-1/2}{W}_{c,in}([0])$.

Discussion and consequences—Both the exact hierarchy for the $g([n])$’s and the closed approximate one for them after the long-term approximation are genuinely different from the non-equilibrium classical BBGKY hierarchy [7,8]. In fact, in the latter, in the equation for the distribution function for n particles, one leaves unintegrated their position vectors and momenta, while those for the remaining $N-n$ particles are integrated over. Moreover, such an equation also depends on the distribution function for $n+1$ particles but not on that for $n-1$ ones, a feature which, beyond the approximate framework of the standard Boltzmann equation (for $n=1$) [2,7,8], does not seem to shed much light on the long-term approach to thermal equilibrium for larger n. By contrast, in the equation for $g([n])$ in the actual non-equilibrium hierarchy based upon $W_{c,eq}$, the contributions from $g([n+1])$’s are neatly different from (for large n, approximately smaller by a factor $\simeq n^{-1/2}$ than) those coming from $g([n-1])$’s, at least in the long-term approximation [25].

The structure of (30), with $D[[n=1];ϵ]$ replaced by a constant is similar to that of the linear Smoluchowski equation in the standard Rouse model for polymer dynamics [38].

One leading purpose of this section is to introduce the t-dependent Liapunov function, for generic $n_0\ge 1$:

$$\begin{array}{c}L=\frac{1}{2}\sum _{n=0}^{n_0-1}\int \prod _{i=1}^N\prod _{\alpha =1}^{3}d{x}_{i,\alpha }g{([n])}^2\hfill \end{array}$$

(31)

The integration over any $x_{i,\alpha }$ is performed in $-\infty <x_{i,\alpha }<+\infty$. One has:

$$\begin{array}{c}\partial L/\partial t\le 0\hfill \end{array}$$

(32)

This follows readily for $n_0=1$ by using Equation (30) and the fact that the operator $D[[n=1];ϵ]$ is Hermitian, with non-negative eigenvalues. The generalization for $n_0>1$, under the long-term approximation, proceeds by employing Sections 4 and 5 in [25]: that generalization is direct, although rather cumbersome (N particles in three dimensions and higher order moments being involved), and will be omitted. Then, the actual L is a nondecreasing function in the time evolution. Equations (31) and (32) constitute other results. At equilibrium:

$$\begin{array}{c}L=\frac{1}{2}\int \prod _{i=1}^N\prod _{\alpha =1}^{3}d{x}_{i,\alpha }{q}_{eq}^{3N}exp[-\beta V]\equiv L_{eq}\hfill \end{array}$$

(33)

which is proportional to the standard exponential of the equilibrium entropy. There has been much interesting work aimed at defining, in different frameworks (thermodynamical, statistical...), a t-dependent entropy $S(t)$ in non-equilibrium phenomena, in such a way that the inequality $dS(t)/dt\ge 0$ would hold. The aim would be to characterize in a general setting, the notion of evolution. For updated accounts, see [20]. It appears that at present, a different number of non-equilibrium entropies have been proposed (and do the job) in various frameworks separately [39]. However, to the best of the present author’s knowledge, no general agreement has been reached on a unique non-equilibrium entropy $S(t)$ valid for different frameworks, and an accepted definition of it based on statistical mechanics appears also to be lacking, so far. Both in case that a general definition of non-equilibrium entropy does not exist or in case that it could indeed be formulated at the end (although it be unknown at present), a partial way out towards the characterization of an arrow of time is the following. Once the long-term approximation has been carried out as indicated above, we have obtained a t-dependent function L through Equation (31), the variation of which seems to be adequate to define an arrow of time, by virtue of Equation (32).

4. Quantum Particles: One Particle

In recent times, various rewarding results have been obtained for a system subject to a larger “heat bath” in quantum mechanics, under various conditions. Some of those results (which are related to and support our present work) are: (i) derivation of the canonical density operator for the system for long time [40], and for a suitable overwhelming majority of wave functions (typicality) [41]; (ii) the system will approach an equilibrium state and remain close to it for almost all times, independently both on the $hb$ and, under suitable conditions, also on the initial state of the system [42]. For equilibration of an isolated system, see [43]. Those results required to consider the quantum states of the system and of the $hb$ altogether and rely on their reciprocal entanglement. For further works strengthening those results, see [44,45,46,47] and references therein. For general aspects on quantum open systems, see also [2,8,9,13,14,15,17,18,19] and references therein.

General assumption and simplification—Throughout our quantum 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$ at absolute temperature T, as studied in [40,41,42,43,44,45,46,47]. To simplify matters, we will focus exclusively on the quantum states of the system considered.

Our study here of the quantum one particle case subject to an external potential will be addressed towards displaying analogies and genuine differences with respect to the preceding classical case (Section 2, mostly) and to provide some hints for the quantum N particle case (Section 5), and for those with $N=2,3$ (Section 6 and Section 7).

4.1. General Aspects

We shall consider one non-relativistic quantum spinless Brownian particle of mass $m(>0)$ and momentum operator $-i\hslash (\partial /\partial x)$, in one spatial dimension x, with (Hermitian) quantum Hamiltonian (ℏ being Planck’s constant):

$$\begin{array}{c}H=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}+V\hfill \end{array}$$

(34)

Assumptions—The real potential $V(x)$ fulfills: $V(x)\to 0$ quickly, as $\mid x\mid \to +\infty$ and $V(x)$ and all $d^{n}V(x)/d{x}^n$, $n=1,2,3,\dots$, are continuous for any x. Also, having in mind possible generalizations, we suppose that $V(x)=V(-x)$ and allow for $V(x)<0$ and, hence, the possibility of bound states of the particle by the potential. The particle is also separately subject to the $hb$.

The non-equilibrium statistical evolution for $t>0$ is given by the density operator $\rho =\rho (t)$ (a statistical mixture of quantum states), with the initial condition $\rho (t=0)={\rho }_{in}$. $\rho (t)$ for $t>0$ and ${\rho }_{in}$ are Hermitian and positive-definite linear operators acting in the Hilbert space spanned by the set of all eigenfunctions ${\phi }_j(x)$ of H. Unless otherwise stated, we shall not impose that $\rho (t)$ be normalized. One has ($[H,\rho ]=H\rho -\rho H$):

$$\begin{array}{c}\frac{\partial \rho }{\partial t}=\frac{1}{i\hslash }[H,\rho ]\hfill \end{array}$$

(35)

We consider the matrix element $\langle x-y|\rho (t)|x+y\rangle$ of $\rho (t)$ in generic eigenstates, $|x-y\rangle$, $|x+y\rangle$, of the quantum position operator. The quantum Wigner function $W=W(x,q;t)$, determined by $\rho$, is [48,49,50,51,52]:

$$\begin{array}{c}W(x,q;t)=\frac{1}{\pi \hslash }{\int }_{-\infty }^{+\infty }dyexp[\frac{i2qy}{\hslash }]\langle x-y|\rho (t)|x+y\rangle \hfill \end{array}$$

(36)

The initial non-equilibrium Wigner function at $t=0$ is $W_{in}$, given by Equation (36) if $\rho ={\rho }_{in}$. For $t>0$, the exact dissipationless quantum master equation for W [48,49] is:

$$\begin{array}{c}\frac{\partial W(x,q;t)}{\partial t}=-\frac{q}{m}\frac{\partial W(x,q;t)}{\partial x}+M_{Q}W\hfill \end{array}$$

(37)

$$\begin{array}{c}{M}_{Q}W={\int }_{-\infty }^{+\infty }d{q}^{\prime }W(x,q^{\prime };t){\int }_{-\infty }^{+\infty }\frac{idy}{\pi {\hslash }^2}[V(x+y)-V(x-y)]\hfill \\ \times exp[\frac{i2(q-q^{\prime })y}{\hslash }]=\frac{dV}{dx}\frac{\partial W}{\partial q}-\frac{{\hslash }^2}{3!2^2}\frac{d^{3}V}{d{x}^3}\frac{{\partial }^{3}W}{\partial q^3}+\dots \hfill \end{array}$$

(38)

As $\hslash \to 0$, Equation (37) becomes formally, by dropping all ℏ-dependent terms (containing ${\partial }^{n}W/\partial q^n$, $n=3,5,\dots$) in Equation (38), the classical Liouville Equation (1), with $W\to W_c$ [48,49]. All terms in the series in Equation (38) contribute, in general, for the assumed V’s.

Assumption. As $\mid q\mid \to +\infty$, $W(x,q;t)\to 0$ quickly, for fixed x and t, so that ${\int }_{-\infty }^{+\infty }dqW(x,q;t)q^n$ converges, for any integer $n\ge 0$ [32,34].

Under the latter assumption, Equation (37) readily implies that $(\partial /\partial t){\int }_{-\infty }^{+\infty }dx{\int }_{-\infty }^{+\infty }dqW(x,q;t)=0$.

See [19] for theorems and constructive methods to find stationary solutions of Equation (35).

We shall now consider the equilibrium Wigner function $W_{eq}$, fulfilling Equations (39) and accounting for thermal equilibrium at T with the $hb$. Like in the classical case, the solutions of Equations (37) and (38) are not expected to approach $W_{eq}$ exactly, unless some approximation be made. $W_{eq}$ arises from the canonical (t-independent) density operator ${\rho }_{eq}=exp[-\beta H]$, through Equation (36). $W_{eq}(x,q)$ is neither Gaussian in q nor known in closed form in general [48,49] and their dependences on q and x do not factorize. One has:

$$\begin{array}{c}-\frac{q}{m}\frac{\partial W_{eq}}{\partial x}+M_Q{W}_{eq}=0,\frac{\partial W_{eq}}{\partial t}=0\hfill \end{array}$$

(39)

4.2. Novel Features: $W_{eq}$ as Quasi-Definite Functional in Momentum and Orthogonal Polynomials

We now remind the known quantum difficulty: neither W nor $W_{eq}$ can be warranted to be $\ge 0$ (negativity), in general [48,49]. A necessary and sufficient condition for the Wigner function associated with a Schrodinger wave function be $\ge 0$ is that the latter be the exponential of a quadratic polynomial [53]. The domain in which $W_{eq}<0$ may occur has be consistent with ${\int }_{-\infty }^{+\infty }dx{W}_{eq}\ge 0$ and ${\int }_{-\infty }^{+\infty }dq{W}_{eq}\ge 0$.

Despite $W<0$, we shall now invoke, in outline, an acceptable mathematical framework based upon the theory of orthogonal polynomials [54]. Let us consider a kernel $K=K(y)$ (which could be $\le 0$), a set of functions $f=f(y)$ and the following functional $L_K$ determined by the kernel K: $f\to L_K[f]={\int }_{-\infty }^{+\infty }dyK(y)f(y)$. We shall suppose that all integrals over y are convergent. Let us consider, successively: ${\mu }_n=L_K[y^n]$, $n=0,1,2,3,\dots$, the set of all $(S+1)\times (S+1)$ matrices $M_S$ ($S=0,1,2,3,\dots$) with $(i,j)$-th element equal to ${\mu }_{i+j}$ ($i,j=0,\dots ,S$), and their determinants: $Det[M_S]$. By definition, the functional $L_K$ is quasi-definite if $Det[M_S]\ne 0$ for any $S=0,1,2,3,\dots$ [54]. If $L_K$ is a quasi-definite functional, then a theorem [54] implies the existence of a family of orthogonal polynomials, named here and below $H_{Q,n}=H_{Q,n}(y)$, with weight function K (even if $K<0$ in some domain in y).

Assumption—Let $y=q/q_{eq}$ ($q_{eq}={(2m/\beta )}^{1/2}$). Recall that we have assumed that V and all its derivatives exist and are continuous. Then (by extending directly [32]), one can accept that regarding their y-dependences, the Wigner function W and $W_{eq}$ determine, respectively, quasi-definite functionals $L_W$ (for any x and t) and $L_{W_{eq}}$ (for any x). We shall suppose the validity of the latter property in all that follows. Its interest is obvious: if it holds (as we suppose), it implies the existence of orthogonal polynomials. To prove that assumption for any W and $W_{eq}$, in general, lies outside our scope: its validity for $L_{W_{eq}}$ has been checked (and confirmed) in one interesting case in [32].

We shall introduce the (unnormalized) polynomials in y ($=q/q_{eq}$), $H_{Q,n}=H_{Q,n}(y)$ ($n=0,1,2,3,\dots$), orthogonalized in y (for fixed x) by using $W_{eq}$ as weight function. By choosing $H_{Q,0}(q)=1$, for $n\ne n^{\prime }$ and any x (left unintegrated), we impose:

$$\begin{array}{c}{\int }_{-\infty }^{+\infty }dy{W}_{eq}{H}_{Q,n}(y)H_{Q,n^{\prime }}(y)=0\hfill \end{array}$$

(40)

The $H_{Q,n}$’s, depending parametrically on x for $n\ge 1$, will be used for the time evolution [32]. The normalized polynomials are $H_{Q,n}/h_{Q,n}^{1/2}$ ($h_{Q,n}={\int }_{-\infty }^{+\infty }dy{W}_{eq}{H}_{Q,n}{(y)}^2$). We shall look for the $H_{Q,n}(y)$’s as:

$$\begin{array}{c}{H}_{Q,n}(y)=y^n+\sum _{j=1}^n{ϵ}_{Q,n,n-j}{y}^{n-j}\hfill \end{array}$$

(41)

${ϵ}_{Q,n,n-j}$, y-independent (but x-dependent, in general). One has ${ϵ}_{Q,n,n-j}=0$ for odd j, so that $H_{Q,n}(-y)={(-1)}^n{H}_{Q,n}(y)$. The simplest ${ϵ}_{Q,n,n-j}=0$ is: ${ϵ}_{Q,2,0}=-[{\int }_{-\infty }^{+\infty }dy{W}_{eq}{y}^2]/[{\int }_{-\infty }^{+\infty }dy{W}_{eq}]$. We omit further details of ${ϵ}_{Q,n,n-j}$: see Subsection 4.4 and Appendix A in [32].

4.3. Non-Equilibrium Moments and Hierarchy

We shall analyze general off-equilibrium situations by using Equations (37) and (38). Using for convenience the unnormalized $H_{Q,n}(y)$’s in Equation (41) ($n=0,1,2,\dots$), we introduce the new non-equilibrium moments:

$$\begin{array}{c}{W}_n=W_n(x;t)={\int }_{-\infty }^{+\infty }dy{H}_{Q,n}(y)W\hfill \end{array}$$

(42)

The initial condition $W_{in,n}$ for $W_n$ is obtained by replacing W by $W_{in}$ in Equation (42). One has the following (formal) expansion for W, which generalizes Equation (4):

$$\begin{array}{c}W=W_{eq}(x,q)\sum _{n=0}^{+\infty }{W}_n(x;t)\frac{H_{Q,n}(y)}{h_{Q,n}}\hfill \end{array}$$

(43)

For $W=W_{eq}(x,q)$, Equation (42) yields $W_{eq,n}=0$ if $n>0$, and $W_{eq,0}=h_{Q,0}$. The transformation of Equations (37) and (38) into a linear hierarchy for the new moments $W_n$ can be carried out as in the classical case, but it requires far more work. The structure of the general non-equilibrium hierarchy is:

$$\begin{array}{c}\frac{\partial W_n}{\partial t}=-M_{n,n+1}{W}_{n+1}-\sum _{n^{\prime }=1}^n{M}_{n,n-n^{\prime }}{W}_{n^{\prime }}\hfill \end{array}$$

(44)

$$\begin{array}{c}{M}_{n,n+1}{W}_{n+1}\equiv \frac{q_{eq}}{m}\frac{\partial W_{n+1}}{\partial x}\hfill \end{array}$$

(45)

$$\begin{array}{c}{M}_{n,n-1}{W}_{n-1}=\frac{q_{eq}}{m}\left[({ϵ}_{Q,n+1,n-1}-{ϵ}_{Q,n,n-2})\frac{\partial W_{n-1}}{\partial x}+(\frac{\partial {ϵ}_{Q,n,n-2}}{\partial x})W_{n-1}\right]+\frac{n}{q_{eq}}\frac{\partial V}{\partial x}{W}_{n-1}.\hfill \end{array}$$

(46)

$M_{n,n^{\prime }=0}=0$ for any n, except for $n=1$ (with $n^{\prime }=0$). In [32,34], we have obtained the first five equations in that quantum hierarchy. For any $V\ne 0$, the equation for $\partial W_n/\partial t$ for any n contains a contribution from $W_{n+1}$ (always with the same structure, with n-independent coefficients) given in Equation (45), but does not on higher order moments $W_{n^{\prime }}$, with $n^{\prime }>n+1$. On the other hand, the contribution from $W_{n^{\prime }}$ ($0<n^{\prime }\le n-1$) does carry n-dependent coefficients, which increase with n. The quantum hierarchy is not a three-term one, due to moments $W_{n^{\prime }}$ with $n^{\prime }<n$: the lowest order at which this can be displayed occurs for $n=4$, where the full quantum equation does contain a term of order ${\hslash }^3$ in $W_1$. Equations (44)–(46) constitute another result. No approximation has been carried out, so far, in this section for brevity, as they will be performed in the following sections. See [32,34] for details and further analysis and approximations in the one-dimensional case (in particular, for the long term).

5. $\mathit{N}$(>1) Quantum Particles: Equilibrium and Non-Equilibrium Statistical Distributions

5.1. General Aspects and Factoring Out the Center of Mass

We shall now consider N non-relativistic spinless quantum particles of masses $m_j$ (with $j=1,\dots ,N$) in three spatial dimensions, inside a large finite volume $\mathsf{\Omega }$. Those individual particles can be either atoms or small molecules. In cartesian coordinates, the position vector and the momentum operator of the j-th particle are $x_j=(x_{j,1},x_{j,2},x_{j,3})$ and $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$. We shall denote $(\mathbf{x})=(x_1,\dots ,x_N)$.

Assumptions—We emphasize that $\mathsf{\Omega }$ (outside which, an infinitely repulsive potential is assumed) is large only at the microscopic scale. We assume that there are no external potentials and that all interactions acting on the particles are described, in principle, by time-independent and velocity-independent instantaneous potentials V among the particles. $V=V((\mathbf{x}))$ is the sum of two-body potentials $v_{ij}=v_{ij}(|x_j-x_i|)$, as in Section 3. We assume that $v_{ij}(r)$ and all $d^n{v}_{ij}(r)/d{r}^n$, for $n=1,2,3,\dots$, are continuous for all $r>0$: compare with Section 4. The total (Hermitian) quantum Hamiltonian is

$$\begin{array}{c}H=-\sum _{j=1}^N\sum _{\alpha =1}^3\frac{{\hslash }^2}{2m_j}\frac{{\partial }^2}{\partial x_{j,\alpha }^2}+V((\mathbf{x})).\hfill \end{array}$$

(47)

Let the N-particle system be immersed in a $hb$ at thermal equilibrium, at absolute temperature T.

We resort to the non-equilibrium Wigner function [2,8,48,49,50,52] for the N particles, which reads (in terms of the N-particle density operator $\rho (t)$):

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

(48)

where $(\mathbf{q})=(q_1,\dots ,q_N)$ and $({\mathbf{dx}}^{\prime })={\prod }_{i=1}^N{\prod }_{\alpha =1}^{3}d{x}_{i,\alpha }^{\prime }$.

The time evolution of the Wigner function is furnished by

$$\begin{array}{c}\hfill \frac{\partial W((\mathbf{x}),(\mathbf{q});t)}{\partial t}=-\sum _{j=1}^N\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,\end{array}$$

(49)

with

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

(50)

where $({\mathbf{dq}}^{\prime })={\prod }_{i=1}^N{\prod }_{\alpha =1}^{3}d{q}_{i,\alpha }^{\prime }$.

Approximation—As $\mathsf{\Omega }$ is large, we shall approximate spatial integrals by those for an infinite volume when such approximations are harmless, unless some specific discussion be required. Strictly speaking, as $\mathsf{\Omega }$ is not infinite, the $({\mathbf{q}}^{\prime })$ are discretized momenta and $\int ({\mathbf{dq}}^{\prime })$ in Equation (50) should be interpreted as a $3N$-fold series. However, as $\mathsf{\Omega }$ is large, we shall disregard the small spacings in that $3N$-fold series and approximate it as a $3N$-fold integral. We shall accept that all integrals (or all series) over momenta converge for large values of the latter: explicit computations do support it.

We shall suppose a quantum initial condition, determined by some density matrix at $t=0$, ${\rho }_{in}$, eventually suitable for chemical reactions. The initial Wigner function $W_{in}((\mathbf{x}),(\mathbf{q}))$ is determined by ${\rho }_{in}$ through Equation (48). The canonical density operator describing the N-particle system at thermal equilibrium with the $hb$ is ${\rho }_{eq}=exp[-\beta H]$ (with $\partial {\rho }_{eq}/\partial t=0$). The equilibrium Wigner function $W_{eq}((\mathbf{x}),(\mathbf{q}))$ is given by Equation (48) for the actual equilibrium density operator ${\rho }_{eq}$.

Let $\mathit{P}={\sum }_{j=1}^N{\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).

For convenience, the center of mass (CM) degrees of freedom will be factored out from the relative ones from the outset. The standard CM position vector is $\mathit{X}=M^{-1}[{\sum }_{j=1}^N{m}_j{\mathit{x}}_j]$ with $M={\sum }_{j=1}^N{m}_j$. Then, one has

$$\begin{array}{c}H=H_{CM}+H_{N,re},H_{CM}=-\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial {\mathit{X}}^2},\left[H_{CM},H_{N,re}\right]=0,\hfill \end{array}$$

(51)

Here, $H_{CM}$ is the free CM Hamiltonian and $H_{N,re}$ is the Hamiltonian for the remaining (internal) relative position vectors: see Section 6 and Section 7 for $N=2$ and $N=3$, respectively. ${\rho }_{eq}$ factorizes as ${\rho }_{eq}={\rho }_{CM,eq}\otimes {\rho }_{N,re,eq}$, where ${\rho }_{CM,eq}=exp[-\beta H_{CM}]$ and ${\rho }_{N,re,eq}=exp[-\beta H_{N,re}]$. The actual $W_{eq}=W_{eq}((\mathbf{x}),(\mathbf{q}))$ factorizes as

$$\begin{array}{c}\hfill W_{eq}((\mathbf{x}),(\mathbf{q}))=W_{CM,eq}{W}_{N,re,eq},\end{array}$$

(52)

$W_{N,re,eq}$ being the relative equilibrium Wigner function. Spatial integrations over X is carried out inside the volume $\mathsf{\Omega }$. $W_{CM,eq}$ is proportional to $exp\left[-\beta ({\mathit{Q}}^2/2M)\right]$ ($\mathit{Q}={\sum }_{j=1}^N{q}_j$). Correspondingly, by using the assumptions below on the initial conditions, there is also factorization off-equilibrium: $W((\mathbf{x}),(\mathbf{q});t)=W_{CM}(\mathit{X},\mathit{Q};t)W_{N,re}$, with $\partial W_{CM}/\partial t=-(\mathit{Q}/M)\partial W_{CM}/\partial \mathit{X}$, and the (t-reversible) Wigner function $W_{N,re}$ of the relative degrees of freedom. Both $W_{N,re,eq}$ and $W_{N,re}$ depend on $N-1$ relative position vectors and momenta, denoted collectively as ${(\mathbf{x})}_{N-1,re}$ and ${(\mathbf{q})}_{N-1,re}$. In general, for given N, there are several interesting choices for the sets ${(\mathbf{x})}_{N-1,re}$ and ${(\mathbf{q})}_{N-1,re}$ out of the $(\mathbf{x})$’s and $(\mathbf{q})$’s, respectively: see, in particular, Section 7.

Assumptions—It is a reasonable assumption that the initial conditions ${\rho }_{in}$ and $W_{in}((\mathbf{x}),(\mathbf{q}))$ factorize into CM and relative ones. The initial condition for $W_{N,re}$ is $W_{N,re,in}$. Neither $W_{eq}$ nor $W_{N,re,eq}$ can be warranted to be non-negative in general [48,49,52]. $W_{\mathrm{CM},\mathrm{eq}}$, Gaussian in momenta, determines an infinite family of orthogonal polynomials (the standard Hermite ones) in momenta. We shall assume that regarding its ${(\mathit{q})}_{N-1,re}$-dependence, $W_{N,re,eq}$ determines quasi-definite functionals (for any ${(\mathit{x})}_{N-1,re}$), by extending Section 4. This assumption implies the existence of an infinite family of orthogonal polynomials generated by $W_{N,re,eq}$.

5.2. Orthogonal Polynomials, Non-Equilibrium Moments and Hierarchy

We shall consider one choice ${(\mathbf{q})}_{N-1,re}$, and let ${\mathbf{q}}_{N-1,re,j}$, $j=1,\dots ,N-1$ be the momenta constituting the set ${(\mathbf{q})}_{N-1,re}$. Also, let $m_{re,j}$ be some suitable masses, associated to the choice ${(\mathbf{q})}_{N-1,re}$ and defined out of all $m_1,\dots ,m_N$. See Section 6 and Section 7. Let $q_{eq,j}={(2m_{re,j}/\beta )}^{1/2}$ (independent of $(\mathbf{x})$) and let ${(y)}_{N-1,re}$ be the set formed by all $q_{eq,j}^{-1}{q}_{N-1,re,j}$, $j=1,\dots ,N-1$.

We consider the equilibrium distribution $W_{N,re,eq}({(\mathbf{x})}_{N-1,re},{(\mathbf{q})}_{N-1,re})$ and the orthogonal polynomials $H_{([n])}({(\mathbf{y})}_{N-1,re}))$ (depending parametrically on ${(\mathbf{x})}_{N-1,re}$), generated by the former as weight function. Here, $([n])$ denotes $N-1$ (in general, different) triplets of non-negative integers $[n]=(n_1,n_2,n_3)$: the j-th triplet $[n]$ is associated with the j-th momentum in the set ${(\mathit{q})}_{\mathrm{N}-1,\mathrm{re}}$. Use will also be made of additional obvious short-hand notations: $([0])=([0],[0],\dots .)$, ${([1])}_{j,\alpha }$ with 1 in the $\alpha$-th (=1,2,3) position of the j-th triplet with all other components of the remaining triplets being zero and ${([2])}_{j,\alpha }$ similar to ${([1])}_{j,\alpha }$ but with a 2 in the $\alpha$-th position of the j-th triplet.

The orthogonal polynomials $H_{([n])}$ are chosen with $H_{([0])}({(\mathit{y})}_{N-1,re})=1$ and:

$$\begin{array}{c}\hfill \int {(d\mathit{y})}_{N-1,re}{W}_{N,re,eq}({(\mathbf{x})}_{N-1,re},{(\mathbf{q})}_{N-1,re})H_{([n])}({(\mathit{y})}_{N-1,re})H_{([n^{\prime }])}({(y)}_{\mathit{N}-1,re})=0,\end{array}$$

(53)

${(\mathbf{dy})}_{N-1,re}={\prod }_{i=1}^{N-1}{\prod }_{\alpha =1}^{3}d{y}_{i,\alpha }$, $y_{i,\alpha }$ being the $\alpha$-th (=1,2,3) component of the i-th $q_{eq,j}^{-1}{\mathit{q}}_{N-1,re,i}$. $H_{{([1])}_{j,\alpha }}=H_{{([1])}_{j,\alpha }}(y_{j,\alpha })=y_{j,\alpha }$ and $H_{{([2])}_{j,\alpha ,l,\gamma }}=y_{j,\alpha }{y}_{l,\gamma }+{ϵ}_{{([2])}_{j,\alpha ,l,\gamma }}$, ${ϵ}_{{([2])}_{j,\alpha ,l,\gamma }}$ being ${(\mathit{y})}_{N-1,re}$ but dependent on ${(\mathit{x})}_{N-1,re}$: see Appendix A and Appendix B. The non-equilibrium moments of the actual non-equilibrium relative Wigner function read:

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

(54)

and so on for the initial non-equilibrium condition for the actual relative degrees of freedom, to be denoted as $W_{in,([n])}$ for simplicity. The transformation of Equations (49) and (50) (after having factored out the CM) and Equation (54) into an infinite linear hierarchy for the non-equilibrium moments $W_{([n])}$ can be carried out as in previous cases. The structure of the general (t-reversible) linear non-equilibrium hierarchy, a genuine consequence of quantum mechanics, reads:

$$\begin{array}{c}\hfill \frac{\partial W_{([n])}}{\partial t}=-\sum _{j=1}^{N-1}\sum _{\gamma =1}^3{M}_{([n]),{([n+1])}_{j,\gamma }}{W}_{{([n+1])}_{j,\gamma }}-\sum _{([n^{\prime }])}{M}_{([n]),([n-n^{\prime }])}{W}_{([n-n^{\prime }])},\end{array}$$

(55)

with

$$\begin{array}{c}\hfill M_{([n]),{([n+1])}_{j,\gamma }}{W}_{{([n+1])}_{j,\gamma }}\equiv \frac{q_{eq,j}}{m_{re,j}}\frac{\partial W_{{([n+1])}_{j,\gamma }}}{\partial x_{j,\gamma }}.\end{array}$$

(56)

Equation (55) is a result which generalizes Equation (44). The subscript ${([n+1])}_{j,\gamma }$ in $W_{{([n+1])}_{j,\gamma }}$ denotes that one unity has been added to the integer in the $(j,\gamma )$ location, and analogously for the corresponding interpretation of the sum ${\sum }_{([n^{\prime }])}$ in Equation (44): notice that the non-negative integers in $([n^{\prime }])$ grow until ${\sum }_{j=1}^{N-1}{\sum }_{\alpha =1}^3{n}_{j\alpha }^{\prime }={\sum }_{j=1}^{N-1}{\sum }_{\alpha =1}^3{n}_{j\alpha }$. Further explanations are omitted, see [34]. Use will be made of:

$$\begin{array}{c}{M}_{{([1])}_{j,\alpha },([0])}{W}_{([0])}=-\frac{q_{eq,j}}{m_{re,j}}\sum _{l=1}^{N-1}\sum _{\gamma =1}^3\frac{\partial }{\partial x_{l,\gamma }}\left({ϵ}_{{([2])}_{j,\alpha ;l,\gamma }}{W}_{([0])}\right)+\frac{1}{q_{eq,j}}\frac{\partial V}{\partial x_{j,\alpha }}{W}_{([0])}.\hfill \end{array}$$

(57)

$$\begin{array}{c}{ϵ}_{{[2]}_{j,\alpha ,l,\gamma }}=-\frac{\int {(d\mathit{y})}_{N-1,re}{W}_{N,re,eq}({(\mathbf{x})}_{N-1,re},{(\mathbf{q})}_{N-1,re})y_{j,\alpha }{y}_{l,\gamma }}{\int {(d\mathit{y})}_{N-1,re}{W}_{N,re,eq}({(\mathbf{x})}_{N-1,re},{(\mathbf{q})}_{N-1,re})},\hfill \end{array}$$

(58)

Equations (A1) and (A25) will be particular cases of Equation (58). The quantum Equation (55) is also a counterpart of the classical Equation (26). We emphasize three crucial difference between them. First, in the classical Equation (26), the factoring out of $W_{c,eq,0}{(x)}^{-1/2}$ enabled to work with the symmetrized moments $g_n$, while no similar factorization is possible in Equation (55), so that one has to deal with the $W_{([n])}$’s. Second, the CM has not been factored out in the classical Equation (26), while it has indeed been in the quantum Equation (55). Third, in Equation (26) the orthogonal polynomials were independent of $(\mathbf{x})$ and, consequently, the only dependences on the two-body potentials are those displayed explicitly in the right-hand-side of Equation (26). In the quantum Equation (55), the orthogonal polynomials are ${(\mathbf{x})}_{\mathrm{N}-1,\mathrm{re}}$-dependent, and so are ${ϵ}_{{([2])}_{j,\alpha },l,\gamma }$ and, in general, the coefficients of $H_{([n])}$ contained in $M_{([n]),([n-n^{\prime }])}$. All coefficients in (55) are expressed in terms of V and of quantities computed out of the equilibrium solution $W_{N,re,eq}$. When $W_{N,re}=W_{N,re,eq}$, all equilibrium moments $W_{eq,([n])}$ with $([n])\ne ([0])$ vanish, except $W_{eq,([0])}$. Moreover, $W_{eq,([0])}(\ne 0)$ and $W_{eq,([n])}=0$ for $([n])\ne ([0])$ solve the hierarchy Equation (55), with all $\partial W_{eq,([n])}/\partial t=0$.

5.3. Approximations: Small Thermal Wavelength and Long-Term

In the t-evolution, not far from thermal equilibrium, and recalling [40,41,42,43,44,45,46,47], one expects that $W_{([0])}$ would be dominant and that any $W_{([n])}$, for any $([n])\ne ([0])$, be small (the more negligible the larger t and ${\sum }_{j=1}^{N-1}{\sum }_{\alpha =1}^3{n}_{j\alpha }$ are).

Assumption and approximation—We consider that T is not high (so as not to enter fully in the classical regime). We shall also avoid the very low temperature regime. Therefore, we shall treat the regime of typical chemical reactions: small thermal wavelength. Then, the general non-equilibrium hierarchy Equation (55) can be approximated by a three-term hierarchy, namely:

$$\begin{array}{c}\hfill \frac{\partial W_{([n])}}{\partial t}=-\sum _{j=1}^{N-1}\sum _{\gamma =1}^3{M}_{([n]),{([n+1])}_{j,\gamma }}{W}_{{([n+1])}_{j,\gamma }}-\sum _{}{M}_{([n]),([n-1])}{W}_{([n-1])},\end{array}$$

(59)

The short-hand notation ${\sum }_{}$ in Equation (59) means that only $[n^{\prime }]=[1]$ contributes (compared to ${\sum }_{([n^{\prime }])}$ in Equation (55), where various $[n^{\prime }]$ have to be considered). That short-hand notation omits, for brevity, dependences on certain indices: compare with the more explicit Equation (67). The justification of the above small thermal wavelength, leading from Equations (55)–(59) will be omitted, as it is a straightforward, but cumbersome, generalization of the one to be given in Section 6.2 for $N=2$. We emphasize that Equation (59) is still a quantum hierarchy, since it contains the quantum-mechanical $M_{([n]),{([n-1])}_{j,\gamma }}$.

The solution of the resulting approximate non-equilibrium three-term hierarchy Equation (59) 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[([n]);s]$, which constitute extensions of Equations (8) and (29).

Assumption and approximation—We proceed to the long-term approximation (thus introducing irreversibility) for t longer than a certain largest effective evolution time. We shall omit details, which constitute a formal extension of those in Section 2.5. We shall assume the initial condition $W_{in,([0])}\ne 0$ and $W_{in,([n])}=0$ for $([n])\ne ([0])$.

Upon taking the inverse Laplace transform, one gets (for fixed and small $s=ϵ>0$) the following irreversible approximate quantum (Smoluchowski-like) equation for the lowest non-equilibrium moment $W_{([0])}(x;t)$

$$\begin{array}{c}\hfill \frac{\partial W_{([0])}}{\partial t}=\sum _{j=1}^{N-1}\sum _{\alpha =1}^3\frac{q_{eq,j}}{m_{re,j}}\frac{\partial }{\partial x_{j,\alpha }}D[M_{{([1])}_{j,\alpha },([0])}{W}_{([0])}],\end{array}$$

(60)

with D being a suitable approximation for the corresponding combination of quantum operator-continued fractions (omitting various dependences and indices in it) [32,34] and assuming approximately the initial condition $W_{in,([0])}$. This procedure and approximations are, in principle, independent of the specific properties of ${ϵ}_{{([2])}_{j,\alpha ,l,\gamma }}$ (which displays quantum effects: compare with Section 6 and Section 7). The result in Equation (60) 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.

6. $\mathit{N}$ = 2 Quantum Particles: Chemical Reactions

We shall apply the general formalism in Section 5 for N=2.

6.1. Assumptions on Two-Body Potential and Non-Equilibrium Moments and Hierarchy

Here: $V=v(|x_1-x_2|)$. Let us introduce the relative position vector $x=x_2-x_1$, the reduced mass m through $m^{-1}=m_1^{-1}+m_2^{-1}$ and $q_{eq}={(2m/\beta )}^{1/2}$.

$$\begin{array}{c}{H}_{2,re}=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}+v(|x|).\hfill \end{array}$$

(61)

Spatial integrations over x is carried out inside the volume $\mathsf{\Omega }$.

Assumptions—The eigenfunctions corresponding to both $H_{CM}$ and $H_{2,re}$ will vanish, by assumption, at the surfaces enclosing $\mathsf{\Omega }$. We suppose that $v(r)$ is an effective potential between two atoms or small molecules (which, in particular, includes and averages over Coulomb interactions). We shall suppose that such generic two-body potential $v(r)$ satisfies the following conditions:

• $v(r)$ 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(r)$ 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(r)$ and all $d^{n}v(r)/d{r}^n$, for $n=1,2,3,\dots$, are continuous for all $r>0$. Recall Section 4.1 and Section 5.1.
• $v(r)$ 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.

Let ${\phi }_j={\phi }_j(x)$ denote a suitably normalized eigenfunction of $H_{2,re}$ with corresponding eigenvalue $E_j$ and j denoting a set of labels. $H_{2,re}$ has both a discrete spectrum, corresponding to just one bound state and, above it and disjoint from it, to a denumerably infinite number of discrete states. For the bound state $E_j=E_d<0$ and ${\phi }_j={\phi }_d(x)$. The denumerably infinite discrete spectrum (above the bound state) 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 continuous spectrum (CS), both if either the small ${\mathsf{\Omega }}^{-1}$ remains positive or if ${\mathsf{\Omega }}^{-1}\to 0$. The eigenfunctions corresponding to the CS of $H_{2,re}$ are ${\phi }_j={\phi }_{\mathit{k}}(\mathit{x})$, with $j\equiv \mathit{k}$ being an almost continuous wavevector, and eigenvalues $E_j=E_{\mathit{k}}\ge 0$. The CS eigenfunctions are normalized, if ${\mathsf{\Omega }}^{-1}$ remains positive, through: $({\phi }_{\mathit{k}},{\phi }_{{\mathit{k}}^{\prime }})=\int d^{3}x{\phi }_{\mathit{k}}^{*}{\phi }_{{\mathit{k}}^{\prime }}={\delta }_{\mathit{k}{\mathit{k}}^{\prime }}$ (a Kronecker delta). Also, $({\phi }_d,{\phi }_d)=\int d^3\mathit{x}{\phi }_d^{*}{\phi }_d=1$ (normalized) and $({\phi }_d,{\phi }_{\mathit{k}})=\int d^3\mathit{x}{\phi }_d^{*}{\phi }_{\mathit{k}}=0$. Hence, ${\phi }_d$ and all CS ${\phi }_{\mathit{k}}$ span two separate Hilbert subspaces ${\mathcal{H}}_d$ and ${\mathcal{H}}_{CS}$. Let ${\sum }_j$ be a short-hand notation for summation (and, eventually, integration) which includes the contribution of both the bound state eigenfunction plus, if ${\mathsf{\Omega }}^{-1}$ remains positive, that of a three-fold infinite summation over the whole CS. For the CS only: ${\sum }_j\to (\mathsf{\Omega }/{(2\pi )}^3)\int d^3\mathit{k}$ as ${\mathsf{\Omega }}^{-1}\to 0$. See [34]. Let $\mathit{q}=(-m_2{\mathit{q}}_1+m_1{\mathit{q}}_2)/(m_1+m_2)$ be the relative momentum vectors. The equilibrium Wigner function for the relative degrees of freedom reads:

$$\begin{array}{c}{W}_{2,re,eq}\left(\mathit{x},\mathit{q}\right)=\frac{1}{{(\pi \hslash )}^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}$$

(62)

We shall consider the infinite family of orthogonal polynomials $H_{[n]}(\mathit{y})$ generated by $W_{2,re,eq}$, as indicated in Section 5.2. We define $[n]=(n_1,n_2,n_3)$, with $[0]=(0,0,0)$, ${[1]}_{\alpha }$ representing a vector with a 1 in the $\alpha$ position with the rest of the components being zero, and ${[2]}_{\alpha }$ similar to ${[1]}_{\alpha }$ but with a 2 in the $\alpha$ position. We shall also employ: $H_{{[1]}_{\alpha }}=H_{{[1]}_{\alpha }}(\mathit{y})=y_{\alpha }$, and $H_{{[2]}_{\alpha ,\gamma }}=y_{\alpha }{y}_{\gamma }+{ϵ}_{{[2]}_{\alpha ,\gamma }}$. ${ϵ}_{{[2]}_{\alpha ,\gamma }}$ is given in Appendix A. The derivation of an infinite linear hierarchy for the non-equilibrium moments $W_{[n]}$ ($=\int d^3\mathit{y}{H}_{[n]}(\mathit{y})W_{2,re}\left(\mathit{x},\mathit{y};t\right)$) is a particular case of that in Section 5.2. The initial condition $W_{in,[n]}$ for $W_{[n]}$ is obtained by replacing $W_{2,re}$ by $W_{2,re,in}$. The general (t-reversible) hierarchy for $N=2$ and any $[n]$ is ($m_{re,j}=m$, $q_{eq,j}=q_{eq}$), as a result:

$$\begin{array}{c}\frac{\partial W_{[n]}}{\partial t}=-\sum _{\gamma =1}^3{M}_{[n],[{(n+1)}_{\gamma }]}{W}_{[{(n+1)}_{\gamma }]}-\sum _{[n^{\prime }]}{M}_{[n],[n-n^{\prime }]}{W}_{[n-n^{\prime }]},\hfill \end{array}$$

(63)

with

$$\begin{array}{c}{M}_{[n],[{(n+1)}_{\gamma }]}{W}_{{[n+1]}_{\gamma }}\equiv \frac{q_{eq}}{m}\frac{\partial W_{{[n+1]}_{\gamma }}}{\partial x_{\gamma }}.\hfill \end{array}$$

(64)

$$\begin{array}{c}{M}_{{[1]}_{\alpha },[0]}{W}_{[0]}=-\frac{q_{eq}}{m}\sum _{\gamma =1}^3\frac{\partial }{\partial x_{\gamma }}\left({ϵ}_{{[2]}_{\alpha ,\gamma }}{W}_{[0]}\right)+\frac{1}{q_{eq}}\frac{\partial v}{\partial x_{\alpha }}{W}_{[0]}.\hfill \end{array}$$

(65)

Bearing in mind that $W_{[n]}=W_{n_1,n_2,n_3}$, the subscript ${[n+1]}_{\gamma }$ in $W_{[{(n+1)}_{\gamma }]}$ 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 so on for the corresponding interpretation of the sum ${\sum }_{[n^{\prime }]}{M}_{[n],[n-n^{\prime }]}{W}_{[n-n^{\prime }]}$ in (63), which is performed as the non-negative integers in $[n^{\prime }]$ grow until ${\sum }_{\alpha =1}^3{n}_{\alpha }^{\prime }={\sum }_{\alpha =1}^3{n}_{\alpha }$. We shall omit, for simplicity, the $M_{[n],[n-n^{\prime }]}$’s. By invoking rotational invariance, the operator $M_{{[1]}_{\alpha },[0]}$ can be expected to behave as a three-dimensional vector regarding the index $\alpha$ ($=1,2,3$). Equation (65) corrects Equation (23) in [34], by adding what may well be called an off-diagonal correction (namely ${\sum }_{\gamma =1,\gamma \ne \alpha }^3$). Equation (A2) for $\alpha =\gamma$ versus $\alpha \ne \gamma$ and the fact that $v=v(r)$ illustrate that due to their angular dependences (and upon eventual averaging or integrations over the angles involved in $x_{\alpha }$ and $x_{\gamma }$), the contributions of the off-diagonal correction are overcome, in practice, by those due to ${ϵ}_{{[2]}_{\alpha ,\alpha }}$. Such off-diagonal correction does not alter the developments and conclusions in [34], which relied on the diagonal contribution, namely that with $\gamma =\alpha$.

Let $W_{2,re}=W_{2,re,eq}$. By using the expression for ${ϵ}_{{[2]}_{\alpha ,\gamma }}$, given in Appendix A, together with Equation (57), we obtain that

$$\begin{array}{c}\hfill M_{{[1]}_{\alpha },[0]}{W}_{eq,[0]}=0.\end{array}$$

(66)

6.2. Chemical Reactions: Assumptions, Order of Magnitude Estimates and Approximations

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 Å). Let $m=n\times m_{ne}$, with $m_{ne}$ being the neutron mass. 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 1 and 10 electronvolts. Notice that $v_0$ has a similar order of magnitude as typical energies $|E_d|$ for the bound state (see [5]). $\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,\dots$, 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 }_{th}=\hslash {[\beta /2m]}^{1/2}=\hslash /q_{eq}$ be some suitable thermal wavelength. We shall assume that n and T are such that the following conditions are fulfilled: (a) ${\lambda }_{th}$ is smaller than $\delta x$, which sets a lower limit on $nT$; and (b) $(3/2)k_{B}T$ is smaller than or, at most, about $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 }_{th}$). Thus, even if the relative particle is in the quantum regime (but not in the classical high-temperature one), ${\lambda }_{th}$ is, on average, smaller than $\delta x$ and a.

As zeroth-order approximations (see [34]), we shall take ${ϵ}_{{[2]}_{\alpha ,\alpha }}\simeq -(1/2)$ for $\left|x\right|>a+r_0$ (dominated by the CS and, in turn, approximated through the classical Boltzmann equilibrium distribution and ${ϵ}_{{[2]}_{\alpha ,\alpha }}\simeq -({\lambda }_{th}^2/{(\delta x)}^2)$ for $r_0<\left|x\right|<a+r_0$ (dominated by the bound state), both independent of the index $\alpha$ and on the angular coordinates. The order of magnitude of ${ϵ}_{{[2]}_{\alpha ,\gamma }}$, $\alpha \ne \gamma$ can be estimated to be smaller than those of ${ϵ}_{{[2]}_{\alpha ,\alpha }}$.

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

The orders of magnitude of $M_{[n],[n-1]}$ can be expected to be respectively similar to the various contributions of $M_{{[1]}_{\alpha },[0]}$. On the other hand, the specific quantum contributions $M_{[n],[n-n^{\prime }]}$, $n^{\prime }>1$ have structures proportional to $({\hslash }^{2n^{″}}/q_{eq}^{2n^{″}+1})$ times spatial partial derivatives of V of order $2n^{″}+1$ (for various $n^{″}$, $2n^{″}+1$ growing as $2^{-1}{\sum }_{\alpha =1}^3{n}_{\alpha }-1$). Such structures yield contributions smaller than that of $(1/q_{eq})(\partial V/\partial x_{\alpha })$ by factors ${({\lambda }_{th}/\delta x)}^{2n^{″}}\times ({\delta }^{2n^{″}}V/\delta V)$, by virtue of assumption a) above. Moreover, one could expect $|({\delta }^{2n^{″}}v/\delta v)|$ to be smaller than unity and to decrease as $n^{″}$ grows. The quantum contributions $M_{[n],[n-n^{\prime }]}$, $n^{\prime }>1$, will be neglected from now on, compared to that of $M_{[n],[n-1]}$. Under the above assumption a), the overall order of magnitude involved upon applying $M_{[n],[n-1]}$ in (44) is about the same as for $M_{{[1]}_{\alpha },[0]}$. Namely, the sum of (i) plus (ii) plus (iii). For further details, see [34].

6.3. Approximate Equation for Lowest Non-Equilibrium Moment

By using the above assumptions and estimates, one analyzes the non-equilibrium hierarchy Equation (63) for the $W_{[n]}$’s, and proceeds to the small thermal wavelength quantum regime. In the latter, x varies by units of order $\delta x>{\lambda }_{th}$. Then, by neglecting all $M_{[n],{[n-n^{\prime }]}_{\gamma }}$’s with ${\sum }_{\alpha =1}^3{n}_{\alpha }^{\prime }>1$, Equation (63) becomes, as a result, the approximate (t-reversible) three-term hierarchy (actually, the $N=2$ counterpart of Equation (59)):

$$\begin{array}{c}\hfill \frac{\partial W_{[n]}}{\partial t}=-\sum _{\gamma =1}^3{M}_{[n],[{(n+1)}_{\gamma }]}{W}_{[{(n+1)}_{\gamma }]}-\sum _{\gamma =1}^3{M}_{[n],[n-1]}{W}_{[n-1]}.\end{array}$$

(67)

Assumption and approximation—We shall assume the initial condition $W_{in,[0]}\ne 0$ and $W_{in,[n]}=0$ for $[n]\ne [0]$ We solve the approximate Equation (67) and proceed to the long-term approximation as in Section 5.3. One gets the following irreversible approximate quantum (Smoluchowski-like) equation for the lowest non-equilibrium moment $W_{[0]}(x;t)$:

$$\begin{array}{c}\hfill \frac{\partial W_{[0]}}{\partial t}=\frac{q_{eq}}{m}\sum _{\alpha =1}^3\frac{\partial }{\partial x_{\alpha }}D[M_{{[1]}_{\alpha },[0]}{W}_{[0]}],\end{array}$$

(68)

this result (with a constant D) is the $N=2$ counterpart of Equation (60). At a later stage, the constant D could be supposed to be positive. See [34] for a detailed analysis and applications to a binary chemical reaction. To the best of our knowledge, no general agreement has been reached to define a unique non-equilibrium chemical potential based on non-equilibrium statistical mechanics, so far [20]. Within such a limitation, approximate chemical equilibrium and kinetic equations for binary reactions have been analyzed in [34]. For binary chemical reactions in the quantum domain, mean first passage methods enabled to study approximately the times required in the transitions between bound state and unbound configurations and to exhibit typical Arrhenius exponential factors [34].

6.4. Comparison: Binary Chemical Reactions in Classical Statistical Mechanics

We focus on the classical statistical evolution of the relative degrees of freedom, described by the relative position vector $\mathit{x}={\mathit{x}}_2-{\mathit{x}}_1$ and momentum q and by $W_c(\mathit{x},\mathit{q};t)$ (with factored out the CM). Here we shall work with the moments $W_{c,re}([n])=W_{c,re}(\mathit{x};[n];t)$ for the relative degrees of freedom (equivalent to and more advantageous than the $g([n])$’s in Section 3 for dealing with chemical reactions).

Assumption and approximation—We shall assume the initial condition $W_{c,re,in}([0])\ne 0$ and $W_{c,re,in}([n])=0$ for $[n]\ne [0]$. The developments in Section 3 can be directly applied to the present case and yield an infinite hierarchy for all moments $W_{c,re}([n])$. The long-term approximation for the $g([n])$’s can be recast as a corresponding one for the $W_{c,re}([n])$.

We arrive, as a result, at the approximate irreversible equation for the lowest non-equilibrium moment $W_{c,re}([0])$, with initial condition $W_{c,in}([0])$:

$$\begin{array}{c}\frac{\partial W_{c,re}([0])}{\partial t}=\frac{q_{eq}}{m}\sum _{\alpha =1}^3\frac{\partial }{\partial x_{\alpha }}D[M_{{[1]}_{\alpha },[0],c}{W}_{c,re}([0])]\hfill \end{array}$$

(69)

$$\begin{array}{c}{M}_{{[1]}_{\alpha },[0],c}{W}_{c,re}([0])=\frac{q_{eq}}{2m}\frac{\partial }{\partial x_{\alpha }}\left(W_{c,re}([0])\right)+\frac{1}{q_{eq}}\frac{\partial v}{\partial x_{\alpha }}{W}_{c,re}([0]).\hfill \end{array}$$

(70)

The main differences between Equations (65) and (70) are the vanishing of ${ϵ}_{{[2]}_{\alpha ,\gamma }}$ for $\alpha \ne \gamma$ and the replacement of $-{ϵ}_{{[2]}_{\alpha ,\alpha }}$ by $1/2$, so that $M_{{[1]}_{\alpha },[0]}$ becomes $M_{{[1]}_{\alpha },[0],c}$. The equilibrium solution $W_{c,re,eq}([0])$ ($M_{{[1]}_{\alpha },[0],c}{W}_{c,re,eq}([0])=0$) is proportional to $exp(-\beta V)$. Equation (69) is a standard classical diffusion equation: so, we have rederived Equation (7.30) in [55]. The interest of the latter equation stems from the fact that for suitable v (for instance, the well-known Morse potential) and even in the absence of explicit quantum effects, it has been used as a zeroth-order model for a binary chemical reaction. Specifically, it has enabled to estimate the recombination rate of a diffusing molecule in three spatial dimensions, through mean first-time passage methods [55].

7. $\mathit{N}$ = 3 Quantum Particles: Chemical Reactions

We shall apply the general formalism in Section 5 for N = 3.

7.1. Some General Aspects

Assumptions—Here:

$$\begin{array}{c}V=v_{12}(|{\mathit{x}}_2-{\mathit{x}}_1|)+v_{23}(|{\mathit{x}}_3-{\mathit{x}}_2|)+v_{31}(|{\mathit{x}}_1-{\mathit{x}}_3|)\hfill \end{array}$$

(71)

We assume that any of the two-body potentials in Equation (71) fulfill the same assumptions as $v(r)$ in Section 6.1 (r denoting here the distance between any two particles, say, any of $|{\mathit{x}}_2-{\mathit{x}}_1|$, $|{\mathit{x}}_3-{\mathit{x}}_2|$, $|{\mathit{x}}_1-{\mathit{x}}_3|$). Then, the two particle 1 and 2 have one bound state with energy $E_{d(12)}(<0)$ and so on for 2 and 3 and for 3 and 1. We assume $E_{d(12)}<E_{d(23)}<E_{d(31)}(<0)$. For simplicity, we shall assume that there are no bound states of the three particles.

We shall also factor out the CM. We have three natural choices or options for the relative degrees of freedom, denoted as $((12),3)$, $((23),1)$ and $((31),2)$. See Appendix B for the corresponding relative position vectors and momenta and for several useful formulae employed in our study, without further reference to them. We shall consider the $((12),3)$ choice. Then, $H_{3,re}$ becomes:

$$\begin{array}{c}{H}_{3,re}=-\frac{{\hslash }^2}{2m_{12}}\frac{{\partial }^2}{\partial x_{12}^2}-\frac{{\hslash }^2}{2m_{3(12)}}\frac{{\partial }^2}{\partial x_{3(12)}^2}+V((\mathbf{x})).\hfill \end{array}$$

(72)

and so on for the representations of $H_{3,re}$ in terms of the other two choices. Having discussed in Section 6.1 the spectrum of $H_{2,re}$ for small ${\mathsf{\Omega }}^{-1}$ and for ${\mathsf{\Omega }}^{-1}\to 0$, we shall outline the spectrum of $H_{3,re}$ directly for ${\mathsf{\Omega }}^{-1}\to 0$. General studies on the quantum-mechanical three-body problem [56] imply that $H_{3,re}$ has: (i) an infinite CS (i)) in $(E_{d(12)},+\infty )$, associated with the choice $((12),3)$, corresponding to particle 3 interacting with the bound state of (1,2), so that the former could be far from the latter, other two similar CS corresponding respectively to the choices (ii) $((23),1)$ (CS (ii)) and (iii) $((31),2)$ (CS (iii)), and (iv) an infinite CS (CS (iv)) in $(0,+\infty )$ corresponding to the three unbound particles, so that any of them could be far from the other two, but still interacting among themselves. Thus, CS (ii) corresponds to particle 1 interacting with the bound state of (2, 3), so that the former could be far from the latter and so on for CS (iii). Here, ${\sum }_j$ will denote a short-hand notation for summations and integrations over the whole 4 infinite CS’s in (i), (ii), (iii) and (iv) above. Below, ${\phi }_j={\phi }_j(x_{12},x_{3(12)})$ will denote an eigenfunction of $H_{3,re}$ (suitably normalized as ${\mathsf{\Omega }}^{-1}\to 0$, in analogy with Section 6.1) with continuous eigenvalue $E_j$. The equilibrium Wigner function for the relative degrees of freedom in the $((12),3)$ choice reads:

$$\begin{array}{c}{W}_{3,re,eq}\left({\mathit{x}}_{12},{\mathit{x}}_{3(12)},{\mathit{q}}_{12},{\mathit{q}}_{3(12)}\right)=\frac{1}{{(\pi \hslash )}^6}\int d^3{\mathit{x}}_{12}^{\prime }{d}^3{\mathit{x}}_{3(12)}^{\prime }exp\left(\frac{2i({\mathit{q}}_{12}{\mathit{x}}_{12}^{\prime }+{\mathit{q}}_{3(12)}{\mathit{x}}_{3(12)}^{\prime })}{\hslash }\right)\times \hfill \\ \sum _{j}exp[-\beta E_j]{\phi }_j({\mathit{x}}_{12}-{\mathit{x}}_{12}^{\prime },{\mathit{x}}_{3(12)}-{\mathit{x}}_{3(12)}^{\prime }){\phi }_j^{*}({\mathit{x}}_{12}+{\mathit{x}}_{12}^{\prime },{\mathit{x}}_{3(12)}+{\mathit{x}}_{3(12)}^{\prime }).\hfill \end{array}$$

(73)

${\mathit{q}}_{12}$ and ${\mathit{q}}_{3(12)}$ being defined in Appendix B, and so on for the other two choices. Notice that $W_{3,re,eq}$ is the same for the three choices.

7.2. Approximate Equation for Lowest Non-Equilibrium Moment

Here, $q_{eq,1}=q_{eq,12}={(2m_{12}/\beta )}^{1/2}$ and $q_{eq,2}=q_{eq,3(12)}={(2m_{3(12)}/\beta )}^{1/2}$.

Assumptions and approximations—For $N=3$, we shall follow the general strategy in Section 5.2 and Section 5.3 and make assumptions similar to those for $N=2$ in Section 6.2. Then, the estimates in Section 6.2 can be readily extended for $N=3$, $q_{eq}$ being now replaced by $q_{eq,12}$ and $q_{eq,3(12)}$ (all of which with more or less similar orders of magnitude) and v denoting now any of the two-body potentials (no estimates being given explicitly for ${ϵ}_{{([2])}_{j,\alpha ,l,\gamma }}$, in our short analysis). We shall consider the regime of small thermal wavelength and the long-term approximation.

Then, as a result, one gets the following approximate irreversible quantum (Smoluchowski-like) equation for the lowest non-equilibrium moment $W_{([0])}({\mathit{x}}_{12}{\mathit{x}}_{3(12)};t)$ for the relative degrees of freedom in the ((12), 3) choice: (${\mathit{x}}_{12}=({\mathit{x}}_{12,\alpha })$, ${\mathit{x}}_{3(12)}=(x_{3(12),\alpha })$, $\alpha =1,2,3$):

$$\begin{array}{c}\frac{\partial W_{([0])}}{\partial t}=\frac{q_{eq,12}}{m_{12}}\sum _{\alpha =1}^3\frac{\partial }{\partial x_{12,\alpha }}{D}_{12}[M_{{([1])}_{12,\alpha },([0])}{W}_{([0])}]+\hfill \\ \frac{q_{eq,3(12)}}{m_{3(12)}}\sum _{\alpha =1}^3\frac{\partial }{\partial x_{3(12),\alpha }}{D}_{3(12)}[M_{{([1])}_{3(12),\alpha },([0])}{W}_{([0])}]\hfill \end{array}$$

(74)

$$\begin{array}{c}{M}_{{([1])}_{12,\alpha },([0])}{W}_{([0])}=-\frac{q_{eq,12}}{m_{12}}\sum _{\gamma =1}^3\frac{\partial }{\partial x_{12,\gamma }}\left({ϵ}_{{([2])}_{12,\alpha ,12,\gamma }}{W}_{([0])}\right)\hfill \\ -\frac{q_{eq,3(12)}}{m_{3(12)}}\sum _{\gamma =1}^3\frac{\partial }{\partial x_{3(12),\gamma }}\left({ϵ}_{{([2])}_{12,\alpha ,3(12),\gamma }}{W}_{([0])}\right)+\frac{1}{q_{eq,12}}\frac{\partial V}{\partial x_{12,\alpha }}{W}_{([0])}\hfill \end{array}$$

(75)

$$\begin{array}{c}{M}_{{([1])}_{3(12),\alpha },([0])}{W}_{([0])}=-\frac{q_{eq,12}}{m_{12}}\sum _{\gamma =1}^3\frac{\partial }{\partial x_{12,\gamma }}\left({ϵ}_{{([2])}_{3(12),\alpha ,12,\gamma }}{W}_{([0])}\right)\hfill \\ -\frac{q_{eq,3(12)}}{m_{3(12)}}\sum _{\gamma =1}^3\frac{\partial }{\partial x_{3(12),\gamma }}\left({ϵ}_{{([2])}_{3(12),\alpha ,3(12),\gamma }}{W}_{([0])}\right)+\frac{1}{q_{eq,3(12)}}\frac{\partial V}{\partial x_{3(12),\alpha }}{W}_{([0])}\hfill \end{array}$$

(76)

with the initial condition $W_{in,([0])}\ne 0$ ($W_{in,([n])}=0$ for $([n])\ne ([0])$). The $ϵ$’s are defined in Appendix B and obtained near the classical regime (to next to leading order). We have approximated the corresponding combinations of quantum operator-continued fractions by the constants $D_{12}$ and $D_{3(12)}$. This procedure and approximations are, in principle, independent of the specific properties of the $ϵ$’s. Notice that $W_{([0])}$ is the same for the three choices $((12),3)$ $((23),1)$ and $((31),2)$, by virtue of Equation (A23). On the other hand, this unique $W_{([0])}$ satisfies other two equations similar to Equation (74), for the other two choices $((23),1)$ and $((31),2)$.

Notice that for $W_{3,re}=W_{3,re,eq}$, the corresponding moment $W_{eq,([0])}$ fulfills both exact equations $M_{{([1])}_{12,\alpha },([0])}{W}_{eq,([0])}=0$ and $M_{{([1])}_{3(12),\alpha },([0])}{W}_{\mathit{eq},([0])}=0$, which generalize Equation (66). Suitable combinations of the last two exact equations for the ((12),3) choice have to transform into the corresponding exact two equations for each of the ((23), 1) and ((31), 2 choices): there are three equivalent pairs of equations for the same $W_{eq,([0])}$. We omit details.

7.3. Comparison: $N=3$ Particles in Classical Statistical Mechanics

We shall treat the $N=3$ counterpart of Equations (69) and (70), for the classical $W_{c,re}([0])$, for the relative degrees of freedom. Through the corresponding assumptions and approximations, we get as a result, for the ((12), 3) choice:

$$\begin{array}{c}\frac{\partial W_{c,re}(([0]))}{\partial t}=\frac{q_{eq,12}}{m_{12}}\sum _{\alpha =1}^3\frac{\partial }{\partial x_{12,\alpha }}{D}_{12,c}[M_{{([1])}_{12,\alpha },([0]),c}{W}_{c,re}(([0]))]+\hfill \\ \frac{q_{eq,3(12)}}{m_{3(12)}}\sum _{\alpha =1}^3\frac{\partial }{\partial x_{3(12),\alpha }}{D}_{3(12)c}[M_{{([1])}_{3(12),\alpha },([0]),c}{W}_{c,re}(([0]))]\hfill \end{array}$$

(77)

$$\begin{array}{c}{M}_{{([1])}_{12,\alpha },([0]),c}{W}_{c,re}(([0]))=\frac{q_{eq,12}}{2m_{12}}\frac{\partial }{\partial x_{12,\alpha }}\left(W_{c,re}(([0]))\right)+\frac{1}{q_{eq,12}}\frac{\partial V}{\partial x_{12,\alpha }}{W}_{c,re}(([0]))\hfill \end{array}$$

(78)

$$\begin{array}{c}{M}_{{([1])}_{3(12),\alpha },([0]),c}{W}_{c,re}(([0]))=\frac{q_{eq,3(12)}}{2m_{3(12)}}\frac{\partial }{\partial x_{3(12),\alpha }}\left(W_{c,re}(([0]))\right)+\frac{1}{q_{eq,3(12)}}\frac{\partial V}{\partial x_{3(12),\alpha }}{W}_{c,re}(([0]))\hfill \end{array}$$

(79)

These turn out to come from Equations (74)–(76) with simplifications in the classical regime, namely, $-{ϵ}_{{([2])}_{12,\alpha ,12,\alpha }}$ and $-{ϵ}_{{([2])}_{3(12),\alpha ,3(12),\alpha }}$ are approximated by 1/2, while all the remaining $-{ϵ}_{{([2])}_{j,\alpha ,l,\gamma }}$’s vanish. We have approximated the corresponding constants $D_{12}$ and $D_{3(12)}$ by the constants $D_{12,c}$ and $D_{3(12),c}$ (to be supposed positive, eventually).

Let us also consider, avoiding writing them in detail: (i) the t-dependent counterpart of Equation (77) for the ((23), 1) choice, with corresponding constants $D_{23,c}$ and $D_{1(23),c}$, and (ii) the t-dependent counterpart of Equation (77) for the ((31), 2) choice, with corresponding constants $D_{31,c}$ and $D_{2(31),c}$.

Assumption—We shall complete the approximation by assuming that:

$$\begin{array}{c}{D}_{12,c}=D_{3(12),c}=D_{23,c}=D_{1(23),c}=D_{31,c}=D_{2(31),c}\hfill \end{array}$$

(80)

Then, the following results follow: (1) by using Equations (A21) and (A22), it follows that the above t-dependent Equation (77) for the ((12), 3) choice transforms exactly into the corresponding equation for the ((23), 1) choice and, similarly; (2) the t-dependent Equation (77) for the ((12), 3) choice transforms exactly into the corresponding equation for the ((31), 2) choice.

It is natural to conjecture that the t-dependent quantum Equations (74)–(76) for the ((12), 3) choice transform exactly into the corresponding (also t-dependent) quantum equations for the ((23), 1) and ((31), 2) choices, by using assumptions regarding the quantum coefficients $D_{12}$, $D_{3(12)}$ and the corresponding $D_{23}$, $D_{1(23)}$, $D_{31}$ and $D_{2(31)}$ for the other two choices. The latter assumptions would be the counterparts of the classical ones implemented through Equation (80). The derivations necessary to implement such a conjecture are not transparent and, so, they will lie out of our scope here.

The equilibrium solution of Equation (77) is proportional to $exp[-{(k_{B}T)}^{-1}(v_{12}+v_{23}+v_{31})]$. One arrives at the same equilibrium solution by using the counterparts of Equation (77) for the other two choices.

It may well be that the classical Equation (77), for the three different choices and for suitable $v_{12}+v_{23}+v_{31}$, could be useful as a zeroth-order model to estimate quantities characterizing specific transitions among the different particles, as its classical $N=2$ counterpart (69) did for $N=2$ reactions [55]. A study of this possibility for $N=3$ lies outside our scope here.

7.4. Thermal and Chemical Equilibria: An Approximate (Semi-Quantitative) Discussion

Notice that no chemical potentials have been introduced, which does not appear to be straightforward in our Wigner function framework. The semi-quantitative discussion below will try to provide some bypass to discuss chemical equilibrium in a semi-quantitative way.

Let us now return to Equation (73). Upon integration we get (${\mathit{y}}_{12}=q_{eq,12}^{-1}{\mathit{q}}_{12}$, ${\mathit{y}}_{3(12)}={\mathit{q}}_{eq,3(12)}^{-1}{\mathit{q}}_{3(12)}$)

$$\begin{array}{c}\int d^3{\mathit{x}}_{12}{d}^3{\mathit{x}}_{3(12)}{d}^3{\mathit{y}}_{12}{d}^3{\mathit{y}}_{3(12)}{W}_{3,re,eq}\left({\mathit{x}}_{12},{\mathit{x}}_{3(12)},{\mathit{q}}_{12},{\mathit{q}}_{3(12)}\right)=\hfill \\ {(\frac{1}{q_{eq,12}{q}_{eq,3(12)}})}^3\sum _{j}exp[-\beta E_j]\int d^3{\mathit{x}}_{12}{d}^3{\mathit{x}}_{3(12)}\mid {\phi }_j({\mathit{x}}_{12},{\mathit{x}}_{3(12)}){\mid }^2\hfill \end{array}$$

(81)

Notice that $\int d^3{\mathit{x}}_{12}$ and $\int d^3{\mathit{x}}_{3(12)}$ are extended over the volume $\mathsf{\Omega }$ (microscopically large but macroscopically small).

For sufficiently low $k_{B}T$, the contribution of the CS (i) is the dominant one to Equation (81). As $k_{B}T$, varying initially in the range where only the system formed by particle 3 interacting with the bound state of (1, 2) is possible (lower part of the CS (i)), increases and approaches the range where the system formed by particle 1 interacting with the bound state of (2, 3) is also possible (CS (ii)), the transitions (chemical reactions) $(12)+3\to (23)+1$ and its inverse $(23)+1\to (12)+3$ occur. As such $k_{B}T$’s, other chemical reactions involving the states corresponding to CS (iii) and CS (iv) are not produced. However, transitions involving the states corresponding to CS (iii) and CS (iv) are indeed possible as $k_{B}T$ continues to increase.

We shall focus on $k_{B}T$ś such that only $(12)+3\to (23)+1$ and its inverse $(23)+1\to (12)+3$ occur (1 and 2 being bound to each other for (12) and 2 and 3 being bound to each other for (23)). Then, the dominant contributions to (81) come from such states in CS (i) and CS (ii) (the contributions from CS (iii) and CS (iv) being exponentially damped compared to the former).

The above two sets of states $(12)+3$ and $(23)+1$ in the CS (i) and CS (ii) correspond to thermal equilibrium, but there may be, in general, an imbalance between their contributions to $W_{3,re,eq}$, unless some further requirement be fulfilled. To grasp the latter, in a semi-quantitative framework, we advance that there exist key ranges of intermediate Tś and volumes (microscopically large but macroscopically small) such that the orders of magnitude of the contributions of both sets $(12)+3$ and $(23)+1$ are not significantly different from each other at thermal equilibrium. The order of magnitude of the contribution of the states $(12)+3$ to $W_{3,re,eq}$ in Equation (81) will be estimated, as a zeroth-order approximation, by approximating it as the product of the bound state wave function for the pair (12) (${\phi }_{d(12)}({\mathit{x}}_{12}$)) times the classical equilibrium distribution for particle 3 treated as a free particle with any momentum ${\mathit{q}}_{3(12)}$ (thereby omitting the interaction of the bound state of (1, 2) with 3). Accordingly, we must integrate over ${\mathit{x}}_{12}$, ${\mathit{x}}_{3(12)}$ and ${\mathit{q}}_{3(12)}$. We get:

$$\begin{array}{c}{(\frac{1}{q_{eq,12}{q}_{eq,3(12)}})}^{3}exp\left(-\beta E_{d(12)}\right)\int d^3{\mathit{x}}_{12}{\left|{\phi }_{d(12)}({\mathit{x}}_{12})\right|}^2\int d^3{\mathit{x}}_{3(12)}\times \hfill \\ \int d^3{\mathit{q}}_{3(12)}\frac{1}{{(2\pi )}^3}exp[-\beta \frac{{\mathit{q}}_{3(12)}^2}{2m_{3(12)}}]\simeq {(\frac{1}{q_{eq,12}{q}_{eq,3(12)}})}^{3}exp\left(-\beta E_{d(12)}\right)\frac{{\mathsf{\Omega }}_3{q}_{eq,3(12)}^3}{8{\pi }^{3/2}{\hslash }^3}.\hfill \end{array}$$

(82)

The effective integration region for $x_{12}$ should be smaller than ${\mathsf{\Omega }}_3$. Using the equivalent counterpart of (73) for the choice ((23), 1), a similar estimate for the bound state of (2, 3) times the classical equilibrium distribution for particle 1, treated as free, yields:

$$\begin{array}{c}\hfill {(\frac{1}{q_{eq,23}{q}_{eq,1(23)}})}^{3}exp\left(-\beta E_{d(23)}\right)\frac{{\mathsf{\Omega }}_1{q}_{eq,1(23)}^3}{8{\pi }^{3/2}{\hslash }^3}.\end{array}$$

(83)

For chemical equilibrium in the transition under study, it seems natural that there be a (at least semi-quantitative) similarity of the orders of magnitude of both Equations (82) and (83). That could be possible, for the given $E_{d(12)}$ and $E_{d(23)}$ and masses involved, for suitable ${\mathsf{\Omega }}_1$, ${\mathsf{\Omega }}_3$, for the T’s under consideration.

We shall now discuss the size of the microscopically large $\mathsf{\Omega }$ and, hence, that of ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_3$. One could argue physically that our three-particle system is contained in some macroscopic gaseous system at T, containing a very large number N of similarly interacting triplets (upon neglecting the interactions among different or disjoint triplets). That is, our three-particle system is just one of those triplets. The gaseous system occupies a macroscopically large volume ${\mathsf{\Omega }}^{\prime }\gg \mathsf{\Omega }$. Then, we argue that $\mathsf{\Omega }\simeq {\mathsf{\Omega }}^{\prime }/N$ and that the microscopically large ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_3$ could be possibly somewhat smaller than $\mathsf{\Omega }$. Notice that ${\mathsf{\Omega }}_3$ could be regarded as a (microscopically large but macroscopically small) allowed volume for particle 3, per triplet ((12), 3) when (12) are in a bound state in the gaseous system and so on for ${\mathsf{\Omega }}_1$ in a corresponding triplet ((23), 1). Such ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_3$ could allow for the compatibility between the orders of magnitude in Equations (82) and (83), for the T’s considered. For a comparison, recall that typical diatomic gases at adequately low pressures and T’s not smaller than room temperature, one has that $\mathsf{\Omega }$ is about ${10}^4$–${10}^5$ Å${}^3$ per pair [4].

Then, when such a (at least semi-quantitative) similarity of the orders of magnitude occurs, there could be simultaneous coexistence of thermal and chemical equilibria for $(12)+3\to (23)+1$ and its inverse $(23)+1\to (12)+3$, in the actual simple framework.

8. Conclusions and Discussion

In short, and to avoid repetitions, the program in this work is the following in all cases. By assumption, the system evolves subject to a $hb$, at equilibrium at absolute temperature T, in all cases, except in Section 3 (in which a large closed system will be treated). Ab initio dissipation is excluded, by assumption, in all cases. We omit any analysis of the states of the $hb$ and, in an effective way, we focus exclusively on the states of the system considered. The analysis is based upon the non-equilibrium classical Liouville probability distribution function $W_c$ and the quantum Wigner function W. Extensive use is made of the corresponding equilibrium functions $W_{c,eq}$ and $W_{eq}$ as weight functions to generate orthogonal polynomials. In turn, the latter are employed to construct non-equilibrium moments of the corresponding non-equilibrium functions, which fulfill non-equilibrium hierarchies (containing coefficients determined by the corresponding equilibrium functions) and with suitable initial conditions. Various approximate non-equilibrium classical and quantum-mechanical Smoluchowski-like equations are obtained for the lowest non-equilibrium moment for long time: they yield evolutions towards thermal equilibrium. Below, we shall emphasize distinct features in each case.

1. Reviewed material.

For the sake of simpler presentations, the various approximations leading to various Smoluchowski-like equations for lowest non-equilibrium moments are outlined for one-dimensional cases in Section 2 and Section 4, in the classical and quantum domains, respectively.

Classical case.

Section 2: classical one-dimensional particle. Here, $W_{c,eq}$ is Gaussian and the orthogonal polynomials are the standard Hermite ones. They are used to construct moments of the non-equilibrium $W_c$. In turn, those moments fulfill non-equilibrium hierarchies, the solutions of which depend on certain operator-continued fractions and combinations thereof. Under suitable long-term approximations performed on those operator-continued fractions, the lowest moment fulfills an irreversible Smoluchowski-like equation.

Section 3: classical closed N-particle system, in three dimensions, without $hb$. A specific analysis is devoted to the initial non-equilibrium state (describing equilibrium at large distances but off-equilibrium at finite distances). The analysis in Section 2 is suitably generalized. The non-equilibrium hierarchy is genuinely different from the non-equilibrium classical BBGKY hierarchy [7,8]. As an extreme approximation, an irreversible Smoluchowski-like equation (resembling that in the standard Rouse model for polymer dynamics [38]) is obtained. A non-increasing Liapunov function characterizes the long-term evolution towards equilibrium.

Quantum case.

Section 4: one-dimensional particle. Difficulties related to the non-positivity of W are bypassed, by invoking a suitable extension of the theory of orthogonal polynomials [54]. Based upon the latter, the equilibrium Wigner function is used to generate orthogonal polynomials which, in turn, lead to non-equilibrium moments. The corresponding development in Section 2 regarding evolutions is extended for the lowest non-equilibrium moment.

Section 6: various general developments in Section 5 (containing new material, as discussed below) are treated specifically for $N=2$ particles. On the other hand, developments in Section 4 are generalized non-trivially to this case. The resulting construction is applied to binary chemical reactions in the quantum domain, with applications for adequately small thermal wavelength. The corresponding classical Smoluchowski-like equation coincides with Equation (7.30) in [55], which has been a useful zeroth-order model for binary chemical reactions.

2. New material.

Classical case.

Section 2.4: an approximation for an operator-continued fraction in non-equilibrium classical statistical mechanics is proposed, in outline.

Quantum case.

Section 5: general formulation for $N>1$ quantum particles. It is a general extension, for N particles, of the case of binary chemical reactions, previously published in [34]. The statistical evolutions of the center of mass and of the relative degrees of freedom are factored out. An approximate irreversible quantum (Smoluchowski-like) equation for the lowest non-equilibrium moment describing the relative degrees of freedom is derived from the Wigner functions W.

Section 7: various general developments in Section 5 are treated specifically for $N=3$ particles in three spatial dimensions. Non-equilibrium chemical reactions involving three particles are analyzed, in typical regimes of short thermal wavelength and long times. Special attention is paid to the different possible choices of variables for the relative degrees of freedom. Based upon the developments in Section 5, an approximate irreversible quantum (Smoluchowski-like) equation for the lowest non-equilibrium moment is obtained. A semi-quantitative analysis of chemical equilibrium involving transitions among three particles is outlined. The Smoluchowski-like equations for chemical reactions involving three particles are analyzed in the classical domain, thereby justifying paying due attention, for the sake of consistency, to the different possible choices of variables for the relative degrees of freedom.

Some open problems deserving further research—We may quote: (op1) the classical three-term hierarchies and the associated operator-continued fractions, (op2) the quantum hierarchies (no longer three-term ones) and the associated operator solutions, (op3) the various approximations regarding (op1) and (op2) carried out throughout this work, (op4) chemical reactions in the quantum domain involving two and, specifically, three particles. In particular, it is suggested that the classical Equation (77) could be useful as a zeroth-order model for specific transitions among three reacting particles, as its classical $N=2$ counterpart in (69) did for $N=2$.